5
\$\begingroup\$

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)()
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.