GibbsNormalGAU: LGP

In JonathanBradley28/CM: The Conjugate Multivariate Distribution

Description

This code implements a standard Gibbs sampler for normal data using a mixed effects model

Usage

GibbsNormalGAU(
  B,
  X,
  G,
  Z,
  etas2 = NA,
  betas2 = NA,
  xi = NA,
  taus2 = NA,
  Hphib2 = NA,
  sig2xi = NA,
  sig2beta2 = NA,
  printevery = 100
)

Arguments

B	The number of iterations of the Gibbs sampler
X	An nxp matrix of covariates.
G	An nxr matrix of basis functions.
Z	An n dimensional vector data.
etas2	An optional argument to define the initialized value for the correlated random effects.
betas2	An optional argument to define the initialized value for the fixed effects.
xi	An optional argument to define the initialized value for the uncorrelated random effects.
taus2	An optional argument to define the initialized value for the variance of the data.
Hphib2	An optional argument to define the initialized value for the covariance of the random effects.
sig2xi	An optional argument to define the initialized value for the variance of the random effects.
printevery	Option to print the iteration number of the MCMC.

Value

ans A list of updated parameter values.

Examples

#load the necessary packages
library(CM)

set.seed(123)
#define a test function
#A non-linear test function
lambda <- function(t) exp(1.1 + sin(2*pi*t))

#define some 1-d locations
points = seq(0,1,length.out=1001)
points=points[2:1001]
m = dim(as.matrix(points))[1]

#get the true mean at these locations
truemean<-matrix(0,m,1)
for (j in 1:length(points)){
 truemean[j,1] = lambda(points[j])
}

#simulate the data
data = matrix(0,m,1)
for (i in 1:m){
 data[i] = rnorm(1,truemean[i],1)
}

#plot the data
plot(data,xlab="time",ylab="Poisson counts",main="Response vs. time")

#covariate intercept-only
X = matrix(1,m,1)
p <- dim(X)[2]

##compute the basis function matrix
#compute thin-plate splines
r = 8
knots = seq(0,1,length.out=r)

#orthogonalize G
G = THINSp(as.matrix(points,m,1),as.matrix(knots,r,1))
outG<-qr(G)
G<-qr.Q(outG)

#orthogonalize X
outX<-qr(X)
X<-qr.Q(outX)

#Run the MCMC algorithm
output<-GibbsNormalGAU(2000,X,G,Z=data)

#trace plots (without burnin)
plot(as.mcmc(output$betas[1000:2000]))
plot(as.mcmc(output$etas[1,1000:2000]))
plot(as.mcmc(output$etas[8,1000:2000]))
plot(as.mcmc(output$xis[10,1000:2000]))
#estimates (remove a burnin)
f_est = apply(X%*%output$betas[,1000:2000]+G%*%output$etas[,1000:2000],1,mean)
f_lower= apply(X%*%output$betas[,1000:2000]+G%*%output$etas[,1000:2000],1,quantile,0.025)
f_upper= apply(X%*%output$betas[,1000:2000]+G%*%output$etas[,1000:2000],1,quantile,0.975)

#plot estimates and truth
plot(1:m,truemean,ylim = c(0,max(f_upper)+1))
lines(1:m,f_est,col="red")
lines(1:m,f_lower,col="blue")
lines(1:m,f_upper,col="blue")

JonathanBradley28/CM documentation built on July 8, 2021, 9:59 a.m.