function margGausDemo(mu,sigma)
% margGausDemo(mu,sigma)
% this function will draw a few ellipses corresponding to a bivariate
% Gaussian with mean mu and covariance matrix sigma; it will then allow
% the user to interactively set a direction (by specifying two points
% othrough which a line would pass) and plot on the right the marginal
% pdf of a projection to that direction.
% 
% To stop, click to the left of the 2D Gaussian subplot.
% 
% If the arguments (sigma, or both) are ommitted, a random ones are
% used. In particular, a random covariance matrix is drawn using
% makeRandomSigma2d.m

% fill in missing arguments if needed
if (nargin < 1)
  mu=zeros(2,1);
end
if (nargin < 2)
  sigma = makeRandomSigma2d(.5,2,5)
end

% range of the axes should fit the covariance spread
range=3*det(sigma);
stepsize=(2*range)/300;
coordx=mu(1)+(-range:stepsize:range);
coordy=mu(2)+(-range:stepsize:range);

% prepare for plotting
figure(1);clf;

subplot(1,2,2);hold on;
title('Marginal $p(z), z=\mathbf{v}^T\mathbf{x}$','FontSize',24,'Interpreter','latex');

subplot(1,2,1);hold on;
title('Bivariate $p(\mathbf{x}$)','FontSize',24,'Interpreter','latex');

% draw some points from the Gaussian...
x=mvnormrnd(mu,sigma,100)';
plot(x(1,:),x(2,:),'k.','MarkerSize',8);

% draw ellipses for 90/50/10 % of max pdf
maxpdf=mvnormpdf(zeros(2,1),mu,sigma);
contlevels=[.9 .5 .1]*maxpdf;
[gridx,gridy]=meshgrid(coordx,coordy);
xy=[reshape(gridx,1,[]); reshape(gridy,1,[])];
p=mvnormpdf(xy,mu,sigma);
psurf=reshape(p,size(gridx));
[c,hc]=contour(coordx,coordy,psurf,contlevels);

% set the axes right
axis equal
axis([mu(1)-range mu(1)+range mu(2)-range mu(2)+range])

% If you make the cols array longer (e.g., try 'rbg') more than a single
% line will remain visible. 
cols='r';

hl={};for c=1:length(cols),hl{c}=[];end
hp={};for c=1:length(cols),hp{c}=[];end
hpts=[];
domore='y';

counter=0; % some bookkeeping in case we want to maintain multiple lines 

while (domore=='y') 
  % enter two points to specify a line
  [px(1),py(1)]=ginput(1);
  if (px(1) < -range)
    break;
  end
  [px(2),py(2)]=ginput(1);
  
  dy=diff(py);dx=diff(px);

  counter=mod(counter,length(cols))+1;
  if (~isempty(hl{counter})), delete(hl{counter}); end

  % calculate slope and draw the line
  slope=dy/dx;
  intercept=mu(2)-slope*mu(1);
  hl{counter}=line([mu(1)-range mu(1)+range],slope*[mu(1)-range,mu(1)+range]+intercept);
  set(hl{counter},'Color',cols(counter),'LineWidth',4);
  
  % unit vector in the direction of 
  v=[dx;dy]; v=v/norm(v);
  % the function in this case is v'*x
  asigma=v'*sigma*v;

  xprime=v'*x-v'*mu;
  
  subplot(1,2,2);
  if (~isempty(hp{counter})), delete(hp{counter}); end
  % set the range for the projections

  vspan=max(xprime)-min(xprime);
  vstep=(vspan*1.5)/300;
  vrange=mean(xprime)-150*vstep:vstep:mean(xprime)+150*vstep;
  hp{counter}=plot(vrange,normpdf(vrange,0,sqrt(asigma)),cols(counter),'LineWidth',4);

  if (~isempty(hpts)), delete(hpts), end
  hpts=plot(xprime(1,:),-.2,'ko','MarkerSize',6,'LineWidth',2);
  axis([-range range -.3 1]);
  
  %domore=input('more?','s');
  
end