parzen.m

Parzen窗法

% Parzen窗法，使用教材中的例子
function Parzen()
    clear;
     
    % 教材中a图，均值为零、方差为1的一维正态分布
    x = linspace(-2, 2, 100);
    X = norm_a(256);
    p = Parzen_norm(x, X, 0.25, 1);     subplot(4,6,1);plot(x, p);title('N =1, h_{1} = 0.25');axis([-2, 2, 0, 1]);
    p = Parzen_norm(x, X, 1, 1);        subplot(4,6,2);plot(x, p);title('N =1, h_{1} = 1');axis([-2, 2, 0, 1]);
    p = Parzen_norm(x, X, 4, 1);        subplot(4,6,3);plot(x, p);title('N =1, h_{1} = 4');axis([-2, 2, 0, 1]);
    p = Parzen_norm(x, X, 0.25, 16);    subplot(4,6,7);plot(x, p);title('N =16, h_{1} = 0.25');axis([-2, 2, 0, 1]);
    p = Parzen_norm(x, X, 1, 16);       subplot(4,6,8);plot(x, p);title('N =16, h_{1} = 1');axis([-2, 2, 0, 1]);
    p = Parzen_norm(x, X, 4, 16);       subplot(4,6,9);plot(x, p);title('N =16, h_{1} = 4');axis([-2, 2, 0, 1]);
    p = Parzen_norm(x, X, 0.25, 256);   subplot(4,6,13);plot(x, p);title('N =256, h_{1} = 0.25');axis([-2, 2, 0, 1]);
    p = Parzen_norm(x, X, 1, 256);      subplot(4,6,14);plot(x, p);title('N =256, h_{1} = 1');axis([-2, 2, 0, 1]);
    p = Parzen_norm(x, X, 4, 256);      subplot(4,6,15);plot(x, p);title('N =256, h_{1} = 4');axis([-2, 2, 0, 1]);
    p = normpdf(x, 0, 1);               
    subplot(4,6,19);plot(x, p);axis([-2, 2, 0, 1]);
    subplot(4,6,20);plot(x, p);title('均值为零、方差为1的一维正态分布');axis([-2, 2, 0, 1]);
    subplot(4,6,21);plot(x, p);axis([-2, 2, 0, 1]);
    
    % 教材中b图，均匀分布双峰概密
    x = linspace(-3, 3, 100);
    X = equal_b(256);
    p = Parzen_norm(x, X, 0.25, 1);     subplot(4,6,4);plot(x, p);title('N =1, h_{1} = 0.25');axis([-3, 3, 0, 1.5]);
    p = Parzen_norm(x, X, 1, 1);        subplot(4,6,5);plot(x, p);title('N =1, h_{1} = 1');axis([-3, 3, 0, 1.5]);
    p = Parzen_norm(x, X, 4, 1);        subplot(4,6,6);plot(x, p);title('N =1, h_{1} = 4');axis([-3, 3, 0, 1.5]);
    p = Parzen_norm(x, X, 0.25, 16);    subplot(4,6,10);plot(x, p);title('N =16, h_{1} = 0.25');axis([-3, 3, 0, 1.5]);
    p = Parzen_norm(x, X, 1, 16);       subplot(4,6,11);plot(x, p);title('N =16, h_{1} = 1');axis([-3, 3, 0, 1.5]);
    p = Parzen_norm(x, X, 4, 16);       subplot(4,6,12);plot(x, p);title('N =16, h_{1} = 4');axis([-3, 3, 0, 1.5]);
    p = Parzen_norm(x, X, 0.25, 256);   subplot(4,6,16);plot(x, p);title('N =256, h_{1} = 0.25');axis([-3, 3, 0, 1.5]);
    p = Parzen_norm(x, X, 1, 256);      subplot(4,6,17);plot(x, p);title('N =256, h_{1} = 1');axis([-3, 3, 0, 1.5]);
    p = Parzen_norm(x, X, 4, 256);      subplot(4,6,18);plot(x, p);title('N =256, h_{1} = 4');axis([-3, 3, 0, 1.5]);
    p = double_equal(x);               
    subplot(4,6,22);plot(x, p);axis([-3, 3, 0, 1.5]);
    subplot(4,6,23);plot(x, p);title('均匀分布双峰概密');axis([-3, 3, 0, 1.5]);
    subplot(4,6,24);plot(x, p);axis([-3, 3, 0, 1.5]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   function name:  norm_a
%   description:    产生教材中a图例子的样本，即均值为零、方差为1的一维正态分布
%   input variable: N(样本个数)
%   output variable:X(样本矩阵，N by 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = norm_a(N)
    X = normrnd(0, 1, N, 1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   function name:  equal_b
%   description:    产生教材中b图例子的样本，即均匀分布双峰概密
%   input variable: N(样本个数)
%   output variable:X(样本矩阵，N by 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = equal_b(N)
    f = rand(N, 1);
    X = zeros(N, 1);
    for i = 1:N
        if f(i)<0.5
            X(i) = f(i) - 2.5;
        else
            X(i) = 4 * f(i) - 2;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   function name:  double_equal
%   description:    教材中b图例子的的概率密度曲线，即均匀分布双峰概密
%   input variable: x(横坐标)
%   output variable:p(概率密度，纵坐标)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = double_equal(x)
    for i = 1:100
        if ( -2.5<x(i)) & ( x(i)<-2 )
            p(i) = 1;
        elseif ( 0<x(i)) & ( x(i)<2 )
            p(i) = 0.25;            
        else
            p(i) = 0;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   function name:  Parzen_norm
%   description:    用正态窗函数作Parzen窗估计
%   input variable: x(横坐标), X(样本), h1, N
%   output variable:p(估计，纵坐标)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = Parzen_norm(x, X, h1, N)
    hN = h1 / sqrt(N);
    sum = zeros(1, 100);
    for i = 1:N
        sum = sum + normpdf((x - X(i))/hN, 0, 1);
    end
    p = sum/(N * hN);
end