第3章
CHAPTER 3


概率密度函数的估计






贝叶斯决策前提条件是已知各类的先验概率和类条件概率，但实际中所得到的只是样本集，需要由样本集计算所需的概率密度函数，即需要进行概率密度函数的估计。先验概率的估计可以由各类样本在总体样本集中所占的比例进行估计，本章主要介绍类条件概率密度估计的方法，即已知ωj,j=1~c类若干样本xi,i=1~N，采用某种规则估计出样本所属类的概率密度函数p(x|ωj)的方法。
3.1基本概念
1. 概率密度函数估计的方法
概率密度函数的估计方法可以分为两大类： 参数估计(Parametric Estimation)和非参数估计(Nonparametric Estimation)。
如果已知类条件总体概率密度函数形式，未知其中部分或全部参数，用样本来估计这些参数，称为参数估计。如例212中，已知样本集服从正态分布，未知正态分布的均值向量和协方差矩阵，根据样本集确定μ和Σ，即参数估计。本章学习最大似然估计、最大后验估计、贝叶斯估计方法。
非参数估计指的是概率密度函数形式也未知，需要求函数本身。本章主要学习直方图、kN近邻以及核密度估计方法。
2. 参数估计的基本概念
参数估计是统计推断的基本问题之一，回顾一下涉及的几个基本概念。
(1) 统计量。假如概率密度函数的形式已知，但表征函数的参数θ未知，则可将θ的估计值构造成样本的某种函数，这种函数称为统计量。
(2) 参数空间。参数θ的取值范围称为参数空间，表示为Θ。
(3) 点估计、估计量和估计值： 构造一个样本的函数θ^，将其观察值作为未知参数θ的估计，称为点估计，θ^称为θ的估计量，由具体样本(x1,x2,…,xN)例题中使用一维数据。作为自变量计算出来的θ^(x1,x2,…,xN)值称为θ的估计值。
(4) 区间估计。通过从总体中抽取的样本，根据一定的正确度与精确度的要求，构造出适当的区间，作为未知参数的真值所在范围的估计。这个区间称为置信区间。置信区间是一种常用的区间估计方法，是以统计量的置信上限和下限为上下界构成的区间。
本章要学习的参数估计属于点估计问题，点估计有不同的方法，同一参数，用不同的估计方法求出的估计量可能不相同，需要合适的标准来评估估计量。常用的标准有无偏性、有效性、一致性等，此处不再赘述，可以参看数理统计类资料。
3.2参数估计
本节主要学习参数估计方法，包括最大似然估计(Maximum Likelihood Estimation，MLE)、最大后验(Maximum a Posteriori Probability，MAP)估计、贝叶斯估计(Bayesian Estimation)方法。





3.2.1最大似然估计
最大似然估计是通过定义似然函数，对似然函数求最大值，进而确定未知参数的估计值的方法。
1. 前提假设
最大似然估计方法的使用，需要满足下列基本假设： 
(1) 待估计参数θ是确定(非随机)而未知的量。
(2) 样本集分成c类X1,X2,…,Xc，Xj的样本是从概率密度为p(x|ωj)的总体中独立抽取出来的，即满足独立同分布条件。
(3) 类条件概率密度p(x|ωj)具有某种确定的函数形式，其参数θj未知，可以把p(x|ωj)表达为p(x|ωj,θj)或p(x|θj)，表示与θj有关。
(4) 不同类别的参数在函数上独立，可以分别对每一类单独处理。
在这些假设的前提下，参数估计的问题描述为： 分别从概率密度为p(x|ωj)的总体中独立抽取构成c个样本集X1,X2,…,Xc，利用Xj的样本估计p(x|ωj)的参数θj。
2. 似然函数
样本集X 包含N个样本，即X={x1,x2,…,xN}，各样本相互独立，获得样本集的概率l(θ)即各个样本的联合概率，定义l(θ)为样本集X 的似然函数(Likelihood Function)，如式(31)所示。

l(θ)=p(X|θ)=p{x1,x2,…,xN|θ}=∏Ni=1p(xi|θ)(31)

由于各样本相互独立抽取，N个随机变量x1,x2,…,xN的联合概率密度p{x1,x2,…,xN|θ}等于各自概率密度的乘积。
独立抽取时，样本集中的样本最有可能来源于概率密度最大的情况。似然函数定义为联合概率密度，样本独立抽取时为概率密度的乘积，所以，已知一组样本，最有可能来自似然函数最大所对应的情况。因此，可以利用似然函数作参数估计。令l(θ)为样本集的似然函数，如果θ^是参数空间中使l(θ)最大化的值，那么θ^为θ的最大似然估计量，如式（32）所示。

θ^=argmaxl(θ)(32)

至此，估计问题转化为求极值的问题。
为便于分析，可以采用对数函数将连乘转化为累加，定义对数似然函数为

h(θ)=lnl(θ)=∑Ni=1lnp(xi|θ)(33)

由于对数函数的单调性，使对数似然函数h(θ)最大的θ^也必然使得似然函数l(θ)最大。

3. 求解估计量
通过定义似然函数，将估计问题转化为求极值的问题。根据未知量θ的情况，最大似然估计的求解可以分为如下几种情况。
1) 未知参数为一元情况
只有一个未知参数，最大似然估计量是式(34)或式(35)的解。

dl(θ)dθ=0(34)

或

dh(θ)dθ=0(35)

【例31】设样本服从指数分布，密度函数为p(x)=λe-λx，x≥0

0，x＜0，求λ的最大似然估计量。
解： 设样本集X={x1,x2,…,xN}，定义似然函数为

l(θ)=∏Ni=1λe－λxi

λ是未知参数。定义对数似然函数为


h(θ)=lnl(θ)=∑Ni=1lnλe－λxi=∑Ni=1lnλ－λxi=Nlnλ－∑Ni=1λxi

因为


θ=λ


所以


dhdλ=Nλ－∑Ni=1xi=0


所以


λ^=N∑Ni=1xi=1x-


λ的最大似然估计量为1/x-，即根据样本均值的倒数。
2) 未知参数为多元情况

θ=θ1θ2…θsT是由多个未知参数组成的向量，似然函数对θ各分量分别求导，组成方程组，求最值，即解

l(θ)θ1=0,l(θ)θ2=0,…,l(θ)θs=0(36)

或

h(θ)θ1=0,h(θ)θ2=0,…,h(θ)θs=0(37)

【例32】设x服从正态分布N(μ,σ2)，其中参数μ、σ2未知，求它们的最大似然估计量。
解： 设样本集X={x1,x2,…,xN}，定义似然函数为

l(θ)=∏Ni=1p(xi|θ)

定义对数似然函数并化简为


h(θ)=lnl(θ)=∑Ni=1lnp(xi|θ)=∑Ni=1ln12πσexp－xi－μ22σ2

=∑Ni=1－12ln2π－12lnσ2－xi－μ22σ2

=－N2ln2π－N2lnσ2－∑Ni=1xi－μ22σ2


