R/checkFunctions.R

In eifer4/stochasticSampling: Stochastic partitioning sampler.

multiT <- function(mk,nk, h, h0, beta0, alpha0, nu, sigma02, N, X, Y){
  
  intercept <- matrix(diag(nk), nrow=N*nk, ncol=nk, byrow=TRUE)
  Xmat <- matrix(NA, nrow=N*nk, ncol=mk, byrow=TRUE)
  n<-0
  xidx <- 1
  while(n < N*nk){
    Xmat[(1+n):(nk+n),] <- rep(X[xidx,], each=nk)
    n <- n + nk
    xidx <- xidx + 1
  } 
  
  Ytemp <- matrix(t(Y), nrow=N*nk, ncol=1)
  
  W <- cbind(intercept, Xmat)
  H0 <- diag(c(h0 * rep(1,nk), h * rep(1,mk)))
  Sigma <- nu/sigma02 * (W %*% tcrossprod(H0 ,W) + diag(nk*N)) 
  cholSigma <- chol(Sigma)
  ldetSigma <- 2*sum(log(diag(cholSigma)))
  theta0 <- c(alpha0 * rep(1,nk), beta0*rep(1,mk))
  Yresid <- Ytemp - W %*% theta0
  
  ldens <- lgamma(sigma02 + N * nk * 0.5) - lgamma(sigma02) - 
    (nk * N * 0.5) * log(2 * pi * sigma02) -
    0.5 * ldetSigma - 
    (0.5 * nk * N + sigma02) * 
    log(1 + 1/(2*sigma02) * crossprod(Yresid, chol2inv(cholSigma) %*% Yresid))
  
  return(ldens)
}