function [X,ml,cl,acc,ph] = snorm(th,n,method,s,x_start,n_init,th2,p);
% Generates random number from a bivariate normal distribution or from a
% mixture of two bivariate normals, with format
%
%  [X,ml,covloss,acc,ph]
%    = snorm(th,n,method,s,x_start,n_init,th2,p) 
%
% ARGUMENTS:
% th=[m1 m2 sig1 sig2 rho]    :  The parameter vector of a bivariate normal
%                                 with mean m = [m1 m2], variances
%                                 [sig1^2 sig^2] and correlation rho
% n                            :  number of bivariate vectors generated
% method                       :  Either of
%                                 'Ch' = Choleski
%                                 'Gi' = Gibbs
%                                 'MH' = Metropolis-Hastings,random walk kernel
% s                            :  Specifies the random walk kernel
%                                 N_2([0 0],s^2*I_2) for MH
% x_start                      :  Starting point of Gibbs and MH
%                                 (default = [0 0])
% n_init                       :  Number of iterations for Gibbs and MH
%                                 before generating output (default = 0)
% th2                          :  Parameter vector of second component in
%                                 mixture (default = one component)
% p                            :  mixture probability:
%                                 f = (1-p)*f_th + p*f_th2
% OUTPUTS:
% X                            :  Output sample = (n x 2)-matrix with
%                                 rows x_1,...,x_n  
% ml                           :  Gives the affine invariant loss function
%                                 ||V^(-1/2)(mean(X)-m)|| for mean, with
%                                 V=covariance (short for mean_loss)
% cl                           :  Gives the affine invariant loss function
%                                 ||V^(-1/2)(cov(X)-V)V^(-1/2)|| for covariance
%                                 matrix (short for cov_loss)
% acc                          :  Relative frequency of accepted moves
%                                 (only for MH)
% ph                           :  Proportion of time visiting the second
%                                 component (short for p^hat)
m1 = th(1);
m2 = th(2);
sig1 = th(3);
sig2 = th(4);
rho = th(5);
V(1,1) = sig1^2;
V(2,2) = sig2^2;
V(1,2) = rho*sig1*sig2;
V(2,1) = V(1,2);
[Q,D] = eig(V);
d = diag(D);
Dsqrt = diag(sqrt(d));
Vsqrt = Q*Dsqrt*Q';
m_vec = [m1 m2];
acc=0 ;
I1 = 0;  % This initialization is just to avoid warnings
I2 = 0;  % This initialization is just to avoid warnings

if nargin >= 7
  m21 = th2(1);
  m22 = th2(2);
  sig21 = th2(3);
  sig22 = th2(4);
  rho2 = th2(5);
  V2(1,1) = sig21^2;
  V2(2,2) = sig22^2;
  V2(1,2) = rho2*sig21*sig22;
  V2(2,1) = V2(1,2);
  [Q,D] = eig(V2);
  d = diag(D);
  Dsqrt = diag(sqrt(d));
  V2sqrt = Q*Dsqrt*Q';
  m2_vec = [m21 m22];
end;

if nargin <= 5, n_init = 0; end;
if (nargin <= 4 & ~strcmp(method,'Ch')), x_start=[0 0]; end;

if strcmp(method,'Ch')

  X = ones(n,1)*m_vec + randn(n,2)*Vsqrt;
  if nargin >= 7
    X2 = ones(n,1)*m2_vec + randn(n,2)*V2sqrt;
    I = (rand(n,1)<p*ones(n,1));
    ph = sum(I)/n;
    I = [I I];
    X = I.*X2 + (ones(n,2)-I).*X;
  end

elseif strcmp(method,'Gi')

  x = x_start;

  ph = 0;
  for i = 1:n_init+n

    x1 = normrnd(m1 + rho*sig1*(x(2)-m2)/sig2,sig1*sqrt(1-rho^2));
    if nargin >= 7
      x21 = normrnd(m21 +rho2*sig21*(x(2)-m22)/sig22,sig21*sqrt(1-rho2^2));
      px2 =p*normpdf(x(2),m22,sig22)/(p*normpdf(x(2),m22,sig22)+(1-p)*normpdf(x(2),m2,sig2));
      I2 = (rand < px2);
      x1 = x21*I2 + x1*(1-I2);
    end;

    x2 = normrnd(m2 + rho*sig2*(x1-m1)/sig1,sig2*sqrt(1-rho^2));
    if nargin >= 7
      x22 = normrnd(m22 +rho2*sig22*(x1-m21)/sig21,sig22*sqrt(1-rho2^2));
      px1 =p*normpdf(x1,m21,sig21)/(p*normpdf(x1,m21,sig21)+(1-p)*normpdf(x1,m1,sig1));
      I1 = (rand < px1);
      x2 = x22*I1 + x2*(1-I1);
    end;

    x = [x1 x2];

    if i>n_init
      X(i-n_init,:) = x;
      ph = ph + I1+I2;
    end

  end;
  ph = ph/(2*n);

elseif strcmp(method,'MH')

  x = x_start;
  acc = 0;
  if nargin >= 7, ph = 0; end;

  for i = 1:n_init+n

    y = x + s*[randn randn];
    qxy = bnormpdf(y,[0,0,s,s,0]);
    qyx = qxy;
    if nargin <= 6
      pix = bnormpdf(x,th);
      piy = bnormpdf(y,th);
    else
      pix = bnormpdf(x,th,th2,p);
      piy = bnormpdf(y,th,th2,p);
    end;
    prob = min(piy*qyx/(pix*qxy),1);
    I = (rand < prob);
    x = y*I + x*(1-I);

    if i>n_init
      X(i-n_init,:) = x;
      acc = acc + I;
      if nargin >= 7
        ph = ph + p*bnormpdf(x,th2)/((1-p)*bnormpdf(x,th)+p*bnormpdf(x,th2));
      end;
    end;

  end;

  acc = acc/n;
  if nargin >= 7, ph = ph/n; end;

else

  disp('This method is not implemented')

end;

if nargin >= 7
  V = p*V2 + (1-p)*V + p*(1-p)*(m_vec-m2_vec)'*(m_vec-m2_vec);
  m_vec = p*m2_vec + (1-p)*m_vec;
  [Q,D] = eig(V);
  d = diag(D);
  Dsqrt = diag(sqrt(d));
  Vsqrt = Q*Dsqrt*Q';
end;

ml = sqrt((mean(X)-m_vec)*inv(V)*(mean(X)-m_vec)');
Vsqrtinv = inv(Vsqrt);
cl = Vsqrtinv*(cov(X)-V)*Vsqrtinv;
cl = sqrt([1 1]*(cl.*cl)*[1;1]);


