R Code 9: Probit Model
N=200
x=rbinom(N,1,.75)
rtnpos=function(n,b)
{
u=runif(n)
b+qnorm(1-u*pnorm(b))
}
#generates n samples of a normal (b,1) truncated above zero
rtnneg=function(n,b)
{
u=runif(n)
b+qnorm(u*pnorm(-b))
}
M=2000
w=matrix(nrow=M,ncol=N)
z=1:M
sig2.0=1000
mu.0=0
#start values
z[1]=pnorm(.5)
w[1,]=rnorm(N)
above=x>0
A=sum(above)
for (m in 2:M)
{
w[m,above]=rtnpos(A,z[m-1])
w[m,!above]=rtnneg(N-A,z[m-1])
post.v=1/(N+1/sig2.0)
post.c=sum(w[m,])+(mu.0/sig2.0)
z[m]=rnorm(1,post.c*post.v,sqrt(post.v))
}
I=100 #People
J=50 #Trials Per Person
t.p = rbeta(I,5,2)
x=matrix(nrow=I,ncol=J)
for (i in 1:I) x[i,]=rbinom(J,1,t.p[i])
M=2000
w=matrix(nrow=M,ncol=N)
z=1:M
sig2.0=1000
mu.0=0
#start values
z[1]=pnorm(.5)
w[1,]=rnorm(N)
above=x>0
A=sum(above)
for (m in 2:M)
{
#recycle w across iterations
w[m,above]=rtnpos(A,z[m-1])
w[m,!above]=rtnneg(N-A,z[m-1])
post.v=1/(N+1/sig2.0)
post.c=sum(w[m,])+(mu.0/sig2.0)
z[m]=rnorm(1,post.c*post.v,sqrt(post.v))
}