% Data
randn('seed',1) % Used for reproducibility
data = [randn(100,1)-10; randn(100,1)+10]; % Two Normals mixed
% True
phi = @(x) exp(-.5*x.^2)/sqrt(2*pi); % Normal Density
tpdf = @(x) phi(x+10)/2+phi(x-10)/2; % True Density
% Kernel
h = std(data)*(4/3/numel(data))^(1/5); % Bandwidth estimated by Silverman's Rule of Thumb
kernel = @(x) mean(phi((x-data)/h)/h); % Kernel Density
kpdf = @(x) arrayfun(kernel,x); % Elementwise application
% Plot
figure(2), clf, hold on
x = linspace(-25,+25,1000); % Linear Space
plot(x,tpdf(x)) % Plot True Density
plot(x,kpdf(x)) % Plot Kernel Density with Silverman's Rule of Thumb
kde(data) % Plot Kernel Density with Solve-the-Equation Bandwidth
#The same code with R language
#` Data
set.seed(1) #Used for reproducibility
data = c(rnorm(100,-10,1),rnorm(100,10,1)) #Two Normals mixed
#` True
phi = function(x) exp(-.5*x^2)/sqrt(2*pi) #Normal Density
tpdf = function(x) phi(x+10)/2+phi(x-10)/2 #True Density
#` Kernel
h = sd(data)*(4/3/length(data))^(1/5) #Bandwidth estimated by Silverman's Rule of Thumb
Kernel2 = function(x) mean(phi((x-data)/h)/h) #Kernel Density
kpdf = function(x) sapply(x,Kernel2) #Elementwise application
#` Plot
x=seq(-25,25,length=1000) #Linear Space
plot(x,tpdf(x),type="l",ylim=c(0,0.23),col="red") #Plot True Density
par(new=T)
plot(x,kpdf(x),type="l",ylim=c(0,0.23),xlab="",ylab="",axes=F) #Plot Kernel Density with Silverman's Rule of Thumb