function[out]=GwG(x,mu,A,bounds,rg)
%
% Gibbs-within-Gibbs routine to draw from an M-variate normal with one or
% more truncation constraints; M>1
%
% out:     M x 1 vector of draw from this density
% x:       M x 1 vector of current (most recent or starting) set of draws
% mu:      M x 1 vector of unconditional means
% V:       M x M unconditional VCOV
% bounds:  M x 2 matrix of lower and upper bounds for each element of x
% r:       number of reps for GwG (usually in the 20-50 range, speed
%          permitting)
%
% Based on Geweke (1991) and Robert (1995)
% Written by Klaus Moeltner, University of Nevada, Reno
% March 2010
T=length(x);
tdraw=x;

%Get conditional partitioned inverses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
V=inv(A); 
vcon=zeros(T,1);%Will collect conditional variances
mcon=zeros(T,T-1);%will collect vectors to pre-multiply conditional means
for gg=1:T;
    Voth=V;
    Voth(gg,:)=[];
    Voth(:,gg)=[];
    Vcov=V(:,gg);
    Vcov(gg,:)=[];
    Vself=V(gg,gg);
    Aself=A(gg,gg);
    Acov=A(:,gg);
    Acov(gg,:)=[];
    Aothinv=Voth-Vcov*Vcov'/Vself;
    vcon(gg)=Aself-Acov'*Aothinv*Acov;%1 by 1
    mcon(gg,:)=Acov'*Aothinv; %1 by T - 1
end

for ii=1:rg
    oth=tdraw;
    toth=mu;
          
    gg=1; %draw first element of x
          
          oth(gg)=[];%all other x's
          toth(gg)=[]; %all other unconditional means
          tself=mu(gg);%own unconditional mean
          m1=tself + mcon(gg,:)*(oth-toth);
          s1=sqrt(vcon(gg));
          bli=bounds(gg,1);
          bui=bounds(gg,2);
          if isinf(bli)==0 & isinf(bui)==0
              tdraw(gg)=tnormrnd(m1,s1,bli,bui); 
              %for double-bounded, this is more reliable than Robert
          elseif isinf(bli)==0 & isinf(bui)==1
              tdraw(gg)=tnormrnd_robert_single(m1,s1,bli);%lower truncation case
          elseif isinf(bli)==1 & isinf(bui)==0
              tdraw(gg)=-tnormrnd_robert_single(-m1,s1,-bui);%upper truncation case
          else
              tdraw(gg)=normrnd(m1,s1); %no truncation case
          end
          
          for gg=2:T
              oth=tdraw;
              toth=mu;
              oth(gg)=[];%all other thresholds
              toth(gg)=[]; %all other unconditional means
              tself=mu(gg);%own unconditional mean
              m1=tself + mcon(gg,:)*(oth-toth);
              s1=sqrt(vcon(gg));
              bli=bounds(gg,1);
              bui=bounds(gg,2);
              if isinf(bli)==0 & isinf(bui)==0
                  tdraw(gg)=tnormrnd(m1,s1,bli,bui); 
                  %for double-bounded, this is more reliable than Robert
              elseif isinf(bli)==0 & isinf(bui)==1
                  tdraw(gg)=tnormrnd_robert_single(m1,s1,bli);%lower truncation case
              elseif isinf(bli)==1 & isinf(bui)==0
                  tdraw(gg)=-tnormrnd_robert_single(-m1,s1,-bui);%upper truncation case
              else
                  tdraw(gg)=normrnd(m1,s1); %no truncation case
              end
          end
          ii=ii+1;
      end

out=tdraw;