《模式识别原理与应用》非参数估计 第四周

1、模式识别原理与应用非参数估计 第四周一、 密度估计上图为一个N(-5,4)正态分布和N(5,2)正态分布相叠加的分布，后面将围绕这个分布分别使用直方图密度估计和核密度估计。1、基于直方图密度估计基于直方图密度估计的缺点：样本量要大，间隔控量，不光滑，下面为直方图密度估计的代码：clear;% 绘制原始图像X=(-25:0.002:15);Y=(normpdf(X,-5,4)+normpdf(X,5,2)/2;figure(1);plot(X,Y);% 随机生成1000个N(-5,4)正态分布和1000个N(5,2)正态分布rndDistribute=normrnd(-5,4,1000,1);n

2、ormrnd(5,2,1000,1);n=2000; % 直方图生成dif=1; %直方图间隔n=(max(rndDistribute)-min(rndDistribute)/dif; %直方图条数H X=hist(rndDistribute,n);H=H/sum(H)/dif;figure(2);bar(X,H);结果：2、核密度估计概率密度函数为pnx=knnVn=1n i=1n1VnK(x-xihn) ，其中K(u)为核函数，这里给出两种常用的核函数。利用核函数Ku=1u120otherwise得到的结果：其中，hn分别为0.5，0.8，1，2。利用核函数Kz=12exp-z22得到的结

3、果：其中，hn分别为0.5，0.8，1，2。这里程序定义了两个函数：Kernel.m，Density.m。Kernel.m主要是通过传入的method参数来分别调用这两种核函数，Density.m是计算核密度估计公式pnx=1n i=1n1VnK(x-xihn)主程序就是根据不同的参数调用Density并绘制图像。源代码如下：function K = Kernel( ui,method )%KERNEL 01核 高斯核 if strcmp(method,'01') t=zeros(size(ui); t(abs(ui)<=0.5)=1; K=sum(t); elseif

4、strcmp(method,'GAS') t=1/sqrt(2*pi).*exp(-ui.2/2); K=sum(t); endendfunction kn = Density(rndDistribute, x,hn,method )%公式计算 kn=zeros(size(x); for i=1:1:size(x,2); kn(i)=Kernel(rndDistribute-x(i)/hn,method)/hn; end kn=kn/size(rndDistribute,1);endclear;% 绘制原始图像X=(-25:0.002:15);Y=(normpdf(X,-5,4

5、)+normpdf(X,5,2)/2;figure(1);plot(X,Y);% 随机生成1000个N(-5,4)正态分布和1000个N(5,2)正态分布rndDistribute=normrnd(-5,4,1000,1);normrnd(5,2,1000,1);n=2000;% 核密度估计% 01核 和 高斯核% hn 分别为 0.5 0.8 1 2% 01核figure(2);x=-25:0.01:15;hn=0.5;kn=Density(rndDistribute,x,hn,'01');subplot(2,2,1);plot(x,kn);hn=0.8;kn=Density

6、(rndDistribute,x,hn,'01');subplot(2,2,2);plot(x,kn);hn=1;kn=Density(rndDistribute,x,hn,'01');subplot(2,2,3);plot(x,kn);hn=2;kn=Density(rndDistribute,x,hn,'01');subplot(2,2,4);plot(x,kn);% 高斯核figure(3);x=-25:0.01:15;hn=0.5;kn=Density(rndDistribute,x,hn,'GAS');subplot(2

7、,2,1);plot(x,kn);hn=0.8;kn=Density(rndDistribute,x,hn,'GAS');subplot(2,2,2);plot(x,kn);hn=1;kn=Density(rndDistribute,x,hn,'GAS');subplot(2,2,3);plot(x,kn);hn=2;kn=Density(rndDistribute,x,hn,'GAS');subplot(2,2,4);plot(x,kn);二、 k近邻分类器(kNN)根据贝叶斯公式Pwix=pnxwiPwipnx，其中Pwi=nin、pnxwi

8、=kiniVn、pnx=knnVn，得出Pwix=kikn。也即，某点的类别和其邻域内相邻的几个点有关（k1+k2+kc=kn，n1+n2+nc=n）。这里用的是iris数据拿来测试，iris数据每类的前25个为训练数据，每类的后25个为测试数据。clear;% 加载数据load fisheriris;irisFeatures=meas;% 创建训练数据和测试数据irisTrain=irisFeatures(1:25,:);irisFeatures(51:75,:);irisFeatures(101:125,:);irisTrainNum=ones(25,1);ones(25,1)*2;one

9、s(25,1)*3;irisTest=irisFeatures(26:50,:);irisFeatures(76:100,:);irisFeatures(126:150,:);% 创建knn分类器，近邻点个数为15mdl=fitcknn(irisTrain,irisTrainNum,'NumNeighbors',15);% 预测测试点的类别irisTestNum=predict(mdl,irisTest);% 统计结果sum=0;for i=1:1:75 if(irisTestNum(i) = irisTrainNum(i) sum=sum+1; endenddisp(1-su

10、m/75.0);% 显示结果stem(irisTestNum);这里采用的是15近邻分类器，其中测试数据的准确率达到93.33%。下图为具体的结果。三、 MeanShift聚类MeanShift可以用与聚类、图像处理以及运动跟踪上，其核心公式为x=i=1nxiKx-xihi=1nKx-xih通过多次迭代，最终可以收敛到局部最优的位置，其中核函数K()为Kz=12exp-z22算法的大致流程为：第一步、从x*=xj出发第二步、 x=i=1nxiKx*-xihi=1nKx*-xih，若x-x*2<则中止迭代，否则x*=x继续第二步这里我们测试这个算法的性能，如下图为随机生成满足二维正态分布的

11、两组点。我们使用MeanShift算法进行聚类：如上图，可以看到聚类的效果还是比较理想的，这里设置的h=2。h太小或者太大会大大的影响结果，所以这里对于h的取值要适度。附上代码：clear;% 随机生成满足二维正态分布的两组点mu1=5 5;s1=1 0;0 1;n1=100;mu2=10 10;s2=2 0;0 1;n2=100;n=n1+n2;data=mvnrnd(mu1,s1,n1);mvnrnd(mu2,s2,n2);figure(1);plot(data(:,1),data(:,2),'*'); % meanshifth=2; %核大小sortNum=zeros(n

12、,2); %最终迭代点dif=0.01; %设置终止迭代的条件for j=1:1:n x=data(j,:); pre_x=0 0; while norm(pre_x-x)>dif pre_x=x; u=0 0; %分子 d=0; %分母 %公式 分子和分母 累加 for i=1:n ui=(x-data(i,:)/h; k=(1/sqrt(2*pi)*exp(-0.5*norm(ui).2); %高斯核 u=u+data(i,:)*k; d=d+k; end x=u/d; end sortNum(j,:)=x;end% 聚类label=zeros(n,1); %创建标签nowId=1;

13、 %最大类别idlabel(1)=nowId; %设置第一个点的标签为1% 计算所有点的标签for i=2:1:size(sortNum,1) % 将该点的最终迭代坐标与已知标签点的最终迭代坐标比较 % 这里设置两坐标之差的二范数小于2到3倍的dif 则聚为一类 % 否则为新的标签 flag=0; %用来判断是否找到同类的点 for j=1:1:i-1 if norm(sortNum(i,:)-sortNum(j,:)<dif*3 label(i)=label(j); flag=1; break; end end if(flag=0) nowId=nowId+1; label(i)=nowId;