一、 实验目的1. 掌握密度函数监督参数估计方法;2. 掌握贝叶斯最小错误概率分类器设计方法。
二、 实验原理贝叶斯分类器是各种分类器中分类错误概率最小或者在预先给定代价的情况下平均风险最小的分类器。
它的设计方法是一种最基本的统计分类方法。
其分类原理是通过某对象的先验概率,利用贝叶斯公式计算出其后验概率,即该对象属于某一类的概率,选择具有最大后验概率的类作为该对象所属的类。
对于两类分类问题,已知先验概率P (ω1)和 P (ω2),以及类别标号 ω1和ω2,得到相应的类条件概率密度P (x |ω1), P (x|ω2), 由贝叶斯公式:计算得到条件概率P (ωi |x) (i=1,2),又称为后验概率。
如果:P (ωi |x)=max P (ωi |x),x ∈ ωi或者:P (ω1|x) > P (ω2|x),x ∈ ω1P (ω2|x) > P (ω1|x),x ∈ ω2三、 实验内容对于一个两类分类问题,设两类的先验概率相同(12()()P P ωω=),两类的类条件概率密度函数服从二维正态分布,即111(|)~(,)P N ωx μΣ 222(|)~(,)P N ωx μΣ其中,1[3,6]T =μ,10.5002⎡⎤=⎢⎥⎣⎦Σ,1[3,2]T =-μ,12002⎡⎤=⎢⎥⎣⎦Σ。
1.生成两类模式随机样本点并进行分类;2.设计最大似然估计算法对两类类条件概率密度函数进行估计;3.用2中估计的类条件概率密度函数设计最小错误概率贝叶斯分类器,实现对两类样本的分类。
四、实验步骤1.产生训练样本根据实验提供的先验均值向量和协方差矩阵,利用编写的multivrandn函数构造二维正态分布,分别产生N=500及N=1000个样本,所得结果如图1.1及1.2所示。
图1.1两类训练样本(N=500)图1.2两类训练样本(N=1000)2. 参数估计对产生的样本进行最大似然估计,估计出样本二维正态分布的均值向量和协方差矩阵。
其中,。
对于样本N=500估计结果如下:μ1=[3.0575 6.0294],μ2=[2.9404 −1.9881],∑1=[0.47340.00520.0052 2.1152],∑2=[2.1241−0.1233−0.1233 2.0153]对于样本N=1000估计结果如下:μ1=[3.0072 5.9923],μ2=[2.9869 −1.9336],∑1=[0.54990.04880.0488 2.0033],∑2=[1.98410.05250.0525 1.8102]3. 分类器设计根据上面得出的参数估计结果和贝叶斯最大后验概率判决准则设计分类器。
当,则。
设计分类函数,对样本进行分类判决。
例如对类别1中的第一个样本进行分类,结果如图2所示:图2.分类结果对两组样本进行分类,运用matlab 理论分别计算出N=500及N=1000个样本的分界线,结果如图3.1及3.2所示:图3.1两组样本分类结果(N=500)图3.2两组样本分类结果(N=1000)五、实验分析1.在产生样本的过程中,利用二维正态分布函数函数产生大量样本,经过均值和协方差矩阵的估计后可以看出:①随着样本数量N的增加,估计出的均值μ更接近于真实值;②方差相对变化较大,即样本数据的波动较大,不确定性越大。
故样本数越多,分类器将两类样本分离的会更加清楚,分类器的性能越好。
2.参数估计完全按照最大似然估计过程,结果如上所示,由于样本产生较好且数量较大,估计值也比较准确,从反面验证了参数估计过程的正确性。
3.根据最大后验概率判决准则,利用估计出的参数设计分类器,两组样本分类结果如图3.1及3.2所示,可以看出:①N=500和N=1000时均有个别误差,大部分样本分类正确;②随着样本数的增加,分类错误的样本数越少,即错误率越小,分类器的性能越好。
4.添加干扰,检测实验结果在产生样本时,添加均匀分布的一个干扰项,再次验证参数估计和分类结果如下:μ1=[8.5751 11.6213],μ2=[8.4160 3.7158],∑1=[7.4464−0.8548−0.85488.9503],∑2=[9.1159 1.24931.24938.8853]可以看出,得到的均值及方差估计值与真实值差距较大。
分类结果如图4所示:图4样本分类结果(有干扰)从分类结果可知,样本混淆现象严重,分类错误的样本数较多。
因此,在有干扰的情况下该分类器的训练误差较大,错误率高,性能较差。
六、程序代码本次实验程序代码共分为三部分:主程序及两个函数程序。
1.主程序如下:①N=500clear all;close all;clc;d=2; %二维pw1=0.5;pw2=0.5;u1=[3,6];u2=[3,-2]; %均值向量sigma1=[0.5,0;0,2]; %协方差矩阵sigma2=[2,0;0,2];N=500; %训练样本数samples1=multivrandn(u1,sigma1,N);samples2=multivrandn(u2,sigma2,N);figure(1);for i=1:Nplot(samples1(i,1),samples1(i,2),'b*');hold on;plot(samples2(i,1),samples2(i,2),'ro');hold on;endlegend('训练样本1','训练样本2');hold on;u_1=mean(samples1,1); %估计均值u_2=mean(samples2,1);sig1=zeros(2,2); %协方差矩阵的估计for i=1:Ntemp=(samples1(i,:)-u_1)'*(samples1(i,:)-u_1);sig1=sig1+temp;endsig1=sig1/N; %估计协方差矩阵sig2=zeros(2,2);for i=1:Ntemp=(samples2(i,:)-u_2)'*(samples2(i,:)-u_2);sig2=sig2+temp;endsig2=sig2/N;%% 画出实际分类线syms x1 x2f1=(1/(sqrt(2*pi).^d)./det(sig1).*exp(-1/2*([x1 x2]-u_1)*inv(sig1)*([x1 x2]-u_1)')); f2=(1/(sqrt(2*pi).^d)./det(sig2).*exp(-1/2*([x1 x2]-u_2)*inv(sig2)*([x1 x2]-u_2)')); f3=f2/f1-pw1/pw2;ezplot(f3,[-4,10]);hold on;text(-2,4,'分界线');②N=1000clear all;close all;clc;d=2; %二维pw1=0.5;pw2=0.5;u1=[3,6];u2=[3,-2]; %均值向量sigma1=[0.5,0;0,2]; %协方差矩阵sigma2=[2,0;0,2];N=1000; %训练样本数samples1=multivrandn(u1,sigma1,N);samples2=multivrandn(u2,sigma2,N);figure(1);for i=1:Nplot(samples1(i,1),samples1(i,2),'b*');hold on;plot(samples2(i,1),samples2(i,2),'ro');hold on;endlegend('训练样本1','训练样本2');hold on;u_1=mean(samples1,1); %估计均值u_2=mean(samples2,1);sig1=zeros(2,2); %协方差矩阵的估计for i=1:Ntemp=(samples1(i,:)-u_1)'*(samples1(i,:)-u_1);sig1=sig1+temp;endsig1=sig1/N; %估计协方差矩阵sig2=zeros(2,2);for i=1:Ntemp=(samples2(i,:)-u_2)'*(samples2(i,:)-u_2);sig2=sig2+temp;endsig2=sig2/N;%% 画出实际分类线syms x1 x2f1=(1/(sqrt(2*pi).^d)./det(sig1).*exp(-1/2*([x1 x2]-u_1)*inv(sig1)*([x1 x2]-u_1)'));f2=(1/(sqrt(2*pi).^d)./det(sig2).*exp(-1/2*([x1 x2]-u_2)*inv(sig2)*([x1 x2]-u_2)'));f3=f2/f1-pw1/pw2;ezplot(f3,[-4,10]);hold on;text(-2,4,'分界线');2.二维正态分布样本产生函数multivrandnfunction Y = multivrandn(u,R,M)% this function draws M samples from N(u,R)% where u is the mean vector(row) and R is the covariance matrix which must be positive definite n = length(u); % get the dimensionC = chol(R); % perform cholesky decomp R = C'CX = randn(M,n); % draw M samples from N(0,I)Z=unifrnd(1,10,M,n);Y = X*C +Z+ ones(M,1)*u;end3.分类判决函数分类fenlei( samples1(1,:),pw1,pw2,u_1,u_2,sig1,sig2,d )function fenlei(samples1(1,:),pw1,pw2,u_1,u_2,sig1,sig2,d )%UNTITLED2 Summary of this function goes here% Detailed explanation goes heref1=(1/(sqrt(2*pi).^d)./det(sig1).*exp(-1/2*(x-u_1)*inv(sig1)*(x-u_1)')); f2=(1/(sqrt(2*pi).^d)./det(sig2).*exp(-1/2*(x-u_2)*inv(sig2)*(x-u_2)')); if (f1/f2)>(pw2/pw1)disp('样本属于第一类');elsedisp('样本属于第二类');endend。