[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

# st: Gibbs Sampler

 From "Javier Escobal" To 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–2015 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index