因为


θ=μ,σ2T


所以


hμ=－∑Ni=1－2xi－μ2σ2=∑Ni=1xi－μσ2=0

hσ2=－N2σ2+12σ22∑Ni=1xi－μ2=0

μ^=1N∑Ni=1xi=x-

σ^2=1N∑Ni=1xi－x-2



μ^和σ^2是μ、σ2的最大似然估计量。
3) 特殊情况
如果导数方程无解，可以根据具体数据来定。
【例33】设X={x1,x2,…,xN}是来自p(X|θ)的随机样本，当0≤x≤θ时，p(x|θ)=1/θ，否则为0。证明θ的最大似然估计是maxixi。
解： 定义似然函数为

l(θ)=∏Ni=1p(xi|θ)=∏Ni=11θ

对数似然函数为

h(θ)=lnl(θ)=－Nlnθ

令导数为零得

dhdθ=－N·1θ=0

方程的解为无穷大，但实际问题中θ≠∞。已知有N个随机样本，且0≤x≤θ时，p(x|θ)=1/θ，所以，取θ^为样本中的最大值，即θ的最大似然估计是maxixi。
三个例题分别实现了对指数分布、正态分布和均匀分布参数的最大似然估计，其他不同分布参数估计的过程类似： (1)生成数据集； (2)定义似然函数或对数似然函数； (3)对函数求最大值得参数的最大似然估计量。
4. 仿真实现
MATLAB提供了一系列的函数对不同的分布进行最大似然估计，下面分别进行简要介绍。
1) mle函数
(1) ［PHAT, PCI］=mle(DATA)： 对向量DATA中的数据采用最大似然估计方法估计正态分布参数。返回PHAT为一个向量，其元素分别为正态分布的均值和标准差； PCI为95%的置信区间。
(2) ［...］=mle(DATA,'distribution',DIST)： 对字符串DIST指定的分布进行参数估计，DIST可以是'beta'、'discrete uniform'、'exponential'、'normal'、'poisson'、'rayleigh'、'uniform'等27种取值。
(3) ［...］=mle(DATA,'pdf',PDF,'cdf',CDF,'start',START,...)： 对由概率密度函数PDF和累积分布函数CDF指定的分布的参数进行最大似然估计。PDF和CDF是使用@创建的函数句柄，以数据向量和分布参数为输入，分别返回概率密度和累积分布取值。START是参数初始值向量。
2) 分布参数估计类函数
MATLAB统计工具箱中有一系列函数，函数名以fit结尾，用来求常见分布的参数的最大似然估计和置信区间估计，例如normfit、poissfit、unifit、expfit等。这里主要介绍根据样本观测值求正态分布均值和标准差的估计函数normfit，其调用格式如下。
(1) ［MUHAT,SIGMAHAT］=normfit(X)： 估计X中数据服从的正态分布的参数，MUHAT是均值估计，SIGMAHAT是标准差估计。MUHAT等于均值参数的MLE，但SIGMAHAT是方差的无偏估计量的平方根(σ=∑ixi-μ2/N-1)，不等于标准差参数的MLE(σ=∑ixi-μ2/N)，在样本数很大的情况下区别可以忽略不计。
(2) ［MUHAT,SIGMAHAT,MUCI,SIGMACI］=normfit(X)： 同时返回95%的置信区间。
(3) ［MUHAT,SIGMAHAT,MUCI,SIGMACI］=normfit(X,ALPHA)： 同时返回(1-ALPHA)的置信区间。
3) fitdist函数
fitdist函数使用最大似然估计拟合大多数的概率分布，其调用格式如下。
(1) PD=fitdist(X,DISTNAME)： 用DISTNAME指定的概率分布来拟合列向量X中的数据，返回的PD指明分布类型及其参数。DISTNAME取'kernel'时拟合非参数核平滑分布，此外，可以取如'beta'、'binomial'、'exponential'、'logistic'、'normal'、'poisson'、'rayleigh'等23种取值，分别对应不同的分布。
(2) ［PDCA,GN,GL］=fitdist(X,DISTNAME,'BY',G)： 按照成组变量G给定的分组对X中的数据分别进行估计。返回的PDCA为元胞数组，每个元素为对应组的分布类型和参数； GN是表明组别标签的元胞数组； GL是表明组别变量级别的元胞数组，每个分组变量对应于其中一列。
(3) PD=fitdist(..., 'NAME1',VALUE1,'NAME2',VALUE2,...)： 指定参数实现估计。
4) pdf函数
pdf函数获取指定分布的密度函数，其调用格式如下。
(1) Y=pdf(NAME,X,A)： 计算概率密度函数值Y，NAME指定分布类型，有'Beta'、'Exponential'、'Normal'、'Poisson'等26种取值； A为该分布的参数，X是要计算函数值的变量点。
(2) Y=pdf(NAME,X,A,B)或Y=pdf(NAME,X,A,B,C)： 功能同上，概率分布函数有两个或三个参数，由A、B(和C)指定。
5) cdf函数
cdf函数获取指定分布的累积分布函数，其调用格式如下。
(1) Y=cdf(NAME,X,A)： 计算累积分布函数值Y，NAME、A、X同pdf中的参数。
(2) Y=cdf(NAME,X,A,B)或Y = cdf(NAME,X,A,B,C)： 概率分布函数有两个或三个参数情况下计算累积分布函数值。
(3) Y=cdf(NAME,X,A,'upper')、Y=cdf(NAME,X,A,B,'upper')或Y=cdf(NAME,X,A,B,C,'upper')： 计算累积分布函数的上尾概率。
【例34】利用normrnd函数生成正态分布数据，分别用normfit和mle函数对参数进行估计。
程序如下：

clc,clear,close all;

rng('default') 

N=100;

x=normrnd(0,1,［N,1］);                        %生成单变量的正态分布数据

［muHat,sigmaHat］=normfit(x)                 %采用normfit函数估计均值和标准差

sigmaHat_MLE=sqrt((N-1)/N)*sigmaHat       %将标准差换算为MLE估计的标准差

PHAT=mle(x)                                    %调用mle函数进行估计

运行程序，将在命令窗口输出：

muHat = 0.1231

sigmaHat = 1.1624

sigmaHat_MLE = 1.1566

PHAT = 0.12311.1566

即采用normfit函数估计数据的均值muHat为0.1231，标准差sigmaHat为1.1624； mle函数的输出PHAT为均值和标准差组成的向量［0.12311.1566］，两个函数输出的标准差关系如这条命令所示： sigmaHat_MLE=sqrt((N-1)/N)*sigmaHat。
【例35】MATLAB中的hospital数据集中有100个病人的数据，包括性别、年龄、体重、是否抽烟、血压等数据，读取该数据集，按照性别分为两类，采用fitdist函数对体重数据进行正态分布拟合，并绘制概率密度函数曲线。
程序如下：

clc,clear,close all;

load hospital                                          %导入hospital数据集

x=hospital.Weight;                                    %获取体重数据

