1. 一维高斯函数密度函数估计
close all;clear all;
X = normrnd(0,1,1,10000);%从正态分布中产生10000个均值为0,方差为1 的样本
f = -3:0.01:3;%确定横坐标范围
function p = Parzen(X,h_1,N,x)
% X 是所有的样本
% h 是Parzen窗的窗口大小
% N 是采样的总样本大小
% x 是密度估计的点
% 采用高斯窗口大小
p = zeros(length(x),1);
h = (h_1/sqrt(N)); %半径
for i = 1:length(x)
b = 0;
for j = 1:N
u = (x(i) - X(j))/h;
b = b + exp(-u.^2/2)/(sqrt(2*pi)*h); %一维高斯分布
end
p(i) = b/N;
end
end
index = 1;
for N = [10 100 1000]
for h_1 = [0.25,1,4]
p = Parzen(X,h_1,N,f);
subplot(3,3,index); %将当前图窗划分为3 3 网格,在index指定位置创建坐标区
plot(f,p);
plot(X,Y) 创建 Y 中数据对 X 中对应值的二维线图。
如果 X 和 Y 都是向量,则它们的长度必须相同。plot 函数绘制 Y 对 X 的图。
如果 X 和 Y 均为矩阵,则它们的大小必须相同。plot 函数绘制 Y 的列对 X 的列的图。
如果 X 或 Y 中的一个是向量而另一个是矩阵,则矩阵的各维中必须有一维与向量的长度相等。如果矩阵的行数等于向量长度,则 plot 函数绘制矩阵中的每一列对向量的图。如果矩阵的列数等于向量长度,则该函数绘制矩阵中的每一行对向量的图。如果矩阵为方阵,则该函数绘制每一列对向量的图。
如果 X 或 Y 之一为标量,而另一个为标量或向量,则 plot 函数会绘制离散点。但是,要查看这些点,您必须指定标记符号,例如 plot(X,Y,'o')。
hold on;
legend(['h = ', num2str(h_1), ', N = ',num2str(N)]);
plot(f,normpdf(f,0,1),'r-');
index = index + 1;
end
end
2. 二维高斯函数密度函数估计
类似地,我们给出二维高斯密度函数估计的实现过程。
mu = [0 0]; % mu 是需要生成的数据的均值
SIGMA = [1 1]; % SIGMA 是需要生成的数据的自相关矩阵(相关系数矩阵)
% mvnrnd 是用来生成多维正态数据的
X = mvnrnd(mu,SIGMA,100000);
[fx,fy]=meshgrid(linspace(-3,3,100)',linspace(-3,3,100)');
f = [fx(:) fy(:)];
for N = [100 1000]
figure;
index = 1;
for h_1 = [0.01,1,4]
p = Parzen2d(x,h_1,N,f);
subplot(2,3,index);
surf(fx,fy,reshape(p,100,100)); hold on;
index = index+1;
end
end
function p = Parzen2d(X,h_1,N,x)
% X 是所有的样本
% h 是Parzen窗的窗口大小
% N 是采样的总样本大小
% x 是密度估计的点
% 采用高斯窗口大小
[n,d] = size(X);
p = zeros(length(x),1);
h = (h_1/sqrt(N))^(1/d);
for i = 1:length(x)
b = 0;
for j = 1:N
ux = (x(i,1) - X(j,1))/h;
uy = (x(i,2) - X(j,2))/h;
u = [ux,uy];
b = b + exp(-norm(u)^2/2)/(sqrt(2*pi)*h);
end
p(i) = b/N;
end
end