非参数估计
模式识别上机实验3:密度的非参数估计 学号：[1**********] 姓名：寸正雄
给出密度函数的非参数估计公式，并产生1、16、256和16384个服从一维标准正态分布的样本， 1． 分别就窗宽为h 1=0. 25, 1, 2． 分别就k N =1题结果：
4，h N =h 1/
窗函数为高斯函数的情形估计所给样本的密度函数并划出图形。 N ，
N 时用k N 近邻方法估计所给样本的密度函数并划出图形。
Parzen 窗法对个样本的估计图像表
2题结果：
k N
邻近法对各样本的估计图像表
\
模式识别上机实验3:密度的非参数估计 学号：[1**********] 姓名：寸正雄
1. 问题分析
该实验要求利用总体分布的非参数估计对利用随机产生的不同大小的样本进行密度函数估计。利用总体分布的
非参数估计中的两种常用方法-Parzen 窗法和k N 邻近估计，并讨论Parzen 窗法中不同窗宽h i 的估计效果。
1.1. 总体分布的非参数估计
总体分布的分参数估计就是用某些方法直接利用样本来估计总体分布。估计总体分布函数的方法有许多种，但他们的基本思想都是采用某钟办法证明估计的收敛性，一个最根本的出发点是随机向量x 落入到区域R 的概率P 为
P =
⎰
R
p (x ) dx （1.1）
p (x ) 为x 的总体概率密度函数。如有N 个样本x 1, x 2, x N 是从密度为p (x ) 的总体中独立抽取的，则N
个样本中有k 个落入区域R 中符合随机变量的二项分布，更具二项分布的性质可以得到p (x ) 的估计为
∧
p (x ) =
k /N V
（1.2）
其中N 表示赝本数，k 表示落入区域R 的样本个数，V 表示区域R 的度量。为了估计x 点的密度，构造一串包括x 的区域系列R 1, R 2, , R N ，对各个区域进行样本估计
∧
p i (x ) =
∧
k i /N V i
, i =1, 2, , N （1.3）
各区域的估计函数p i (x ) 知道后，总体的密度函数估计值为
∧
N
p (x ) =
∑
i =1
∧
p i (x ) =
1N
N
∑V
i =1
k i
i
（1.4）
∧
其中V i 表示区域R i 的度量，k i 表示落在区域R i 的样本个数。p i (x ) 收敛的必要条件是
lim V N =0，lim k N =0，lim
N →∞
k N N
N →∞N →∞
=0 （1.5）
1.2. Parzen 窗法
1.2.1．Parzen 窗估计
固定（1.3）式中的V i ，就是构造区域度量相等的区域系列。定义一个窗函数ϕ(u )
⎧1, u ≤2
（1.6） ϕ(u )=⎨
⎩0, 其他
ϕ(u )的意思就是：若样本落在以u 为中心，宽为1的区域中值为1，否则值为0。可求出落在区域R i 中的样本个数为
N
k N =
∑ϕ
i =1
⎛x -x i ⎫
⎪ (1.7) ⎪
⎝V x ⎭
把（1.7）式带入到（1.4）式中可得到parzen 窗法的估计的基本公式
∧
p (x ) =
∧
1N
N
∑V
i =1
1
N
ϕ
⎛x -x i ⎫
⎪ （1.8） ⎪
⎝h N ⎭
1.2.2．估计量p N (x ) 的密度函数的条件
要是p N (x ) 满足（1.5）式中的三个条件，只要窗函数ϕ(u )满足以下两个条件
∧
ϕ(u )≥0，⎰ϕ(u )du =1 （1.9）
1.2.3. 窗函数的选择
窗函数有多种选择，只要满足（1.5）式中的三个条件就可，一般有以下三个比较常见的窗函数
⎧1, u ≤1
2
方窗函数 ϕ(u )=⎨
⎩0, 其他
正太窗函数 ϕ(u )=指数窗函数 ϕ(u )=e 它们的图像如下：
12π
-u
e
12-u 2
图2.1
方窗函数 正太窗函数 指数窗函数
1
0.9
0.8
0.3
0.70.60.50.40.3
0.1
0.20.1
0.05
0.20.1
0.25
0.20.15
0.70.60.50.40.3
0.4
0.35
1
0.90.8
1.3. k N 邻近估计
固定（1.3）式中的k i ，就是构造落在某区域的样本数相同，而区域度量不相等的区域系列，这样可得到k N 邻近估计在x 处的密度估计值为
ˆ(x )=p
k /N V i
（1.10）
其中
N =16384k =128
是样本总数，V i 是点x 的影响范围，定义如下
设X i =⎡X i , X i ,... X i ⎤是样本中距离x 点最近的k 的样本点，取R 为这k 个点中距离x 最远的一个点到x
⎣⎦
1
2
k
点的距离，即R =max (X i -x )，V i 就是以x 为中心以R 为半径的区域的度量。
2. 问题解答 2.1．随即产生样本
利用matlab 中的 normrnd()函数可以生成服从正太分布的样本，其格式如下
X =normrnd
(μ, σ, d 1, d 2) （2.1）
其中μ是均值，σ是方差，d 1, d 2表示要产生的样本X 的维数。
按题意要产生1、16、256和16384个服从一维标准正态分布的样本，（2.1）式中取μ=0，σ=1，d 1=1，d 2分别取1、16、256和16384就可以生成随机样本X
d 1⨯d 2
。
2.2．Parzen 窗法
提意要求窗函数选择高斯函数，即正太函数
ϕ(u )=
12π
e
12-u 2
（2.2）
再利用Parzen 窗法的基本公式（1.8）
∧
p N (x ) =
1N
N
∑V
i =1
1
N
ϕ
⎛x -x i ⎫
⎪ （2.3. ） ⎪
⎝h N ⎭
上式中N 分别取1、16、256和16384，h 1分别取0.25,1和4，V N =h N =h 1/
N 可得到结果如下表：
表2.1 Parzen窗法对个样本的估计图像表
2.3．k N 邻近估计
题意要求分别就k N =式（1.10），取k =
N 时用k N 近邻方法估计所给样本的密度函数并划出图形。根据k N 邻近估计的基本公
ˆ(x )=p
/N V i
=
(2.4)
这里是的样本是一维的，取V i =2R =2max (X i -x )可得到其结果如表2.2中。
题中所给出的样本为正态分布，而在k N 临近估计中，当样本数较少时得到的并不像正太图像，同样采用k N 临近的思想，就取相邻的k 个样本进行估计
ˆi (x )=p
k /N V i
, i =1, 2,..., K （2.5）
总体的密度估计公式为
K
ˆ(x )=p
∑
i =1
ˆi (x )=p
1N
∧
K
∑V
i =1
k
i
（2.6）
而在公式中还有一个问题，V i 要如何确定才能使得p (x ) 满足（1.5）式的三个要求？与Parzen 窗法类似，利用k N 邻近法估计时要注意所要包含k 个样本的窗函数的选择，为了计算方便，常采用方窗函数。 采用方窗函数可定义如下
11⎧1, min X -≤u ≤max X +()()⎪i i
ϕ(u )=⎨22（2.7）
⎪0, other ⎩
其中X i =⎡X i , X i ,... X i ⎤是样本中相邻的k 个样本
⎣⎦
1
2
k
ˆ(x )的要求和更好的反应样本，X i 可采用逐次推进的方式，如单N =16, k =4时，为了保证满足p
X 1=⎡, X 12, X 13, X 14⎤, X 2, X 2, X ⎤=⎡(1:)4⎤2⎣X ⎦, X 2=⎡⎣11234⎣X 1⎦=⎡⎣X 2⎦
∧
（2.6）式中K 应该... 。上面（2.5）(X 2):⎤⎦5
取N -k +1才能保证p (x ) 满足要求。
这里也可以采用正态窗函数来进行总体估计。采用正态窗函数时可利用上把正态窗函数中的
μ=m e a n (X i ), i =i 1, i 2, , i
N
这样可得到采用正态窗函数的一般估计公式为：
∧
p N (x ) =
1N -K +1
N -K +1
∑ϕ(mean
i =1
(X i ) ) （2.8）
其结果如下表
表2.2
k N
邻近法对各样本的估计图像表
3. 程序代码
实验结果如上面各表所示，其MA TLAB 代码如下： 3.1. 三个常用窗函数 chuanghanshu.m
function fu=chuanghanshu(u,t) if t==1
if abs(u)
fu=exp(-(u^2)/2)/sqrt(2*pi); elseif t==3
fu=exp(-abs(u)); end
3.2. 绘制三个常用窗函数图像huizhiCHS.m
u=-5:0.1:5;N=size(u); for i=1:N(2)
fu1(i)=chuanghanshu(u(i),1); fu2(i)=chuanghanshu(u(i),2); fu3(i)=chuanghanshu(u(i),3);
end
subplot(1,3,1),plot(u,fu1); subplot(1,3,2),plot(u,fu2); subplot(1,3,3),plot(u,fu3);
3.3.Parzen 窗法 Parzen.m
function px=Parzen(X,h1) syms x ;
N=size(X);N=N(2); hN=h1/sqrt(N); for i=1:N
pNx(i)=chuanghanshu((x-X(i))/hN,2)/hN/N; ezplot(pNx(i));hold on ; plot(X(i),0,'*');hold on ; end
px=sum(pNx);ezplot(px);hold off ;
3.4.
k N 邻近法 KNLJ.m
function KNLJ(X,h1)
x=linspace(min(X)-0.5,max(X)+0.5,100);step=(max(X)-min(X))/100; N=size(X);N=N(2); SX=sort(X);px=0; for i=1:100
absX=sort(abs(SX-x(i)));d=absX(h1); px(i)=1/(sqrt(N)*2*d); end
plot(x,px,'r' );hold off ;
3.5. 方窗k N 邻近法 KNLJ.m
function px=FCKNLJ(X,h1)
x=linspace(min(X)-0.5,max(X)+0.5,100); N=size(X);N=N(2); SX=sort(X);px=0; for i=1:N-h1+1 XI=SX(i:i+h1-1);
u=(x-mean(XI))./(max(XI)-mean(XI)+1/2); for j=1:100 if abs(u(j))
pNx(j)=1/(max(XI)-min(XI)+1)/(N-h1+1); else pNx(j)=0; end end
plot(x,pNx);hold on ;
plot(X(i),0,'*');hold on ; px=px+pNx; disp(i); end
plot(x,px,'r' );hold off ;
3.6. 正态窗k N 邻近法 KNLJ.m
function px=KNLJ(X,K) %kN邻近 syms x ;
N=size(X);N=N(2); SX=sort(X); for i=1:N-K+1
pNx(i)=chuanghanshu((x-mean(SX(i:i+K-1))),2)/(N-K+1); ezplot(pNx(i));hold on ; end
px=sum(pNx);ezplot(px);hold off ;
3.7. 执行过程如下
因为用该方法绘制出来的图像会比较复杂，所以采用每个图一个figure 的方式，就依次执行如下命令：
X1=normrnd(0,1,1,1);
Parzen(X1,0.25)执行后会出图形把图像复制到表中在执行 Parzen(X1,1)…… …………
X2=normrnd(0,1,1,16); Parzen(X2,0.25)…… …………
后面的KNLJ 同样如此 X1=normrnd(0,1,1,1); KNLJ(X1,1)… FCKNLJ （X1，1） ZTKNLJ （X1,1） ………
X2=normrnd(0,1,1,16) KNLJ(X2,4)… ……
若样本数太大ezplot 无法绘出图像可采用下方法 把ezplot 改成
x=linspace(SX(1)-1,SX(N)+1,100 ); plot(x,eval(px));hold off ; 若eval 不支持数组的可以用for 来搞定。