function res =CSN_mcmc(Mu,Sigma,Gamma,Nu,Delta,num,burn)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Usage: [res, t_wi, t_inv] =CSN_mcmc(Mu,Sigma,Gamma,Nu,Delta,num,burn)
%   
%   Mu,Sigma,Gamma,Nu,Delta: parameters in CSN 
%   num: number of realisations
%   burn: burnin time, the first 'burn' realisations are not returned
%
%   res: returned realisations
%
%   notes:   see dahoiv.net/master
%
%   Ex of use res= CSN_mcmc( [2; 4],[5 0.7; 0.7 1],[5 0; 0 -5],[0;0],[1 0; 0 1], 400,20);
%   gives 400 realisations with mean around: [3.5392    3.3367]
%
% Created by Daniel Hoyer Iversen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  gamma_vt=Gamma*Sigma;
  sigma_v=Delta+Gamma*Sigma*Gamma';
  sigma_v_inv=inv(sigma_v);
  sigma_vt_chol=chol(Sigma-gamma_vt'*sigma_v_inv*gamma_vt,'lower');
  
  n_v=length(Nu);
  n_x=length(Mu);
  res=zeros(num,n_x);
  v=-Nu'*0+4;
  sigma_v__sigma_inv_temp=zeros(n_v-1,n_v);
  sd=zeros(n_v,1);
  j=0;
%  V=zeros(num+burn,n_v);
%  V(1,:)=v;
  while j < n_v
      j=j+1;
      index=1:n_v;
      index(j)=[];
      sigma_inv_temp=sigma_v_inv(index,index)-sigma_v_inv(index,j)* ...
          sigma_v_inv(index,j)'/sigma_v_inv(j,j);
      sigma_v__sigma_inv_temp(:,j)=sigma_v(index,j)'*sigma_inv_temp;
      sd(j)=sqrt(sigma_v(j,j)-sigma_v__sigma_inv_temp(:,j)'*sigma_v(index,j));
  end
  i=0;
  while i <num+burn
    i=i+1;
    j=0;
    while j< n_v
      j=j+1;
      index=1:n_v;
      index(j)=[];
      m=-Nu(j)+sigma_v__sigma_inv_temp(:,j)'*(v(index)'+Nu(index));
      v(j)=uni_tr(m,sd(j),0);
    end
%      V(i+1,:)=v;
     res(i,:)=Mu+gamma_vt'*sigma_v_inv*(v'+Nu) +sigma_vt_chol*randn(n_x,1);
  end

 % save 'V1' V
  
  
  res=res(burn+1:num+burn,:);
 %  mean(res)
%  toc
end




function res = uni_tr(m,sd,m_lower)
  m_lower=(m_lower-m)/sd;
  a=(m_lower+sqrt(m_lower^2+4))/2;
  z=-log(rand(1))/a+m_lower;
  if m_lower < a
    while rand(1) >exp(-(z - a)^2/2),
      z=-log(rand(1))/a+m_lower;
    end
  else
    'OBS'
    while rand(1) > exp(-(z - a)^2/2)*exp(-(m_lower-a)^2/2),
      z=exprnd(1/a)+m_lower;
    end
  end
  res=m+z*sd;
  
end