gender=hospital.Sex;                                  %按性别分组

［pdca,gn,gl］=fitdist(x,'Normal','By',gender);    %用正态分布对两组数据进行拟合

female=pdca{1};                                        %第一组对应分布

male=pdca{2};                                          %第二组对应分布

x_values=50:250; 

femalepdf=pdf(female,x_values);                     %计算概率密度函数值

malepdf=pdf(male,x_values);

figure

plot(x_values,femalepdf,'Color','r','LineWidth',2);%绘制概率密度函数曲线

hold on

plot(x_values,malepdf,'Color','b','LineStyle','--','LineWidth',2);

xlabel('x'),ylabel('p(x)');

legend(gn,'Location','NorthEast')

hold off

程序运行结果如图31所示。从工作区窗口可以看到female组呈正态分布，估计的均值为130.4717，标准差为8.3034； male组同样为正态分布，估计的均值为180.5319，标准差为9.1932。


图31体重数据对应的正态概率密度函数曲线


3.2.2最大后验估计
在最大似然估计中，认为参数θ是确定的量，假如把θ看作随机变量，在估计中，需要考虑θ本身服从的分布。因此，利用贝叶斯公式计算θ的后验概率，最大后验概率对应的参数值为参数的估计值，这种方法即是最大后验估计。
1. 基本原理
设参数θ是先验概率为P(θ)的随机变量，其取值和样本集有关。设样本集为X={x1,x2,…,xN}，各样本相互独立，使P(θ|X)取最大值的θ^为θ的最大后验估计。通过贝叶斯公式进行P(θ|X)的计算，如式(38)所示。

P(θ|X)=P(θ)p(X|θ)p(X)=P(θ)p(X|θ)∫P(θ)p(X|θ)dθ(38)

式中，p(X|θ)是样本集X 的似然函数，如式(31)所示。
由于p(X)相对于θ独立，最大后验估计即是求解如下问题

argmaxθP(θ)l(θ)(39)

或

argmaxθlnP(θ)+lnl(θ)(310)

假设θ本身服从均匀分布，即对于所有的θ，P(θ)是常量，比较式(39)和式(32)，两者一致，即最大后验估计和最大似然估计结果一致； 一般情况下，两种估计会得到不同的结果。
【例36】设X={x1,x2,…,xN}是来自总体分布为N(μ,σ2)的样本集，已知μ服从N(μ0,σ20)分布，假设σ2已知，用最大后验估计的方法求μ的估计量μ^。
解： 确定μ的先验分布P（μ）。已知μ服从N(μ0,σ20)分布，得

P（μ）=12πσ0exp－12μ－μ0σ02

求样本联合分布p(X|μ)，即似然函数l（μ）。样本集X 总体分布为N(μ,σ2)，即p(xi|μ)~N(μ,σ2)

所以l（μ）=∏Ni=1pxi|μ=∏Ni=112πσexp－12xi－μσ2


定义准则函数J（μ）为


J（μ）=lnP（μ）+lnl（μ）

=-ln2πσ0－12μ－μ0σ02-Nln2πσ－∑Ni=112xi－μσ2


求导数并令其等于零

dJdμ=－1σ0μ－μ0σ0+∑Ni=11σxi－μσ=0

化简

μ0σ20+1σ2∑Ni=1xi－1σ20+Nσ2μ=0

得μ的最大后验估计

μ^=σ2σ2+σ20Nμ0+σ20σ2+σ20N∑Ni=1xi

2. 仿真实现
【例37】设X={x1,x2,…,xN}是来自总体分布为对数正态分布的样本集，对数正态分布的概率密度函数为p(x)=1σx2πexp－lnx－μ22σ2,(x＞0)。已知μ服从N(0,3)分布，σ=2，用最大后验估计的方法求μ的估计量μ^。
设计思路： 
设定均值μ，标准差σ=2，利用lognrnd函数生成对数正态分布样本x； 设定μ的取值范围为［-5，5］，计算准则J（μ）=lnP（μ）+lnl（μ）的取值，并取最大值对应的μ作为估计值。
程序如下：

clc,clear,close all;

rng('default') 

N=500;

mu0=0; sigma0=3;                 %μ先验分布参数

mu_value=-5:0.01:5;                   %μ的取值范围

mu=2; sd=2;                      %对数正态分布参数

training=lognrnd(mu,sd,［N,1］);      %生成对数正态分布样本500个

Poster=zeros(length(mu_value),1);

for j=1:length(mu_value)

lpdf=lognDist(training,mu_value(j),sd);     %当前均值下对数正态分布的对数似然函数

Poster(j)=lpdf+normalDist(mu_value(j),mu0,sigma0);%计算J(μ)

end

plot(mu_value,Poster);                %绘制J(μ)函数图形

xlabel('μ'),ylabel('p(μ|X)');

［~,pos］=max(Poster);                  %找最大值

mu_value(pos)

function lpdf=lognDist(X,Mu,Sigma)              %计算lnl(μ)

Z=(log(X)-Mu)./Sigma;

lpdf=sum(-log(Sigma*X)- .5*log(2*pi)- .5*(Z.^2));

end

function lpdf=normalDist(X,Mu,Sigma) %计算lnP(μ)

Z=(X-Mu)./Sigma;

lpdf=sum(-log(Sigma)- .5*log(2*pi)- .5*(Z.^2));

end


图32J（μ）函数图形




运行程序，绘制的准则函数图形如图32所示，准则函数的极大值对应μ的估计值。程序在命令窗口输出这个估计值： ans = 1.9500。

3.2.3贝叶斯估计
参数估计寻找最优的参数值，可以看作是在连续的参数空间Θ对参数的取值进行决策的问题，因此，采用决策的思路来进行估计。根据样本集X，找出一个估计量，用来估计X所属总体分布的某个真实参数，使带来的贝叶斯风险最小，这种方法称为贝叶斯估计。
1. 基本原理
样本独立抽取的样本集为X={x1,x2,…,xN}，用估计量θ^近似代替真实参数θ，带来的损失为λθ,θ^，参数θ的先验分布为P(θ)，定义Rθ^|X为给定X条件下估计量θ^的条件风险。

Rθ^|X=∫Θλ(θ,θ^)P(θ|X)dθ(311)

如果θ的估计量θ^*使条件风险Rθ^|X最小，称θ^*是关于θ的贝叶斯估计量，即

θ^*=argminθ^Rθ^|X(312)

条件风险的计算需要定义损失函数λθ,θ^，可以定义为不同的形式，最常用的是平方误差损失函数，即

λθ,θ^=θ－θ^2(313)

可以证明，如果采用平方误差损失函数，在给定样本集X 条件下，θ的贝叶斯估计量为

θ^*=E［θ|X］=∫ΘθP(θ|X)dθ(314)

