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))
}