高斯分布

http://liuxqsmile.blogbus.com/logs/17285781.html

正态分布和椭圆、椭球

二元正态分布（高斯分布）的等概率曲线是一个椭圆，而三元正态分布的等概率曲面是一个椭球。本文的目的是在已知概率分布参数（均值向量μ和协方差阵Σ）的情况下，在MATLAB中绘出等概率椭圆或椭球。

多元正态分布（multivariate normal distribution ）的概率密度：





很显然，其等概率密度曲面方程为



协方差矩阵作特征值分解：[V,D] = eig(Σ);使得






，即


因为

,所以

。

代入等概率方程，可以得到



作变量代换，

，得到

;

二元分布的情况下，y为包含两个元素的列向量，

。D为对角阵，两个对角线元素分布为d1，d2；则上式可写为标准的椭圆方程


。

实际上，我们对原始椭圆作一个平移变换（x-u）和一个旋转变换（V*），可以使其中心平移到原点，长短轴的方向与坐标轴重合。而绘等概率椭圆的过程则与此相反，先用标准的椭圆方程产生组成曲线的离散点，然后经过相反的旋转变换和平移，得到原始的椭圆。

举例来说，首先，设置二维正态分布的参数，均值和协方差，并用mvnrnd产生一组符合此分布的随机数。
Mu = [2 3]';

Sigma = [0.9 0.4;0.4 0.2];

p = mvnrnd(Mu,Sigma,100);

plot(p(:,1),p(:,2),'.','MarkerSize',6)

设置半径，进行特征值分解
r =1;

[V,D] = eig(Sigma);

用linspace产生一个坐标轴（y）上的一组等间隔离散坐标值，再根据标准椭圆方程产生对应的x的坐标。
y = linspace(-sqrt(r^2*D(2,2)),sqrt(r^2*D(2,2)),60);

% compute x

x(1,:) = sqrt((r^2-y(:).^2/D(2,2))*D(1,1));

x(1,:) = real(x(1,:));

这只产生了半个椭圆，还要产生另一半（注意两条曲线的坐标旋转方向要一致），然后旋转，平移，画图：
Ellip = [x,-x(1,:)]; % x

Ellip(2,:) = [y,fliplr(y)]; %y

Ellip = Ellip'*inv(V); % rotate

Ellip(:,1) = Ellip(:,1)+Mu(1); %shift

Ellip(:,2) = Ellip(:,2)+Mu(2);

hold on;

plot(Ellip(:,1),Ellip(:,2));

plot(Mu(1),Mu(2),'+'); %Plot center

最终的效果：



另外，在对原始椭圆做旋转变换时，如果在V的前面再乘以一项，改为

，则椭圆会变为圆。对多元正态分布的随机变量应用此变换，则其分布在各个方向上也变为均匀的。这就是信号处理中的白化变换。

三维正态分布的等概率曲面为椭球，其绘制过程也是类似的。

产生1/8曲面
xhalf = linspace(sqrt(r^2*D(1,1)),0,Nint);

Ninthalf = round(Nint/2);

zsect = zeros(Nint,Ninthalf);

ysect = zeros(Nint,Ninthalf);

for ti = 1:Nint

r2d = r^2 - xhalf(ti).^2/D(1,1);

ysect(ti,:) = linspace(0,sqrt(r2d*D(2,2)),Ninthalf);

zsect(ti,:) = sqrt((r2d - ysect(ti,:).^2/D(2,2) )*D(3,3));

xsect(ti,1:Ninthalf) = xhalf(ti);

end

zsect = real(zsect);



通过镜像产生1/4
%x>0,Z>0

xsect = [xsect,xsect];

ysect = [ysect,fliplr(ysect)];

zsect = [zsect,-fliplr(zsect)];



1/2
%x>0

xsect = [xsect,xsect];

ysect = [ysect,-fliplr(ysect)];

zsect = [zsect,fliplr(zsect)];



1/1
% make it a whole

xsect = [xsect;-flipdim(xsect,1)];

ysect = [ysect;flipdim(ysect,1)];

zsect = [zsect;flipdim(zsect,1)];



旋转、平移，完成。
% rotate

[lr,lc] = size(xsect);

for ti = 1:lr

for tj = 1:lc

newcodi = [xsect(ti,tj),ysect(ti,tj),zsect(ti,tj)]*inv(V);

xsect(ti,tj) = newcodi(1);

ysect(ti,tj) = newcodi(2);

zsect(ti,tj) = newcodi(3);

end

end

% shift

xsect = xsect+Mu(1);

ysect = ysect+Mu(2);

zsect = zsect+Mu(3);

surf(xsect,ysect,zsect);