整理贝叶斯估计的步骤： 
(1) 确定θ的先验分布P(θ)。
(2) 根据式(31)，由样本集X={x1,x2,…,xN}求出样本联合分布p(X|θ)。
(3) 利用贝叶斯公式(38)，求θ的后验分布P(θ|X)。
(4) 结合损失函数计算贝叶斯估计量。
【例38】设X={x1,x2,…,xN}是来自总体分布为N(μ,σ2)的样本集，已知μ服从N(μ0,σ20)分布，采用平方误差损失函数，用贝叶斯估计的方法求μ的估计量μ^。
解： (1) 确定μ的先验分布P（μ）。

P（μ）=12πσ0exp－12μ－μ0σ02

(2) 求样本联合分布p(X|μ)。

pX|μ=∏Ni=1pxi|μ=∏Ni=112πσexp－12xi－μσ2

(3) 求μ的后验分布Pμ|X 。


P(μ|X)=pX|μ·P（μ）∫pX|μ·P（μ）dμ=α∏Ni=1pxi|μ·P（μ）

=α·12πσ0exp－12μ－μ0σ02∏Ni=112πσexp－12xi－μσ2

=α′·exp－12μ－μ0σ02∏Ni=1exp－12xi－μσ2

=α′·exp－12μ－μ0σ02+∑Ni=1xi－μσ2

=α″·exp－12Nσ2+1σ20μ2－21σ2·∑Ni=1xi+μ0σ20μ


P(μ|X)是μ的二次函数的指数函数，所以仍然是一个正态密度函数，把P(μ|X)写成N(μN,σ2N)的形式

P(μ|X)=12πσNexp－12μ－μNσN2

应用待定系数法，令两式对应的系数相等


1σ2N=Nσ2+1σ20

μNσ2N=Nσ2x-+μ0σ20
→
μN=Nσ20Nσ20+σ2x-+σ2Nσ20+σ2μ0

σ2N=σ20σ2Nσ20+σ2


其中，


x-=1N∑Nk=1xk



(4) 求贝叶斯估计量。由式(314)知μ^=∫μP(μ|X)dμ，即给定样本集下的μ的条件期望； 而P(μ|X)~N(μN,σ2N)，所以，μ的贝叶斯估计量为μ^=Nσ20Nσ20+σ2x-+σ2Nσ20+σ2μ0。
2. 推论分析
需要注意到，进行参数估计的根本目的是估计样本的概率密度函数p(x|X)，由于概率密度函数形式已知，才转换为估计概率密度函数参数的问题。因此，在贝叶斯估计中，在得到参数的后验概率后，可以不再去估计参数，直接根据P(θ|X)得到样本的概率密度函数。

p(x|X)=∫Θp(x|θ)P(θ|X)dθ(315)

其中，P(θ|X )如式(38)所示。
式(315)的含义是： 要估计的概率密度p(x|X)是所有可能的参数取值下样本概率密度p(x|θ)的加权平均，加权是在观测样本下估计出的随机变量θ的后验概率P(θ|X)。
假设P(θ|X)已知，由式(315)可知，p(x|X)是关于θ的p(x|θ)的期望，有

p(x|X)=EΘ［p(x|θ)］(316)

如果采样点θm足够多，m=1,2,…,M，可以用式(317)代替式(316)

p(x|X)≈1M∑Mm=1p(x|θm)(317)

如何生成一系列的采样点θm呢？如果Pθ|X已知，可以根据这个概率密度函数进行抽样生成。但Pθ|X的计算一般比较困难，可以使用马尔可夫链蒙特卡罗(Markov Chain Monte Carlo，MCMC)理论的方法，在不计算贝叶斯公式(38)的分母pX的情况下，对参数的后验概率Pθ|X进行随机抽样，生成采样点序列θm，进而实现概率密度函数的估计。Gibbs抽样和MetropolisHastings算法是两种最流行的方法，本书不详细介绍这些方法，有兴趣的读者可以查阅相关资料。
3. 贝叶斯学习
为了反映样本数目，重新标记样本集： XN={x1,x2,…,xN}，θ的贝叶斯估计量为

θ^=∫ΘθPθ|XNdθ(318)

其中，Pθ|XN为θ的后验分布

Pθ|XN=pXN|θ·P(θ)∫pXN|θ·P(θ)dθ(319)

样本独立抽取，当N＞1时，有

pXN|θ=pxN|θpXN-1|θ(320)

将式(320)代入式(319)得递推公式

Pθ|XN=pxN|θPθ|XN-1∫pxN|θPθ|XN-1dθ(321)

把先验概率记作Pθ|X0=P(θ)，表示在没有样本情况下的概率密度估计。随着样本数的增加，得到一系列对概率密度函数参数的估计

P(θ),P(θ|x1),Pθ|x1,x2,…,P(θ|x1,x2,…,xN),…(322)

称作递推的贝叶斯估计。随着样本数的增加，式(322)的后验概率序列逐渐尖锐，逐步趋向于以θ的真实值为中心的一个尖峰，当样本无穷多时收敛于在参数真实值上的脉冲函数，这一过程称作贝叶斯学习。
【例39】设X={x1,x2,…,xN}是来自总体分布为N(μ,σ2)的样本集，已知μ服从N(0,9)分布，σ=2，用贝叶斯学习的方法求μ的估计量μ^。
程序如下：

clc,clear,close all;

rng('default') 

N=300;                                             %样本数

mu0=0;   sigma0=3;                               %μ的先验分布参数

mu_value=-10:0.1:10;                            %μ的取值范围

Prior=pdf('Normal',mu_value,mu0,sigma0);    %μ的先验概率密度函数

plot(mu_value,Prior);

mu=2;   sd=2;                                     %样本服从的正态分布参数

training=normrnd(mu,sd,［N,1］);                %生成样本集

prev=Prior;

hold on

for i=1:N

tempdata=training(i);                       %当前样本xi

npdf=normalDist(tempdata,mu_value,sd);  %当前样本对应的p(xi |μ) 

numerator=npdf.*prev;                       %计算p(Xi|μ)=p(xi |μ)p(Xi-1│μ)P(μ)

prev=numerator;

TotalP=sum(numerator);                      %贝叶斯公式的分母

Poster=numerator/TotalP;                   %计算P(μ|Xi )

if i==10 || i==50 || i==100|| i==200 || i==300

plot(mu_value,Poster);                  %绘制后验概率曲线

end

end

［value,pos］=max(Poster);

plot(［mu_value(pos) mu_value(pos)］,［0 value］,'k-.');%绘制P(μ|XN)的尖峰

mu_value(pos)                                               %P(μ|XN)的尖峰所对应的μ

hold off 

function npdf=normalDist(X,Mu,Sigma)                   %计算p(xi |μ)的函数

Z=(X-Mu)./Sigma;

npdf=exp(-0.5*Z.^2)/(sqrt(2*pi)*Sigma);

end

运行程序，绘制的一系列概率密度函数如图33所示，可以看出，随着样本数N的增大，θ的后验概率密度函数逐渐尖锐。程序在命令窗口输出ans=2，表示后验概率序列最后一个的尖峰在μ=2处出现，即估计的μ为2，和设定一致。受计算机计算能力限制，程序中N取值有限。


