%
%  math495/stat490 -- 02 october 03 -- djm
%
%  w05bivar.m:

figure(1);  clf
rand('state',6543210)

for j=1:4
    switch j
        %  initialize params (-1 < rho < +1)
        case{1}
            sx = 1.0;  sy = 1.0;  rho = 0.0;  alpha = asin(rho)/2;
            Nrvs = 1000;
        case{2}
            sx = 1.0;  sy = 1.0;  rho = 0.5;  alpha = asin(rho)/2;
            Nrvs = 1000;
        case{3}
            sx = 3.0;  sy = 1.0;  rho = 0.0;  alpha = asin(rho)/2;
            Nrvs = 1000;
        case{4}
            sx = 2.0;  sy = 1.0;  rho = 0.5;  alpha = asin(rho)/2;
            Nrvs = 1000;
    end

    %  transformation matrix
    Tmat = [ cos(alpha)/sx -sin(alpha)/sy ;
            -sin(alpha)/sx  cos(alpha)/sy ] / sqrt(1-rho^2);
    Imat = inv(Tmat);
    
    %  make Nrvs random pairs
    XYrvs = randn(2,Nrvs);
    BIrvs = Imat*XYrvs;

    %  scatterplot
    subplot(2,2,j)
    plot(BIrvs(1,:),BIrvs(2,:),'.')
    data = ['\bf (' num2str(sx) ',' num2str(sy) ',' num2str(rho) ')'];
    title(['\bf scatterplot of (\sigma_x,\sigma_y,\rho) = ' ...
           data])
    xlabel('\bf x-axis')
    ylabel('\bf y-axis')
    axis equal
    
    mean_xy = [mean(BIrvs(1,:)) mean(BIrvs(2,:))]
    cov_xy =   cov(BIrvs(1,:),BIrvs(2,:))
end