
function []=normal2d(mx,my,dx,dy,ro)
%
%  normal2d(mx,my,dx,dy,ro)
%
%  Plottar den simultana täthetsfunktionen för (X,Y) om
%
%          X är N(mx,dx)
%
%          Y är N(my,dy)
%
%          (X,Y) har korrelationskoeficient ro.
%
%  Betingade väntevärdesfunktionen E(Y|X=x) plottas tillsammans
%  med nivåkurvor för den simultana fördelningen.
%  
subplot(1,1,1)

if ro<=-1|ro>=1
  error('    korrelation ska uppfylla: -1 < ro < 1')
end
if dx<=0|dy<=0
  error('    standardavvikelse ska uppfylla: d > 0')
end
d=max(dx,dy);
a1=3;
a2=10;
xmin = mx-a1*dx;
xmax = mx+a1*dx;
ymin = my-a1*dy;
ymax = my+a1*dy;
bx = (xmax-xmin)/80;
by = (ymax-ymin)/30;

x = [xmin:dx/a2:xmax];
y = [ymin:dy/a2:ymax];
[X,Y] = meshgrid(x,y);

Qxy = (((X-mx)/dx).^2-2*ro*(X-mx).*(Y-my)/(dx*dy)+((Y-my)/dy).^2)/(1-ro^2);
kxy = 1/(2*pi*dx*dy*sqrt(1-ro^2));
fxy = kxy*exp(-Qxy/2);

Qx = ((x-mx)/dx).^2;
kx = 1/sqrt(2*pi)/dx;
fx = kx*exp(-Qx/2);

Qy = ((y-my)/dy).^2;
ky = 1/sqrt(2*pi)/dy;
fy = ky*exp(-Qy/2);

subplot(2,1,1)
 
  %S(:,1)=[zeros(1,63-15),0:(.25-0)/15:.25]';
  S(:,1)=zeros(64,1);
  S(:,2)=zeros(64,1);
  S(:,3)=zeros(64,1);
  S(:,1)=[.2:(1-.2)/63:1]';
  S(:,2)=[.2:(1-.2)/63:1]';
  %S(:,3)=[.08:(1-.08)/63:1]';
  colormap(S)
  brighten(.75)
  meshz(x,y,fxy')
  axis([xmin-bx xmax+bx ymin-by ymax+by 0 max(max(fxy))])
  view(-10,45)
  xlabel('x')
  ylabel('y')
  title('Simultan täthetsfunktion f(x,y)')

subplot(2,2,4)
  contour(x,y,fxy,'y')
  hold on
%  plot(x,my+(x-mx)*ro*dy/dx,'w')
  plot(x,my+(x-mx)*ro*dy/dx,'b')
  child=get(gca,'Children');
  set(child(1),'linewidth',1)
  grid
  axis([xmin-bx xmax+bx ymin-by ymax+by])
  hold off
  title('E(Y|X=x) och nivåkurvor för f(x,y)')
  xlabel('x')
  pos=get(gca,'position');
  set(gca,'position',[pos(1:3) pos(4)-0.08]);

subplot(2,2,3)
  plot(x,fx,'.');
barn=get(gca,'children');
set(barn,'markersize',8);
  hold on
  plot(y,fy);
  hold off
  axis([min(xmin,ymin) max(xmax,ymax) 0 max(max(fx,fy))*1.2])
  title('Marginella tätheter f(x) (..) och f(y) (-)')
  pos=get(gca,'position');
  set(gca,'position',[pos(1:3) pos(4)-0.08]);






