% This code generates a single frame of the above animation. 
% 10 frames are generated and put together to create the gif animation.
function densityEstimation

N = 100;
d = 3;
pp = 0.4;  % Mixture ratio
p = rand(N, 1) > pp;

% 1) Plot the Actual density
x = -4:0.01:6.5;
y = (1-pp)*normpdf(x,0,1) + pp*normpdf(x,d,1);
plot(x, y, 'LineWidth', 2.5)

% 2) Get samples from the density
samp = genSamp(0, d, p);
hold on
stem(samp, 0.05*ones(size(samp)), 'r', 'LineWidth',1.5, 'MarkerSize', 6)

% 3) Put small Gaussian kernels around each sample
lambda = 0.35;
for i = 1:N
    x = (samp(i)-3*lambda):0.01:(samp(i)+3*lambda);
    y = normpdf(x,samp(i),lambda);
    h = plot(x, 0.03*y, ':g', 'LineWidth',1.5);
    set(h,'Color',0.5*[1,1,1])
end

% 4) Find the kernel density estimate
x = -4:0.01:6.5;
y = zeros(size(x));
for i = 1:N
    y = y + normpdf(x,samp(i),lambda);
end
y = y/N;
plot(x, y, '--k', 'LineWidth', 2.5)

% Format axes
set(gca,'LineWidth',2, 'FontSize', 20)
axis([-4 8 0 0.3])

end

function samp = genSamp(x1, x2, p)

samp = p.*(x1+ randn(length(p), 1)) + (1-p).*(x2+ randn(length(p), 1));
end