Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: Gibbs Sampler


From   "Javier Escobal" <jescobal@grade.org.pe>
To   <statalist@hsphsun2.harvard.edu>
Subject   st: Gibbs Sampler
Date   Mon, 1 Dec 2003 20:30:19 -0500

Hello All

I would like to know if somebody has programmed a Gibbs Sampler.

I am interested in doing some iterative simulation. My idea is to reproduce
in STATA the program that I show below (it runs well in MATLAB). However my
programming skills are not good enough to translate it into stata.

Any help in this matter will be greatly appreciated

Javier

----MATLAB PROGRAM-----

%FICTIOUS DATA
n=100;
x=normrnd(0,1,n,k);
b=ones(k,m);
y=x*b+normrnd(0,1,n,m);
i=find(y(:,2)<=0);
j=find(y(:,2)>=0);
ni=length(i);
nj=length(j);

%SAMPLES
cens=1;
reps=2;
gibbs=2000;

%CENSORED
for c=1:cens
   
   %START
	z=y;
	b=inv(x'*x)*x'*z;
   t=0;
  
   %BURNIN AND GIBBS
   for r=1:reps
      for g=1:gibbs
         if c==1
            t=0;
         else
            t=unifrnd(max(z(i,2)),min(z(j,2)));
         end
 
u=mvnrnd(zeros(m,1),inv((z-x*inv(x'*x)*x'*z)'*(z-x*inv(x'*x)*x'*z)),n-k+m+1)
';
   		v=inv(u*u');
         s=v./v(1,1);
 
b=reshape(mvnrnd(reshape(inv(x'*x)*x'*z,k*m,1),kron(s,inv(x'*x)),1)',k,m);
         m1=x*b(:,1)+s(1,2)*inv(s(2,2))*(z(:,2)-x*b(:,2));
         s1=sqrt(s(1,1)-s(1,2)*inv(s(2,2))*s(2,1));
         m2=x*b(:,2)+s(2,1)*inv(s(1,1))*(z(:,1)-x*b(:,1));
         s2=sqrt(s(2,2)-s(2,1)*inv(s(1,1))*s(1,2));
 
z(i,1)=norminv(unifrnd(0,1,ni,1).*normcdf(ones(ni,1)*t,m1(i,:),s1),m1(i,:),s
1);
 
z(j,1)=norminv(unifrnd(0,1,nj,1).*(ones(nj,1)-normcdf(ones(nj,1)*t,m1(j,:),s
1))
                +normcdf(ones(nj,1)*t,m1(j,:),s1),m1(j,:),s1);
 
z(i,2)=norminv(unifrnd(0,1,ni,1).*normcdf(ones(ni,1)*t,m2(i,:),s2),m2(i,:),s
2);
         thetas(g,:)=t;
      	sigmas(g,:)=reshape(s,m*m,1)';
         betas(g,:)=reshape(b,k*m,1)';
         z1s(g,:)=z(:,1)';
         z2s(g,:)=z(i,1)';
         rg=[r g]
      end
   end
   
%RESULTS
	figure(1)
	plot(thetas)
	title('thetas')  


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index