Here is a naive implementation of a kernel density regression algorithm in Julia.
The code works, but there is quite a bit of room for improvement though. First, there is not an optimal bandwidth selection implemented and second; it's really slow. Nevertheless, it took me a while to implement and I am really looking forward to any suggestions on how to improve the precision and performance of the code.
function kernelRegression(x, y, h) #length of the actual data n = length( x ) #if length of the data is smaller then 100, #use 100 for the length of the sequence to make sure the curve #can be plotted smoothly l = min( n, 100 ) seq = linspace(minimum(x), maximum(x), l) #bandwidth for the gaussioan kernel approx to: h = 1.06 * std( x ) * n ^(-.2) #store the result in est = zeros(Real, l) #gaussian kernel function kernel( u ) exp ( -.5 * u^2 ) end #compute the estimate over the entire sequence of values #basic algorithm: #for every x in sequence from min ( data ) to max ( data ) # for every observation in data # calulate the kernel density function at this particular point # which is the sum of all the kernels on the x axis devided by n * binwidth # multiply the kernel at this point with the respective y value # and devide by the sum of kernels for (index, sequence) in enumerate( seq ) #reset the weights weight = 0.0 for ( i, ob ) in enumerate( x ) temp = 0.0 #inner loop: calculate the mean #of the kernels at X_i = x for obs in x u = ( sequence - obs ) / h temp += kernel ( u ) end temp /= ( n * h ) #use temp to calculate the estimate of y at X_i = x v = ( sequence - ob ) / h weight += ( kernel ( v ) / temp ) * y[i] end est[index] = weight / ( n * h ) end est #scatter(x, y) #plot(seq, est, color="red") end (() -> begin x = rand(Uniform(-6, 6), 1000) y = cos(x) + rand(Normal(0, 0.2), 1000) kernelRegression(x, y) scatter(x,y) plot(linspace(-6,6,1000), y_est) end)()