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/