图33一系列的后验概率估计图


3.3非参数估计
很多情况下对样本的分布并没有充分的了解，无法事先给出密度函数的形式，而且有些样本分布的情况也很难用简单的函数来描述。在这种情况下，就需要进行非参数估计，直接用样本估计出整个函数，即直接计算p(x)的值。
3.3.1直方图方法
1. 基本原理

以一维为例，解释直方图(Histogram)方法估计概率密度函数的含义。将数据x的取值范围分为若干个长度相等的区间，统计每个区间内的样本数，以样本数占总数的比例作为该区间内的概率密度值，即直方图方法。例如，图像的灰度直方图，就是统计图像中一维灰度数据的概率密度。
如果推广到n维，采用直方图方法估计概率密度函数的做法是： 
(1) 把n维样本的每个分量在其取值范围内分成m个等间隔的小窗，则这种分割会得到mn个小舱，每个小舱的体积记作V。
(2) 统计落入每个小舱内的样本数目k。
(3) 把每个小舱内的概率密度看作常数，并用k/NV作为其估计值(N为样本数)，即小区域范围内概率密度函数为

p^(x)=kNV(323)

2. 仿真实现
MATLAB提供了统计、绘制数据直方图的函数histcounts、histogram、histcounts2和histogram2，主要调用格式如下： 
(1) histogram(X,M)： 创建并以柱状形式显示X的直方图。M指定直方图中柱的数目，每个柱宽度相等，省略M时根据X中元素的范围和分布形状计算； 柱的高度表明该柱内元素的数目。X可以是向量、矩阵或多维矩阵，但把矩阵作为列向量X(:)，绘制一个直方图。
(2) histogram(X,EDGES)： 利用向量EDGES指定柱状图各个柱的边界。EDGES中指定每个柱的左边界，但不一定有右边界，最后一个柱拥有左右边界。
(3) histogram('BinEdges', EDGES, 'BinCounts', COUNTS)： 指定每个柱的高度COUNTS和柱的边界EDGES实现直方图的绘制。EDGES和COUNTS可以是histcounts函数的返回值。
(4) histogram(…,NAME,VALUE)： 设定参数统计并显示直方图。主要参数如表31所示。
(5) ［N,EDGES］=histcounts(X,M)： 将X中的数据均分成M柱，并且返回每个柱中的数目N以及边界EDGES。省略M时由算法自动确定。
(6) ［N,EDGES］=histcounts(X,EDGES)： 通过指定柱状图各个柱的边界进行统计。
(7) ［N,EDGES］=histcounts(…, NAME,VALUE)： 设定参数统计直方图,主要参数如表31所示。
(8) histogram2(X,Y)： 绘制二维直方图，X和Y大小一致，其余参数同histogram函数。
(9) ［N,XEDGES,YEDGES］=histcounts2(X,Y)： 统计二维直方图，返回各柱内数目N以及X和Y方向的边界，其余参数同histcounts函数。


表31histcounts和histogram函数主要参数



参数名称取值及含义
BinWidth指定每个柱的宽度。最多65536个柱，若BinWidth太小，65536个柱不够，则调整柱的宽度
BinLimits二维向量［BMIN,BMAX］，指定参与统计的元素范围，即统计X(X＞=BMIN & X＜=BMAX)的数目
Normalization指定每个柱的高度表示的含义： 取'count'，默认值，表示柱内元素的数目； 取'probability'，表示概率，柱内元素数/总元素数； 取'countdensity'，表示柱内元素数/柱的宽度； 取'pdf'，表示概率密度估计，柱内元素数/(总元素数×柱的宽度)； 取'cumcount'，表示累积数，当前柱及之前柱内元素数； 取'cdf'，表示累积分布函数，当前柱及之前柱内元素数/总元素数
DisplayStylehistogram函数参数，指定绘制直方图风格： 取'bar'，默认，以柱状图显示直方图； 取'stairs'，以阶梯状折线绘制直方图外轮廓线，但不填充内部


【例310】设定参数生成正态分布数据集，利用直方图方法估计概率密度函数并绘制函数曲线。
程序如下：

clc,clear,close all;

rng('default')

N=200;

mu=0;sd=0.8;                          %指定正态分布参数

x=normrnd(mu,sd,［N,1］);                   %生成数据集

x_values=-3:0.01:3;

px=pdf('Normal',x_values,mu,sd);        %计算概率密度函数

plot(x_values,px,'Color','k');          %绘制概率密度函数曲线

［hist1,edge1］=histcounts(x,10,'Normalization','pdf');%设定10个区间进行估计，大间隔

［hist2,edge2］=histcounts(x,30,'Normalization','pdf');%设定30个区间进行估计，小间隔

hold on

histogram('BinEdges',edge1,'BinCounts',hist1,'FaceColor','w');%绘制柱状图

hold off

figure,plot(x_values,px,'Color','k','LineWidth',2);

hold on

histogram(x,30,'Normalization','pdf');                             %统计并绘制柱状图

hold off

运行程序，绘制的概率密度函数曲线如图34(a)和图34(b)所示； 修改样本数为1000，修改倒数第二条命令为histogram(x,30,'Normalization','pdf','DisplayStyle','stairs')，绘制阶梯状折线，如图34(c)所示； 样本数为2000时，绘制的阶梯状折线如图34(d)所示。同等间隔，随着样本数的增多，估计的效果逐渐提高。所以，估计概率密度函数本身需要充足的样本。


图34利用直方图方法估计概率密度函数


3. 方法分析
直方图方法估计效果与小舱的选择密切相关。当小舱较大时，每个小舱内设概率密度为常数，估计出的密度函数粗糙，如图34(a)所示； 当小舱较小时，有些小舱内可能会没有样本或者样本很少，导致估计出的概率密度函数不连续，如图34(b)所示。所以，小舱的大小要合理选择。
小舱的选择应与样本总数相适应。直方图方法在每个小舱内用常量表示概率密度函数值，当样本数N→∞时，近似值收敛于真实值的条件是： 随样本数的增加，小舱体积应该尽可能小，同时必须保证小舱内有充分多的样本，但每个小舱内的样本数又必须是总样本数中很小的一部分，用公式表示为

limN→∞VN=0limN→∞kN=∞limN→∞kNN=0(324)

VN、kN表示小舱体积V和小舱内样本数k与样本数N有关。
小舱内的样本数与样本分布有关。在有限样本数目下，如果小舱体积相同，样本密度大的地方小舱内有较多的样本，密度小的地方小舱内样本很少甚至没有，将导致密度的估计在样本密度不同的地方表现不一致。所以，最好能够根据样本分布情况调整小舱体积。
小舱数目随着样本维数的增加急剧增多，实际中很难有足够的样本数据用于估计。例如，如果样本为100维，每维分为4个间隔，将有4100≈1.6×1060个小舱，即使要保证一个小舱内有一个样本，也需要1.6×1060个样本，这是很难达到的，所以，估计时很多小舱内将会没有样本，方法失效。但在低维问题中，直方图方法还是一个用来进行建模和可视化的有效工具。
3.3.2Parzen窗法
Parzen窗法又称核密度估计(Kernel Density Estimation，KDE)，采用滑动小舱估计每一点的概率密度。
1. 基本原理
样本x∈Rn，假设每个小舱是一个超立方体，在每一维的棱长都为h，则小舱的体积是

