function [betamat]= ...
    gs_MNP_I(X,y,k,M,n,r1,r2,mu0,V0,rg,E,sig);

R = r1+r2;
betamat=zeros(k,r2); 

Udraw=zeros(n*(M-1),1); %starting vector for latent utility
U=mat2cell(Udraw,(M-1)*ones(n,1),1);

%create differencing matrices to adjust the error variance
%note: in this simple case w/o RCs A,B, and C below will be the same, i.e.
%[1 2;2 1];
D=-eye(M);
D1=D;
D1(:,1)=1;
D1(1,:)=[];
D2=D;
D2(:,2)=1;
D2(2,:)=[];
D3=D;
D3(:,3)=ones(M,1);
D3(3,:)=[];

% Case 1: yi==1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=sig*D1*D1'; %2 by 2
V=inv(A); %see Robert(1995), p. 124
%variance for U1|[U2]
i=1;
Voth=V;
Voth(i,:)=[];
Voth(:,i)=[];
Vcov=V(:,i);
Vcov(i,:)=[];
Vself=V(i,i);
Aself=A(i,i);
Acov=A(:,i);
Acov(i,:)=[];
Aothinv=Voth-Vcov*Vcov'/Vself;
vA1_2=Aself-Acov'*Aothinv*Acov;
mA1_2=Acov'*Aothinv; %need this for conditional mean below
%variance for U2|[U1]
i=2;
Voth=V;
Voth(i,:)=[];
Voth(:,i)=[];
Vcov=V(:,i);
Vcov(i,:)=[];
Vself=V(i,i);
Aself=A(i,i);
Acov=A(:,i);
Acov(i,:)=[];
Aothinv=Voth-Vcov*Vcov'/Vself;
vA2_1=Aself-Acov'*Aothinv*Acov;
mA2_1=Acov'*Aothinv;

% Case 2: yi==2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=sig*D2*D2'; %2 by 2
V=inv(B); %see Robert(1995), p. 124
%variance for U1|[U2,U3]
i=1;
Voth=V;
Voth(i,:)=[];
Voth(:,i)=[];
Vcov=V(:,i);
Vcov(i,:)=[];
Vself=V(i,i);
Bself=B(i,i);
Bcov=B(:,i);
Bcov(i,:)=[];
Bothinv=Voth-Vcov*Vcov'/Vself;
vB1_2=Bself-Bcov'*Bothinv*Bcov;
mB1_2=Bcov'*Bothinv; %need this for conditional mean below
%variance for U2|[U1,U3]
i=2;
Voth=V;
Voth(i,:)=[];
Voth(:,i)=[];
Vcov=V(:,i);
Vcov(i,:)=[];
Vself=V(i,i);
Bself=B(i,i);
Bcov=B(:,i);
Bcov(i,:)=[];
Bothinv=Voth-Vcov*Vcov'/Vself;
vB2_1=Bself-Bcov'*Bothinv*Bcov;
mB2_1=Bcov'*Bothinv;

% Case 3: yi==3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C=sig*D3*D3'; %2 by 2
V=inv(C); %see Robert(1995), p. 124
%variance for U1|[U2,U3]
i=1;
Voth=V;
Voth(i,:)=[];
Voth(:,i)=[];
Vcov=V(:,i);
Vcov(i,:)=[];
Vself=V(i,i);
Cself=C(i,i);
Ccov=C(:,i);
Ccov(i,:)=[];
Cothinv=Voth-Vcov*Vcov'/Vself;
vC1_2=Cself-Ccov'*Cothinv*Ccov;
mC1_2=Ccov'*Cothinv; %need this for conditional mean below
%variance for U2|[U1,U3]
i=2;
Voth=V;
Voth(i,:)=[];
Voth(:,i)=[];
Vcov=V(:,i);
Vcov(i,:)=[];
Vself=V(i,i);
Cself=C(i,i);
Ccov=C(:,i);
Ccov(i,:)=[];
Cothinv=Voth-Vcov*Vcov'/Vself;
vC2_1=Cself-Ccov'*Cothinv*Ccov;
mC2_1=Ccov'*Cothinv;

thoucount=100;
% start main loop
for i=1:R
   %%%%%%%%%%%%%%%%%%%%%%
   % draw beta
   %%%%%%%%%%%%%%%%%%%%%
   % CAREFUL: differenced errors no longer have identity VCOV!
   % need to work with PERSON-specific D1D1', etc
   int1=0; 
   int2=0;
   j=1;
   for j=1:n
       Xi=X{j};
       Ui=U{j};
       Ei=E{j};
       int1=int1+Xi'*inv(Ei)*Xi;
       int2=int2+Xi'*inv(Ei)*Ui;
   end
   V1=inv(inv(V0)+int1);
   mu1=V1*(inv(V0)*mu0+int2);
   betadraw=mvnrnd(mu1,V1)';
   if i>r1
       betamat(:,i-r1)=betadraw;
   end
 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % draw latent utilities (actually utility-differences...)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  for j=1:n
      Xi=X{j};
      Ui=U{j};
      yi=y(j);
      ui=Xi*betadraw;%2 by 1 - already differenced!
      ui1=ui(1);
      ui2=ui(2);
      int2=ui2; %starting draw for second conditional
      if yi==1
          jj=1;
          for jj=1:rg
              m1=ui1 + mA1_2*(int2-ui2);
              s1=sqrt(vA1_2);
              int1=tnormrnd_robert_single(m1,s1,0);
              m2=ui2 + mA2_1*(int1-ui1);
              s2=sqrt(vA2_1);
              int2=tnormrnd_robert_single(m2,s2,0);
              jj=jj+1;
          end
          U{j}=[int1;int2];
      elseif yi==2
          jj=1;
          for jj=1:rg
              m1=ui1 + mB1_2*(int2-ui2);
              s1=sqrt(vB1_2);
              int1=tnormrnd_robert_single(m1,s1,0);
              m2=ui2 + mB2_1*(int1-ui1);
              s2=sqrt(vB2_1);
              int2=tnormrnd_robert_single(m2,s2,0);
              jj=jj+1;
          end
          U{j}=[int1;int2];
      else
          jj=1;
          for jj=1:rg
              m1=ui1 + mC1_2*(int2-ui2);
              s1=sqrt(vC1_2);
              int1=tnormrnd_robert_single(m1,s1,0);
              m2=ui2 + mC2_1*(int1-ui1);
              s2=sqrt(vC2_1);
              int2=tnormrnd_robert_single(m2,s2,0);
              jj=jj+1;
          end
          U{j}=[int1;int2];
      end
      
  end
          
  Udraw=cell2mat(U);
     
  if i== thoucount
      i
      thoucount=thoucount+100; 
  end
      
  
        
    
end


