function [x_g, smooth_f_x, mean_Y, std_Y] = kde(X,Y, sigma_factor)

if ~exist('Y','var')
    do_Y = false;
else
    do_Y = true;
end

if ~exist('sigma_factor','var')
    sigma_factor = 1;
end

if do_Y
    assert(size(X,1)==size(Y,1));
    n_Y = size(Y,2);
end

N = length(X);
sigma_pop_x = std(X);
sigma_x = sigma_factor*(4/(3*N))^(1/5)*sigma_pop_x;

x_min = min(X);
x_max = max(X);

n_x = 50;

x_g = linspace(x_min, x_max,n_x);

%%
if do_Y
    mean_Y = zeros(n_Y, n_x);
    meansq_Y = zeros(n_Y, n_x);
end

smooth_f_x = zeros(1,n_x);
delta_x = x_g(2)- x_g(1);

for k=1:N
    current_smooth_f_x = normpdf(x_g,X(k),sigma_x).*delta_x;
    smooth_f_x = smooth_f_x + current_smooth_f_x;

    if do_Y
        mean_Y = mean_Y + (Y(k,:)')*current_smooth_f_x;
        meansq_Y = meansq_Y + ((Y(k,:).^2)')*current_smooth_f_x;
    end
end

if do_Y
    mean_Y = multiprod(mean_Y,1./smooth_f_x,1);
    meansq_Y = multiprod(meansq_Y,1./smooth_f_x,1);
    var_Y = meansq_Y - mean_Y.^2;
    std_Y = sqrt(var_Y);
else
    mean_Y = [];
    std_Y = [];

end
%%
%close all
% no_samps = smooth_f_x < 2;
% mean_Y(:,no_samps) = NaN;
% std_Y(:,no_samps) = NaN;

%figure
%errorbar(x_g, mean_Y, std_Y)
%hold all
%scatter(X, log(Y),'.','r','SizeData',5)

end
