Parzen窗

Halcom

非参数检验方法，考虑有N个数据点，x属于R ，i=1,2,……，N ，采用Parzen窗来估计未知的密度函数，具体的表达式如下：

             p(x) = sum(fai( (x-xi)/h )) / N/ h

其中N足够大，h充分小，通常由用户自己设定。fai( x ) 函数为核密度估计函数。通常核密度估计函数采用高斯分布函数，因此有：

p(x) = 1/N * sum( exp(-(x-xi)^2/2/h/h )/h/sqrt(2*pi) )

MATLAB代码如下：

1. clc,clear,close all                            % 清屏、清工作区、关闭窗口
   
2. warning off                                 % 消除警告
   
3. feature jit off                               % 加速代码执行
   
4. % 概率密度函数实际是一个混合高斯分布函数
   
5. % 采用 generate_gauss_classes函数产生所需要的数据
   
6. m=[1; 4]';   % 初始化
   
7. S(:,:,1)=[0.3];
   
8. S(:,:,2)=[0.3];
   
9. P=[2/3 1/3];
   
10. N=1000;
    
11. randn('seed',0);
    
12. [X]=generate_gauss_classes(m,S,P,N);
    
13. 
    
14. % 绘图pdf
    
15. x=-5:0.1:5;
    
16. pdfx=(2/3)*(1/sqrt(2*pi*0.2))*exp(-.5*((x-1).^2)/0.2)+(1/3)*(1/sqrt(2*pi*0.2))*exp(-.5*((x-4).^2)/0.2);
    
17. plot(x,pdfx); hold on;
    
18. 
    
19. %Parzon窗计算，h = 0.1 and x in [-5, 5]
    
20. h=0.1;
    
21. pdfx_approx=Parzen_gauss_kernel(X,h,-5,5);
    
22. plot(-5:h:5,pdfx_approx,'r');
    
23. legend('原始分布函数','Parzen窗逼近效果')
    

复制代码

相应的函数如下：

1. function [px]=Parzen_gauss_kernel(X,h,xleftlimit,xrightlimit)
   
2. % 函数调用格式
   
3. %   [px]=Parzen_gauss_kernel(X,h,xleftlimit,xrightlimit)
   
4. % Parzen 使用高斯基逼近一维PDF
   
5. %输入:
   
6. %   X:           数据点
   
7. %   h:           步长
   
8. %   xleftlimit:     x的最小估计值
   
9. %   xrightlimit:    x的最大估计值
   
10. %输出:
    
11. %   px:          p(x)的估计值
    
12. 
    
13. [l,N]=size(X);
    
14. xstep=h;
    
15. k=1;
    
16. x=xleftlimit;
    
17. while x<xrightlimit+xstep/2
    
18.     px(k)=0;
    
19.     for i=1:N
    
20.         xi=X(:,i);
    
21.         px(k)=px(k)+exp(-(x-xi)'*(x-xi)/(2*h^2));
    
22.     end
    
23.     px(k)=px(k)*(1/N)*(1/(((2*pi)^(l/2))*(h^l)));
    
24.     k=k+1;
    
25.     x=x+xstep;
    
26. end

复制代码

产生服从高斯分布的数据集

1. function [X,y]=generate_gauss_classes(m,S,P,N)
   
2. % 产生服从高斯PDF分布的数据集
   
3. %输入:
   
4. %   m:  均值向量
   
5. %   S:  正态分布的协方差矩阵
   
6. %   P:  先验概率密度
   
7. %   N:  总数据个数
   
8. %输出:
   
9. %   X:  待分类的数据
   
10. %   y:  分类的标签，对应于X
    
11. [l,c]=size(m);
    
12. X=[];  % 初始化
    
13. y=[];
    
14. for j=1:c
    
15. % 产生每一个分布下的 the [p(j)*N] 向量
    
16.    t=mvnrnd(m(:,j),S(:,:,j),fix(P(j)*N))';
    
17.    % 由于采用的是fix函数，总数据数可能小于N
    
18.    X=[X t];
    
19.    y=[y ones(1,fix(P(j)*N))*j];
    
20. end

复制代码