MATLAB核密度估计「终于解决」

核密度估计

在很多统计问题中，需要由样本去估计总体的概率分布密度，常用的估计方法由参数法和非常数法。参数法是假定总体服从某种已知的分布，即密度函数的形式是已知的，需要由样本估计其中的参数，这种方法依赖于实现对总体分布的假设，而做出这种假设往往是非常困难的。非参数法则不存在这样的“假设”困难，这里介绍的就是一种非参数密度估计法–核密度估计。

核密度估计需要指定核函数和窗宽，但是取不同的核函数对核密度估计影响不大。

（1）常用的核函数

Uniform（或Box）、Triangle、Epanechnikov、Quaritic、Triweight、Gaussian、Cosinus

（2）窗宽对核密度估计的影响

窗宽会影响光滑程度，如果窗宽h去较大的值，图形较为光滑，但同时也丢失了数据所包含的一些信息；如果窗宽取值较小，则图像是不光滑的曲线，但它能反映出每个数据所包含的信息。

（4）核密度估计的MATLAB实现

MATLAB统计工具箱中提供了ksdensity函数，用来求核密度估计，其调用格式如下：

<1>[f,xi]=ksdensity(x)

求样本观测向量x的核密度估计，xi是在x取值范围内等间隔选取的100个点构成的向量，f是与xi相对应的核密度估计值向量。在在所用的核函数是Gaussian核函数，窗宽也是默认值。

<2>f=ksdensity(x,xi)

根据样本观测向量x计算xi出的核密度估计值f，xi和f是等长的向量。

<3>ksdensity(…….)

不返回任何输出，此时在当前坐标系中绘制出核密度函数图。

<4>ksdensity(ax,…….)

不返回任何输出，此时在句柄在ax对应的坐标系中绘制出核函密度函数图。

<5>[f,xi,u]=ksdensity(…….)

返回窗宽u

<6>[……]=ksdensity（…..，param1，val1，param2，val2，….）

通过可选的成对出现的参数名与参数值来控制核密度估计，可用的参数名与参数值如下表

参数名 参数值 说明

‘censoring’ 与x等长的逻辑向量 指定哪些项是截尾观测，默认是没有截尾

‘normal’ 指定用Gaussian（高斯或正态）核函数，默认情况

‘box’ 指定用Uniform核函数

‘bernel’ ‘triangle’ 指定用Triangle核函数

‘epanechnikov’ 指定用Epanechnikov核函数

函数句柄或函数名 自定义核函数

‘npoints’ 正整数 指定xi中包含的等间隔点的个数，默认100

‘unbounded’ 指定密度函数的支撑集合为全体实数集，默认情况

‘support’ ‘positive’ 指定密度函数的支撑集为正实数集

包含两个元素的向量 指定密度函数的支撑集的上下限

‘weights’ 与x等长的向量 指定x中元素的权重

‘pdf’ 指定对密度函数进行估计

‘cdf’ 指定对累积分布函数估计

’function’ ‘icdf’ 指定对逆概率分布函数估计

‘survior’ 指定对生存函数进行估计

‘cumhazard’ 指定对累积危险函数进行估计

（5）核密度估计的例子

调用ksdensity函数对总成绩数据进行密度函数估计，并通过改变窗宽和核函数类型，来观察窗宽参数对函数平滑程度的影响和观察核函数类型对估计结果的影响。

<1>总成绩数据的核密度估计

采用默认的最佳窗宽和默认的Gaussian核函数，调用ksdensity函数进行核密度估计，并将核密度估计图、频率直方图和前面章节中求出的总成绩的正态分布的密度函数图放在一起对比。

%读取文件中的第一个工作表中的总成绩数据，即G2：G52

score=xlsread(‘成绩.xls’,’G2:G52′);

score=score(score>0);

%调用ecdf函数计算xc处的经验分布函数值f_ecdf

[f_ecdf,xc]=ecdf(score);

%新建图形窗口，然后绘制评论直方图，直方图对应7个小区间

figure;

ecdfhist(f_ecdf,xc,7);

hold on;

xlabel(‘考试成绩’);