V=hn(325)

定义n维单位方窗函数为

φ(［u1u2…un］T)=1，|uj|≤12,j=1,2,…,n

0，其他(326)

该式的含义是： 在以原点为中心的n维单位超立方体内的所有点的函数值为1，其他点的值为0，如图35(a)所示。


图35二维的方窗函数示意图


如果一个样本xi在以x为中心的棱长为h的超立方体内，xi到x的每一维距离都小于或等于h/2，即向量|x－xi|/h的每一维都小于或等于1/2，则φ［(x－xi)/h］=1，否则φ［(x－xi)/h］=0，如图35(b)所示。将统计落入以x为中心的超立方体内的样本数，转化为统计φ［(x－xi)/h］=1的次数，即


k=∑Ni=1φx－xih(327)

将式(327)代入式(323)，得任意一点x的密度估计表达式

p^(x)=1NV∑Ni=1φx－xih(328)

或

p^(x)=1N∑Ni=11Vφx－xih(329)

定义核函数(窗函数)为

Kx,xi=1Vφx－xih(330)

反映了一个观察样本xi对在x处的概率密度估计的贡献与样本xi和x的距离有关。概率密度估计即在每一点上把所有观测样本的贡献进行平均

p^(x)=1N∑Ni=1Kx,xi(331)

所以，这种用窗函数(核函数)估计概率密度的方法被称作Parzen窗方法或核密度估计，参数h也常被称为带宽(Bandwidth)、窗宽或平滑参数。
2. 核函数选择
由式(331)可知，如果通过定义核函数直接估计概率密度，要求估计出的概率密度函数满足： 非负且积分为1，也就是要求核函数满足

K(x,xi)≥0且∫K(x,xi)dx=1(332)

常用的核函数如下。
1) 方窗函数
综合式(326)和式(330)，得方窗核函数

K(x,xi)=1hn，|xj－xji|≤h2,j=1,2,…,n

0，其他(333)

2) 高斯窗函数
一维的单位高斯窗函数为

φ(u)=12πexp－12u2(334)

得一维情况下的高斯核函数为

Kx,xi=12πhexp－x－xi22h2(335)

假设样本x属性条件独立，n维情况下的高斯核函数表示为

K(x,xi)=12πn/2hnexp－x－xiTx－xi2h2(336)

3) 三角窗函数
一维的三角窗函数为

φ(u)=1－|u|，|u|≤1

0，其他(337)

得一维情况下的三角核函数

Kx,xi=1h1-x－xih，x－xi≤h

0，其他(338)

n维情况下的三角核函数为

Kx,xi=1hn1-x－xih，|xj－xji|≤h,j=1,2,…,n

0，其他(339)

4) Epanechnikov窗函数
一维的Epanechnikov窗函数为

φ(u)=3451－u25，|u|≤5

0，其他(340)

得一维情况下的Epanechnikov核函数为

Kx,xi=345h1－x－xi25h2，x－xi≤5h

0，其他(341)

设样本x属性条件独立，n维情况下的Epanechnikov核函数为


Kx,xi=345hn1-x－xiTx－xi5h2，|xj－xji|≤5h,j=1,2,…,n

0，其他(342)

Epanechnikov核函数被证明是最小平方误差意义上的最优核函数，证明可参看非参数统计类书籍； 实践中高斯核应用更广泛。
3. 带宽的选择
带宽h对估计有很大的影响，带宽的选择是核密度估计中一个很重要的问题。
估计的密度函数p^(x)=1VN∑Ni=1φx-xih，由V=hn可知： h较大时，1Vφx-xih的幅度较小，宽度较大，p^(x)是N个宽度较大的缓慢变化函数的迭加，估计值p^(x)受到过度的平均作用，p(x)比较细致的性质不能显露，估计分辨率较低。
如果h较小，1Vφx-xih的峰值较大，宽度较小，p^(x)就变成N个以xi为中心的尖脉冲的迭加，呈现不规则的形状，从而使估计不稳定。
h应随着样本数N增大而缓慢下降，从理论上来说，随着N→∞而趋于零。实际问题中样本数有限，只能进行折中选择。
【例311】设待估计的p(x)是均值为零、方差为1的正态密度函数，用Parzen窗法估计p(x)。
设计思路： 
(1) 确定窗函数： 选择正态窗函数并确定窗宽，窗函数为

Kx,xi=12πhexp－x－xi22h2

(2) 计算估计值

p^(x)=1N∑Ni=1Kx,xi=1hN2π∑Ni=1exp－x－xi22h2

当采集到样本后，可以根据上式计算估计式p^(x)。
程序如下：

clc,clear,close all;

rng('default')

N=1000; mu=0; sigma=1;

x=normrnd(mu,sigma,［N,1］);                               %生成1000个样本

minx = min(x);    maxx = max(x);    dx=(maxx-minx)/N;

x_values=minx:dx:maxx-dx;                                %确定采样点x

px=pdf('Normal',x_values,mu,sigma);                    %确定概率密度函数

plot(x_values,px,'Color','k','LineWidth',2);         %绘制原始的概率密度函数曲线

hold on

h=0.01;                                                      %设置较小的带宽

pxe1=kde(x,x_values,h,N);                                %估计概率密度函数

plot(x_values,pxe1,'r:','LineWidth',2);

h=2;                                                          %设置较大的带宽

pxe3=kde(x,x_values,h,N);

plot(x_values,pxe3,'b--','LineWidth',2);

xlabel('x'),ylabel('p(x)');

legend('p(x)','h=0.01','h=2');

hold off

figure,plot(x_values,px,'Color','k','LineWidth',2);

h=0.3;                                                       %设置较适中的带宽

pxe2=kde(x,x_values,h,N);

hold on

plot(x_values,pxe2,'g--','LineWidth',2);

xlabel('x'),ylabel('p(x)');

legend('p(x)','h=0.3');

hold off

function pxe=kde(x,x_values,h,N)                        %利用高斯窗进行估计的函数

pxe=zeros(1,N);

for j=1:N

for i=1:N

pxe(j)=pxe(j)+exp(-0.5*(x_values(j)-x(i))^2/h^2)/sqrt(2*pi);

end

pxe(j)=pxe(j)/N/h;

end

end

