8/3 Parzen 窗 Matlab实现

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