多元正态分布 Matlab,概率区域
Multivariate Normal Distribution Matlab, probability area
我有 2 个数组:一个是 x 坐标,另一个是 y 坐标。
作为蒙特卡洛模拟的结果,两者都是正态分布。我知道如何找到两个数组的 sigma 和 mu,并获得 95% 的置信区间:
[mu,sigma]=normfit(x_array);
hist(x_array);
x=norminv([0.025 0.975],mu,sigma)
但是,两个数组相互关联。为了绘制组合数组的概率分布,我使用多元正态分布。在 MATLAB 中,这给了我:
[MuX,SigmaX]=normfit(x_array);
[MuY,SigmaY]=normfit(y_array);
mu = [MuX MuY];
Sigma=cov(x_array,y_array);
x1 = MuX-4*SigmaX:5:MuX+4*SigmaX; x2 = MuY-4*SigmaY:5:MuY+4*SigmaY;
[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],mu,Sigma);
F = reshape(F,length(x2),length(x1));
surf(x1,x2,F);
caxis([min(F(:))-.5*range(F(:)),max(F(:))]);
set(gca,'Ydir','reverse')
xlabel('x0-as'); ylabel('y0-as'); zlabel('Probability Density');
到目前为止一切顺利。现在我想计算 95% 概率区域。我正在寻找 mndinv
和 norminv
一样的函数。然而,MATLAB 中不存在这样的函数,这是有道理的,因为有无限的可能性......有人有关于如何获得 95% 概率区域的提示吗?提前致谢。
要获取顶部所在位置 F 的数值,您应该使用 top5=prctile(F(:),95)
。这将 return F
的值,它将底部 95% 的数据限制为顶部 5%。
然后您可以通过
获得前 5%
Ftop=zeros(size(F));
Ftop=F>top5;
Ftop=Ftop.*F;
%// optional: Ftop(Ftop==0)=NaN;
surf(x1,x2,Ftop,'LineStyle','none');
对于双变量情况,您可以添加面积对应于 NORMINV(95%) 的椭圆。这个椭圆是唯一标识的,为了证明,请参阅 link 中的第一个来源。
% Suppose you know the distribution params, or you got them from normfit()
mu = [3, 7];
sigma = [1, 2.5
2.5 9];
% X/Y values for plotting grid
x = linspace(mu(1)-3*sqrt(sigma(1)), mu(1)+3*sqrt(sigma(1)),100);
y = linspace(mu(2)-3*sqrt(sigma(end)), mu(2)+3*sqrt(sigma(end)),100);
% Z values
[X1,X2] = meshgrid(x,y);
Z = mvnpdf([X1(:) X2(:)],mu,sigma);
Z = reshape(Z,length(y),length(x));
% Plot
h = pcolor(x,y,Z);
set(h,'LineStyle','none')
hold on
% Add level set
alpha = 0.05;
r = sqrt(-2*log(alpha));
rho = sigma(2)/sqrt(sigma(1)*sigma(end));
M = [sqrt(sigma(1)) rho*sqrt(sigma(end))
0 sqrt(sigma(end)-sigma(end)*rho^2)];
theta = 0:0.1:2*pi;
f = bsxfun(@plus, r*[cos(theta)', sin(theta)']*M, mu);
plot(f(:,1), f(:,2),'--r')
来源
我有 2 个数组:一个是 x 坐标,另一个是 y 坐标。 作为蒙特卡洛模拟的结果,两者都是正态分布。我知道如何找到两个数组的 sigma 和 mu,并获得 95% 的置信区间:
[mu,sigma]=normfit(x_array);
hist(x_array);
x=norminv([0.025 0.975],mu,sigma)
但是,两个数组相互关联。为了绘制组合数组的概率分布,我使用多元正态分布。在 MATLAB 中,这给了我:
[MuX,SigmaX]=normfit(x_array);
[MuY,SigmaY]=normfit(y_array);
mu = [MuX MuY];
Sigma=cov(x_array,y_array);
x1 = MuX-4*SigmaX:5:MuX+4*SigmaX; x2 = MuY-4*SigmaY:5:MuY+4*SigmaY;
[X1,X2] = meshgrid(x1,x2);
F = mvnpdf([X1(:) X2(:)],mu,Sigma);
F = reshape(F,length(x2),length(x1));
surf(x1,x2,F);
caxis([min(F(:))-.5*range(F(:)),max(F(:))]);
set(gca,'Ydir','reverse')
xlabel('x0-as'); ylabel('y0-as'); zlabel('Probability Density');
到目前为止一切顺利。现在我想计算 95% 概率区域。我正在寻找 mndinv
和 norminv
一样的函数。然而,MATLAB 中不存在这样的函数,这是有道理的,因为有无限的可能性......有人有关于如何获得 95% 概率区域的提示吗?提前致谢。
要获取顶部所在位置 F 的数值,您应该使用 top5=prctile(F(:),95)
。这将 return F
的值,它将底部 95% 的数据限制为顶部 5%。
然后您可以通过
获得前 5%Ftop=zeros(size(F));
Ftop=F>top5;
Ftop=Ftop.*F;
%// optional: Ftop(Ftop==0)=NaN;
surf(x1,x2,Ftop,'LineStyle','none');
对于双变量情况,您可以添加面积对应于 NORMINV(95%) 的椭圆。这个椭圆是唯一标识的,为了证明,请参阅 link 中的第一个来源。
% Suppose you know the distribution params, or you got them from normfit()
mu = [3, 7];
sigma = [1, 2.5
2.5 9];
% X/Y values for plotting grid
x = linspace(mu(1)-3*sqrt(sigma(1)), mu(1)+3*sqrt(sigma(1)),100);
y = linspace(mu(2)-3*sqrt(sigma(end)), mu(2)+3*sqrt(sigma(end)),100);
% Z values
[X1,X2] = meshgrid(x,y);
Z = mvnpdf([X1(:) X2(:)],mu,sigma);
Z = reshape(Z,length(y),length(x));
% Plot
h = pcolor(x,y,Z);
set(h,'LineStyle','none')
hold on
% Add level set
alpha = 0.05;
r = sqrt(-2*log(alpha));
rho = sigma(2)/sqrt(sigma(1)*sigma(end));
M = [sqrt(sigma(1)) rho*sqrt(sigma(end))
0 sqrt(sigma(end)-sigma(end)*rho^2)];
theta = 0:0.1:2*pi;
f = bsxfun(@plus, r*[cos(theta)', sin(theta)']*M, mu);
plot(f(:,1), f(:,2),'--r')