hbmem (version 0.2)

samplePosNorm: Function samplePosNorm

Description

Samples posterior of mean parameters of the positive hierarchical linear normal model with a single Sigma2 $(x = N(exp(mu+alpha_i+beta_j),sigma2))$.

Usage

samplePosNorm(sample, y, sub, item, lag, I, J, R, 
    sig2mu, a, b, met, sigma2, sampLag)

Arguments

sample

Block of linear model parameters from previous iteration.

y

Vector of data

sub

Vector of subject index, starting at zero.

item

Vector of item index, starting at zero.

lag

Vector of lag index, zero-centered.

I

Number of subjects.

J

Number of items.

R

Total number of trials.

sig2mu

Prior variance on the grand mean mu; usually set to some large number.

a

Shape parameter of inverse gamma prior placed on effect variances.

b

Rate parameter of inverse gamma prior placed on effect variances. Setting both s2a AND s2b to be small (e.g., .01, .01) makes this an uninformative prior.

met

Vector of tuning parameter for metropolis-hastings sampling. There is one tuning parameter for mu, each of I alphas, each of J betas, s2alpha,s2beta,and theta. Those for s2alpha and s2beta are placeholders, as these parameters are sampled wi

sigma2

Variance of distribution.

sampLag

Logical. Whether or not to sample the lag effect.

Value

• The function returns a list. The first element of the list is the newly sampled block of parameters. The second element contains a vector of 0s and 1s indicating which of the decorrelating steps were accepted.

References

Not Published yet

See Also

hbmem

Examples

library(hbmem)

I=25
J=40
R=I*J
t.sigma2=3
t.mu=.5
t.sig2alpha=.2
t.sig2beta=.6
t.alpha=rnorm(I,0,sqrt(t.sig2alpha))
t.beta =rnorm(J,0,sqrt(t.sig2beta))
t.theta=-.5
sub=rep(0:(I-1),each=J)
item=rep(0:(J-1),I)
lag=scale(rnorm(R,0,sqrt(t.sigma2)/10))

tmean=1:R
for(r in 1:R) tmean[r]=exp(t.mu+t.alpha[sub[r]+1]+t.beta[item[r]+1]+t.theta*lag[r])
y=rnorm(R,tmean,sqrt(t.sigma2))

M=1000 #Way too low for real analysis!
B=I+J+4
block=matrix(0,nrow=M,ncol=B)
met=rep(.1,B);jump=.0001
b0=rep(0,B)
keep=500:M
for(m in 2:M)
{
  tmp= samplePosNorm(block[m-1,],y,sub,item,lag,I,J,R,100,.01,.01,met,t.sigma2,1)
  block[m,]=tmp[[1]]
  b0=b0+tmp[[2]]

  if(m<keep[1])
  {
   met=met+(b0/m<.3)*-jump +(b0/m>.5)*jump
   met[met<jump]=jump
  }
}

est=colMeans(block[keep,])
b0/M

par(mfrow=c(3,2))
plot(exp(block[keep,1]),t='l',main="Mu",ylab="Mu")
abline(h=exp(t.mu),col="blue")
abline(h=mean(y),col="green")
acf(block[keep,1],main="ACF of Mu")

est.alpha=est[2:(I+1)]
plot(est.alpha~t.alpha,ylab="Est. Alpha",xlab="True Alpha");abline(0,1)
est.beta=est[(I+2):(I+J+1)]
plot(est.beta~t.beta,ylab="Est. Beta",xlab="True Beta");abline(0,1)

plot(block[keep,(I+J+4)],t='l',main="Theta",ylab="Theta")
abline(h=t.theta,col="blue")

plot(density(block[keep,(I+J+2)]),col="red",main="Posterior of Variances",xlim=c(0,1))
abline(v=t.sig2alpha,col="red")
lines(density(block[keep,(I+J+3)]),col="blue")
abline(v=t.sig2beta,col="blue")