samplePosNorm {hbmem}R Documentation

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 with gibbs.
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.

Author(s)

Michael S. Pratte

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

[Package hbmem version 0.2 Index]