Search code examples
matlabprobability-densityhistogram2d

histogram of Bivariate Normal distribution and pdf on Matlab

I have values from a Monte Carlo simulation. now I am supposed to plot the results using a bivariate Gaussian. Since it is experimental data on a specific topic, I was given the empirical formula of the PDF(x,y):

Formula

and this is the plot I should get( obviously it will be different depending on my data):

Plot

I don't understand how to implement the plot

This is the part of the code where I calculate STD, means:

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))

I'll probably have to use the histogram2 command, and I saw that in the properties there is a 'pdf' option. I would like to insert the pdf that I calculated, though. How can I do that?

Solution

  • To reiterate the comments, you’ll need to plot your projections separately from the histogram; this just isn’t a functionality that histogram2 has. Also, your Delta_x and Delta_y should be calculated with the point being plotted, not with the points in your dataset.

    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