# R Code 7, Race Model

library('MCMCpack')
library('msm')
I=200

#SET TRUE VALUES
t.mu1=-.7
t.mu2=-.5
t.s2=.2

#MAKE "TRUE" LATENT DATA
t.z1=rnorm(I,t.mu1,sqrt(t.s2))
t.z2=rnorm(I,t.mu2,sqrt(t.s2))

#MAKE OBSERVABLE DATA
t=exp(pmin(t.z1,t.z2))
x=t.z1 < t.z2
r=log(t)

#######################################################
#ANALYSIS

M=2000
z1=matrix(nrow=M,ncol=I)

z2=z1
mu1=1:M
mu2=1:M
s2=1:M

#prior
mu.0=0
s2.0=100
a=1
b=1

#start values
mu1[1]=-.4
mu2[1]=-.4
s2[1]=1
z1[1,]=rep(0,I)
z2[1,]=rep(0,I)

m=2
for (m in 2:M)
{
z1[m,]=r
z2[m,]=r
z1[m,!x]=rtnorm(sum(!x),mu1[m-1],sqrt(s2[m-1]),r[!x],Inf)
z2[m,x]=rtnorm(sum(x),mu2[m-1],sqrt(s2[m-1]),r[x],Inf)
#sample mu1 given z1 and s2
v=1/((I/s2[m-1])+(1/s2.0))
c=(I*mean(z1[m,]/s2[m-1])+(mu.0/s2.0))
mu1[m]=rnorm(1,v*c,sqrt(v))

#sample mu2 given z2 and s2
c=(I*mean(z2[m,]/s2[m-1])+(mu.0/s2.0))
mu2[m]=rnorm(1,v*c,sqrt(v))

#sample s2 given everything else
s2.scale=sum((z1[m,]-mu1[m])^2+(z2[m,]-mu2[m])^2)/2+b
s2[m]=rinvgamma(1,shape=a+I,scale=s2.scale)
}