• 2007-11-19

    正态分布和椭圆、椭球 - [工作]

    版权声明:转载时请以超链接形式标明文章原始出处和作者信息及本声明
    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);


    收藏到:Del.icio.us




    评论

  • 谢谢博主,这篇博文正是我需要的。我是一名统计专业的研究生,刚刚大体浏览了你的博客,知道你已经工作了,能冒昧地问一下你从事什么工作吗?因为我看到这篇博文属于“工作”类。