add visualization of bivariate normal distribution

content/page/about.md

@@ -33,7 +33,7 @@ Session 1.1: Overview. Filtering and smoothing.
 
 Session 1.2: Iterated maps. Ordinary differential equations (ODEs). Discretization of ODEs.
 
-Session 1.3: Error propagation. Gaussain distribution.
+Session 1.3: Error propagation. Gaussian distribution.
 
 Session 1.4: Stochastic modelling. Statistical dependence and causality. 
 

mfiles/s_mvnormal.m

+pkg load statistics
+
+%This script illustrates a multivariate Gaussian distribution and its
+%marginal distributions
+
+%2-d Mean and covariance matrix
+vx = 1; vy = 1; vxy = -0.5;
+MeanVec = [0 0];
+CovMatrix = [vx vxy; vxy vy];
+
+%Get the sigma ellipses by transform a circle by the cholesky decomp
+L = chol(CovMatrix,'lower');
+t = linspace(0,2*pi,100); %Our ellipse will have 100 points on it
+C = [cos(t) ; sin(t)]; %A unit circle
+E1 = 1*L*C; E2 = 2*L*C; E3 = 3*L*C; %Get the 1,2, and 3-sigma ellipses
+
+%Define limits of plotting
+s_max = sqrt(max(sumsq(5*L*C,1)))
+X = linspace(-s_max, s_max, 30);
+Y = linspace(-s_max, s_max, 30);
+
+%Get the 1-d PDFs for the "walls"
+Z_x = normpdf(X,MeanVec(1), sqrt(CovMatrix(1,1)));
+Z_y = normpdf(Y,MeanVec(2), sqrt(CovMatrix(2,2)));
+
+%Get the 2-d samples for the "floor"
+Samples = mvnrnd(MeanVec, CovMatrix, 10000);
+
+figure; hold on; 
+%Plot the samples on the "floor"
+plot3(Samples(:,1),Samples(:,2),zeros(size(Samples,1),1),'k.','MarkerSize',2)
+%Plot the 1,2, and 3-sigma ellipses slightly above the floor
+%plot3(E1(1,:), E1(2,:), 1e-3+zeros(1,size(E1,2)),'Color','g','LineWidth',2);
+%plot3(E2(1,:), E2(2,:), 1e-3+zeros(1,size(E2,2)),'Color','g','LineWidth',2);
+%plot3(E3(1,:), E3(2,:), 1e-3+zeros(1,size(E3,2)),'Color','g','LineWidth',2);
+
+%Plot the histograms on the walls from the data in the middle
+%[n_x, xout] = hist(Samples(:,1),20);%Creates 20 bars
+%n_x = n_x ./ ( sum(n_x) *(xout(2)-xout(1)));%Normalizes to be a pdf
+%[~,~,~,x_Pos,x_Height] = makebars(xout,n_x);%Creates the bar points
+%plot3(x_Pos, Y(end)*ones(size(x_Pos)),x_Height,'-k')
+
+%Now plot the other histograms on the wall
+%[n_y, yout] = hist(Samples(:,2),20);
+%n_y = n_y ./ ( sum(n_y) *(yout(2)-yout(1)));
+%[~,~,~,y_Pos,y_Height] = makebars(yout,n_y);
+%plot3(X(1)*ones(size(y_Pos)),y_Pos, y_Height,'-k')
+%stem3(X(1)*ones(size(yout)),yout, n_y,'-k')
+
+%Now plot the 1-d pdfs over the histograms
+plot3(X, ones(size(X))*Y(end), Z_x,'-b','LineWidth',2); 
+plot3(ones(size(Y))*X(1), Y, Z_y,'-r','LineWidth',2);
+
+# Plot the pdf surface
+[xx, yy] = meshgrid(X, Y);
+rr = [xx(:) yy(:)];
+dr = rr - MeanVec;
+iCov = cholinv(CovMatrix);
+zz_pdf = exp(-0.5 * sum((dr * iCov) .* dr, 2)) / 2 / pi / sqrt(det(CovMatrix));
+zz = reshape(zz_pdf, size(xx));
+contour3(xx, yy, zz, 25, 'g');
+
+%Make the figure look nice
+grid on; view(45,55);
+axis([X(1) X(end) Y(1) Y(end)])