二元正态分布的直方图和 Matlab 上的 pdf

Question

我有来自蒙特卡洛模拟的值。现在我应该使用双变量高斯绘制结果。由于是特定主题的实验数据，给出了PDF(x,y)的经验公式：

这是我应该得到的情节（显然它会根据我的数据而有所不同）：

看不懂剧情怎么实现

这是我计算 STD 的代码部分，意思是：

N=5000    
x=randn(N,1); %along-track direction [Nx1]
y=randn(N,1); %cross-track direction [Nx1]

mu_x=mean(x);
mu_y=mean(y);
Delta_x= x-mu_x;
Delta_y= y-mu_y;
rho_xy= 0; %termine di correlazione ipotizzato nullo

% STD calcolate dalla sola posizione geometrica dei punti di caduta
sigma_x= sqrt(1/N*sum((x-mu_x).^2)) 
sigma_y= sqrt(1/N*sum((y-mu_y).^2))

 %valutazione pdf
term1 = 1/((2 * pi* sigma_x * sigma_y )*(sqrt(1-rho_xy^2)));
term2 = -1/(2*(1-rho_xy^2));
term3 = (Delta_x / sigma_x) .^ 2;
term4 = (Delta_y / sigma_y) .^ 2;
term5 = -2 *rho_xy*Delta_x.*Delta_y/(sigma_x*sigma_y);
Pdf = term1 * exp(term2 * (term3 + term4 + term5))

我可能不得不使用

histogram2

命令，我看到在属性中有一个

'pdf'

选项。不过，我想插入我计算的 pdf。我该怎么做？

Answer 1

要重申评论，您需要将预测与直方图分开绘制；这不是

histogram2

拥有的功能。此外，您的

Delta_x

和

Delta_y

应该使用绘制的点进行计算，而不是数据集中的点。

N=5000;
x=randn(N,1); %along-track direction [Nx1]
y=randn(N,1); %cross-track direction [Nx1]

mu_x=mean(x);
mu_y=mean(y);
Delta_x= @(x)x-mu_x;
Delta_y= @(y)y-mu_y;
rho_xy= 0; %termine di correlazione ipotizzato nullo

% STD calcolate dalla sola posizione geometrica dei punti di caduta
sigma_x= sqrt(1/N*sum((x-mu_x).^2));
sigma_y= sqrt(1/N*sum((y-mu_y).^2));

 %valutazione pdf
term1 = 1/((2 * pi* sigma_x * sigma_y )*(sqrt(1-rho_xy^2)));
term2 = -1/(2*(1-rho_xy^2));
term3 = @(x)(Delta_x(x) / sigma_x) .^ 2;
term4 = @(y)(Delta_y(y) / sigma_y) .^ 2;
term5 = @(x,y)-2 *rho_xy*Delta_x(x).*Delta_y(y)/(sigma_x*sigma_y);
Pdf = @(x,y)term1 * exp(term2 * (term3(x) + term4(y) + term5(x,y)));

h=histogram2(x,y);

ax=h.Parent;
numPts=100;
xProjected = linspace(ax.XLim(1),ax.XLim(2), numPts);
yProjected = linspace(ax.YLim(1),ax.YLim(2), numPts);

zOverX = Pdf(xProjected,mu_y);
zOverX = zOverX/max(zOverX)*max(h.Values,[],'all');
zOverY = Pdf(mu_x,yProjected);
zOverY = zOverY/max(zOverY)*max(h.Values,[],'all');

hold on
plot3(xProjected,repmat(ax.YLim(2),1, numPts),zOverX)
plot3(repmat(ax.XLim(2),1, numPts),yProjected,zOverY)
hold off