程序运行结果如图36所示。图36(a)中实线为原始的概率密度函数曲线； 点线为带宽h=0.01时估计的概率密度函数曲线，h较小，估计的曲线相对于原始曲线起伏较大，估计不稳定； 虚线是h=2时估计的概率密度函数曲线，h较大，估计的曲线相对于原始曲线较平滑，但跟不上函数p(x)的变化。图36(b)中虚线为带宽h=0.3时估计的概率密度函数曲线，估计的较准确。正如前面分析时所述，带宽h对估计结果影响较大。


图36利用高斯核函数估计概率密度函数


【例312】产生1/64/256/10000个服从一维标准正态分布的样本，用带宽为h=1/N的Epanechnikov核函数估计样本的概率密度函数。
程序如下：

clc,clear,close all;

rng('default')

N=［1 64 256 10000］;                           %设置不同的样本数

mu=0;sigma=1;

x_values=-4:0.01:4;                           %确定采样点

h=1./(N.^0.5);                                 %设置带宽随样本数增多逐渐下降

R=zeros(length(N),N(4));

for m=1:length(N)

R(m,1:N(m))=normrnd(mu,sigma,1,N(m)); %生成样本矩阵，每一行对应一种样本数

end

px=pdf('Normal',x_values,mu,sigma);        %真实的概率密度函数

len=length(x_values);

for m=1:length(N)                              %针对不同样本数分别估计概率密度函数

pxe=zeros(1,len);

for i=1:len

 for j=1:N(m)

if abs(x_values(i)-R(m,j))<=sqrt(5)*h(m) 

pxe(i)=pxe(i)+(1-((x_values(i)-R(m,j))/h(m))^2/5)*3/4/sqrt(5)/h(m);

end                                   %采用Epanechnikov核函数进行估计

end

pxe(i)=pxe(i)/N(m);

end

subplot(1,4,m),plot(x_values,px,'k');  %绘制真实的概率密度函数

hold on    

plot(x_values,pxe,'r--','LineWidth',2),axis(［-3,3,0.001,1.0］);%绘制估计的密度函数

str=strcat('N=',num2str(N(m)));

legend('p(x)',str);

hold off

end

程序运行结果如图37所示，图中实线为原始的概率密度函数曲线； 虚线为估计出的概率密度函数曲线，可以看出，样本数越多估计效果越好。


图37利用Epanechnikov核函数估计概率密度函数


4. 仿真实现
MATLAB提供了进行核密度估计的函数，列举如下。
1) fitdist函数
在最大似然估计中介绍的fitdist函数，其主要调用格式为： PD=fitdist(X,DISTNAME)。参数DISTNAME指定拟合数据的概率分布类型，取'kernel'时拟合非参数核平滑分布。有3个专有参数： 'kernel'指定核函数，值为'normal'(默认)、'box'、'triangle'、'epanechnikov'； 'support'指定是否限制密度值范围，值为'unbounded'(默认)和'positive'，或者二维向量； 'width'指定带宽。
2) ksdensity函数
用于估计一维或二维密度函数，调用格式如下。
(1) ［F,XI,U］=ksdensity(X)： 使用高斯核函数对X中的数据进行概率密度估计，返回采样点XI和对应的密度值F。X为向量时，计算100个覆盖数据范围的点(XI)的密度估计值(F)； X为两列矩阵时，使用meshgrid从每维等间隔采样30个点，即XI含900个点，F为向量。带宽U由X的点数和维数决定。
(2) F=ksdensity(X,XI)： 指定要估计密度值的点向量或两列矩阵XI进行估计。
(3) ksdensity(...)： 无返回值，绘制估计的概率密度函数曲线。
(4) ［...］=ksdensity(...,'PARAM1',val1,'PARAM2',val2,...)： 指定参数实现概率密度函数估计。主要参数如表32所示。


表32ksdensity函数主要参数



参数名称取值及含义
kernel指定核函数，可取'normal'(默认)、'box'、'triangle'和'epanechnikov'
support指定是否限制密度值范围，值为'unbounded'(默认)和'positive'，或者以二维向量为单变量数据、2×2的矩阵为二维变量指定密度上下限

weights和X等长的向量，指定每个元素的权值，默认权值相等
bandwidth指定核函数带宽，对于二维情况，可以是二维向量； 默认为估计正态密度的最优值
function要估计的函数类型，可选'pdf'(默认)、'cdf'和'survivor'； 对于单变量，也可选'icdf'或'cumhazard'

3) mvksdensity函数
用于多元核密度估计，调用格式如下。
(1) F=mvksdensity(X,XI,'Bandwidth',BW)： 基于高斯核函数，根据N×n矩阵X中的采样数据，计算XI指定点的概率密度估计值F，BW为n维向量，指定带宽。
(2) F=mvksdensity(...,'PARAM1',val1,'PARAM2',val2,...)： 指定参数实现估计，参数同ksdensity函数，见表32。
【例313】读取hospital数据集，按照性别分为两类，采用ksdensity函数分别对两类体重数据进行概率密度函数估计，并绘制估计的概率密度函数曲线。
程序如下：

clc,clear,close all;

load hospital

x=hospital.Weight;

gender=hospital.Sex;

x_female=x(gender=='Female');                  %获取两类体重数据

x_male=x(gender=='Male');

［pxe_female,xv_female］=ksdensity(x_female,'kernel','epanechnikov');

［pxe_male,xv_male］=ksdensity(x_male,'kernel','epanechnikov');  %指定核函数进行估计

plot(xv_female,pxe_female,'Color','r','LineWidth',2)

hold on

plot(xv_male,pxe_male,'Color','b','LineStyle','--','LineWidth',2);

xlabel('x'),ylabel('p(x)');

legend('Female','Male')

hold off

程序运行结果如图38所示。


图38采用ksdensity函数对hospital中的体重数据进行概率密度函数估计


3.3.3kN近邻密度估计法
用直方图方法和Parzen窗法估计概率密度函数时，点x周围的体积固定为V=hn，其中的样本点数kN随机变化。kN近邻密度估计法采用可变大小的小舱，固定每个小舱内拥有的样本个数kN，调整包含x的小舱的体积，直到小舱内恰好落入kN个样本，再估算概率密度函数，表示为

p^(x)=kNNV(x)(343)

V(x)依赖于x，在低密度区，体积大，在高密度区，体积小。
实际中，计算所有样本点之间的距离，对于每一个样本点，确定其周围kN个近邻，以到最远近邻点的距离为半径r，得到超球体并计算体积。如果采用Mahalanobis距离，将会得到超椭球体，体积定义为

V=V0|Σ|12rn(344)

其中，V0是半径为1的超球体的体积，定义为

V0=πn2/n2!，n为偶数

2nπn－12n－12!/n!，n为奇数(345)

可以证明，对于二维空间的圆，面积为πr2； 三维空间的球，体积为4πr3/3。
【例314】设定参数生成正态分布数据集，利用kN近邻法估计概率密度函数并绘制函数曲线。
程序如下：

clc,clear,close all;

rng('default')

N=1000;    mu=0;    sigma=0.8;

x=normrnd(mu,sigma,［N,1］);

kn1=10;    kn2=50;    kn3=100;              %设定三个不同的kN