ylabel(‘f(x)’);

%调用ksdensity函数进行核密度估计

[f_ks1,xi1,u1]=ksdensity(score);

%绘制核密度估计图，并设置线条为黑色实线，宽度为3

plot(xi1,f_ks1,’k’,’linewidth’,3);

%计算正态分布的密度函数图

ms=mean(score); %均值

ss=std(score); %方差

%计算xi1处的正态分布密度函数值，正态分布的均值是ms，方差是ss

f_norm=normpdf(xi1,ms,ss);

%绘制正态分布密度函数图，并设置线条颜色为红色点画线，宽3

plot(xi1,f_norm,’r-.’,’linewidth’,3);

%给图形加入标注框

legend(‘频率直方图’,’核密度估计图’,’正态分布密度图’);

%查看默认的窗宽u1

u1

u1 =

5.0474

可以看到，在默认窗宽下，利用Gaussian核函数求出的密度曲线与N（ms，ss）分布的密度曲线非常接近，与总成绩的频率直方图附和的也很好。

<2>固定核函数为Gaussian和函数，让窗宽进行变动，观察不同的窗宽对核密度估计的影响。

%读取文件中的第一个工作表中的总成绩数据，即G2：G52

score=xlsread(‘成绩.xls’,’G2:G52′);

score=score(score>0);

%设置窗宽分别为0.1,1,5，9，调用ksdensity函数进行核密度估计

[f_ks1,xi1]=ksdensity(score,’width’,0.1);

[f_ks2,xi2]=ksdensity(score,’width’,1);

[f_ks3,xi3]=ksdensity(score,’width’,5);

[f_ks4,xi4]=ksdensity(score,’width’,9);

figure;

%分别绘制不同窗对应的核密度估计图，他们对应不同的线型和颜色

plot(xi1,f_ks1,’c-.’,’linewidth’,2);

hold on;

xlabel(‘考试成绩’);

ylabel(‘核密度估计’);

plot(xi2,f_ks2,’r:’,’linewidth’,2);

plot(xi3,f_ks3,’k’,’linewidth’,2);

plot(xi4,f_ks4,’b–‘,’linewidth’,2);

%加标注

legend(‘窗宽为0.1′,’窗宽为1′,’窗宽为5′,’窗宽为9’);

由图可以发现，不同的窗宽下，核密度估计曲线形状差距比较大，对于比较小的窗宽值，核密度估计曲线比较曲折，光滑性很差，但是反映了较多的细节；对于比较大的窗宽值，核密度估计曲线比较光滑，但是掩盖了许多细节。

<3>固定窗宽为默认的最佳窗宽，让核函数变动，观察不同核函数对核密度曲线估计的影响。

%读取文件中的第一个工作表中的总成绩数据，即G2：G52

score=xlsread(‘成绩.xls’,’G2:G52′);

score=score(score>0);

%设置核函数分别为Gaussian，Uniform，Triangle和Epanechnikov

%调用ksdensity函数进行核密度检验

[f_ks1,xi1]=ksdensity(score,’kernel’,’normal’);

[f_ks2,xi2]=ksdensity(score,’kernel’,’box’);

[f_ks3,xi3]=ksdensity(score,’kernel’,’triangle’);

[f_ks4,xi4]=ksdensity(score,’kernel’,’epanechnikov’);

figure;

%分布绘制不同核函数所对应的核密度估计图

plot(xi1,f_ks1,’k’,’linewidth’,2);

hold on;

plot(xi2,f_ks2,’r:’,’linewidth’,2);

plot(xi3,f_ks3,’b-.’,’linewidth’,2);

plot(xi4,f_ks4,’c–‘,’linewidth’,2);

xlabel(‘考试成绩’);

ylabel(‘核密度估计’);

legend(‘Gaussian’,’Uniform’,’Triangle’,’Epanechnikov’);

通过上图可以看出，不同的核函数对核密度估计的影响不大，就光滑性而言，Gaussian和Epanechnikov核函数对应的光滑性较好，Triangle次之，Uniform最差，在应用中，一般用Gaussian核函数。