x_values=-3:0.01:3;

px=pdf('Normal',x_values,mu,sigma);

len=length(x_values);

pxe1=zeros(1,len);   pxe2=zeros(1,len);   pxe3=zeros(1,len);

index=1;

for j=-3:0.01:3

distance=pdist2(j,x);                    %计算当前点和样本点之间的距离

D=sort(distance);                         %距离升序排列

V1=2*D(kn1);V2=2*D(kn2);V3=2*D(kn3);  %第kN个距离作为半径，一维空间，确定长度2r

pxe1(index)=kn1/N/V1;

pxe2(index)=kn2/N/V2;

pxe3(index)=kn3/N/V3;                    %按式(3-43)估计密度值

index=index+1;

end

figure,plot(x_values,px,'Color','k'); 

hold on

plot(x_values,pxe1,'r:');                 %绘制第一种kN情况下的对比图

hold off

figure,plot(x_values,px,'Color','k');

hold on

plot(x_values,pxe2,'g--');                %绘制第二种kN情况下的对比图

hold off

figure,plot(x_values,px,'Color','k');

hold on

plot(x_values,pxe3,'b-.');                %绘制第三种kN情况下的对比图

hold off

程序运行结果如图39所示。可以看出，kN的取值影响到估计的结果，取值大估计相对准确。可以自行更改样本数目或kN值，观察估计结果。


图39采用kN近邻法进行概率密度函数估计


非参数估计的共同问题是对样本数目需求较大，只要样本数目足够大，总可以保证收敛于任何复杂的位置密度，但计算量和存储量都比较大。当样本数很少时，若对密度函数有先验认识，参数估计方法能取得更好的估计效果。
3.4最小错误率贝叶斯决策的实例
【例315】MATLAB中的fisheriris数据集，对3种鸢尾花(setosa，versicolor和virginica)各抽取了50个样本，每个样本含四个特征，分别为花萼长、宽，花瓣长、宽，单位为厘米。对样本x=4.83.51.50.2T进行最小错误率贝叶斯决策。
1. 设计思路
贝叶斯决策需要先了解先验概率和类条件概率密度。题目中仅给出了样本集，数据集中有3类，样本数相同，可以假设各类先验概率相等，但需要估计类条件概率密度函数。对样本进行最小错误率贝叶斯决策，只需要比较样本对应的类条件概率密度的大小即可归类。

由于并不清楚数据服从什么分布，需要对类条件概率密度进行非参数估计。样本为四维，若每维取10个等距点，则四维空间有104个等距点，但每类样本仅50个，样本数量太少，难以支撑四维空间的概率密度函数估计，因此，考虑将样本降维以保证估计的效果。本题中简单取原始样本的第三维和第四维构成新的样本。
综上所述，得设计流程： 首先对于输入的数据降维，再进行非参数概率密度函数估计，最后比较样本对应的各类的概率密度函数值的大小，将样本归入取值最大的类。
2. 程序设计

clc,clear,close all;

load fisheriris                                              %导入数据集

training=meas(:,3:4);                                       %数据降维

X_se=training(ismember(species,'setosa'),:);

X_ve=training(ismember(species,'versicolor'),:);

X_vi=training(ismember(species,'virginica'),:);        %分别表示3类数据

xmin=min(training,［］,1);

xmax=max(training,［］,1);                                    %确定估计点取值范围

dx=0.1;

［x1,x2］=ndgrid(xmin(:,1):dx:xmax(:,1),xmin(:,2):dx:xmax(:,2)); %构建二维网格

dim=size(x1');

x1=x1(:,:)';    

x2=x2(:,:)';    

XI=［x1(:) x2(:)］;                                  %确定二维空间的采样点

Pxe_se=mvksdensity(X_se,XI,'Bandwidth',1);

Pxe_ve=mvksdensity(X_ve,XI,'Bandwidth',1);

Pxe_vi=mvksdensity(X_vi,XI,'Bandwidth',1);    %用核密度估计法对各类概率密度进行估计

Pxe_se=reshape(Pxe_se,dim);

Pxe_ve=reshape(Pxe_ve,dim);

Pxe_vi=reshape(Pxe_vi,dim);                      %将一维概率密度向量表达为二维形式

subplot(131),mesh(xmin(:,1):dx:xmax(:,1),xmin(:,2):dx:xmax(:,2),Pxe_se);

subplot(132),mesh(xmin(:,1):dx:xmax(:,1),xmin(:,2):dx:xmax(:,2),Pxe_ve);

subplot(133),mesh(xmin(:,1):dx:xmax(:,1),xmin(:,2):dx:xmax(:,2),Pxe_vi);

X=［4.8 3.5 1.5 0.2］;

X=X(3:4);

pos=round((X-xmin)/dx+1);

P_se=Pxe_se(pos(2),pos(1));

P_ve=Pxe_ve(pos(2),pos(1));

P_vi=Pxe_vi(pos(2),pos(1));                %获取待决策样本在3类概率密度函数中的取值

maxP=max(［P_se P_ve P_vi］);

if P_se==maxP                                 %比较大小并归类

disp('样本为setosa类');    

elseif P_ve==maxP

disp('样本为versicolor类'); 

else

disp('样本为virginica类');

end

count_se=0;    

count_ve=0;   

count_vi=0;  

for i=1:150                                   %对所有样本进行测试，计算总体识别率

X=training(i,:);

pos=round((X-xmin)/dx+1);

P_se=Pxe_se(pos(2),pos(1));

P_ve=Pxe_ve(pos(2),pos(1));

P_vi=Pxe_vi(pos(2),pos(1));

maxP=max(［P_se P_ve P_vi］);

if P_se==maxP && i<=50

count_se=count_se+1;

elseif P_ve==maxP && i>50 && i<=100

count_ve=count_ve+1;

elseif P_vi==maxP && i>100 

count_vi=count_vi+1;

end

end

ratio=(count_se+count_ve+count_vi)/150

3. 实验结果
运行程序，绘制的概率密度函数图形如图310所示，且将在命令窗口输出：

样本为setosa类

ratio = 0.9600



图310估计的概率密度函数


习题
1. 简述参数估计和非参数估计的含义。
2. 列举参数估计的方法并简述各方法的含义。
3. 列举非参数估计的方法并简述各方法的含义。
4. 设样本集呈现对数正态分布p(x)=1σx2πexp－lnx－θ22σ2,(x＞0)，求θ的最大似然估计。
5. 设正态分布的均值分别为μ1=11T和μ2=1.51.5T，协方差矩阵均为0.2I，先验概率相等，决策表为λ=01

0.50。编写程序，由正态分布生成各1000个二维向量的数据集，利用其中的800个样本，采用最大似然估计方法估计样本分布的参数，利用最小风险贝叶斯决策方法对其余200个样本进行决策，并计算识别率。
6. 打开hospital数据集，分别采用直方图方法、Parzen窗法和kN近邻法对其中的血压数据(二维)进行概率密度函数估计。