Nothing
#####################################################################################
## Author: Daniel Sabanés Bové [daniel *.* sabanesbove *a*t* ifspm *.* uzh *.* ch]
## Project: hypergsplines
##
## Time-stamp: <[getSamples.R] by DSB Mit 25/04/2012 17:42 (CEST)>
##
## Description:
## Get posterior samples for a specific model configuration.
##
## History:
## 14/09/2010 file creation
## 15/09/2010 first complete version
## 27/09/2010 sort the spline coefficients samples according to the
## covariates they belong to
## 28/09/2010 get the samples.linearCoefs proper variable (row)names
## 29/09/2010 correct the precision matrix calculation: *inverse* of diagRho
## is added!
## 21/11/2010 adapt to the generalization dimSplineBasis of nKnots:
## with cubic O-Splines, nKnots is smaller than the dimension of the
## spline basis matrix.
## 11/01/2011 change format of linearCoefs to list with covariate names
## instead of matrix with row names (necessary for compatibility
## of output with that from getBmaSamples).
## 26/07/2011 start revision to allow for non-integer vector config
## 28/07/2011 use new interface for getRhos
## 29/03/2012 use new interface for rshrinkage
## 25/04/2012 fix bug in null model case: bStar was numeric(0) in that case
## before.
#####################################################################################
##' @include rshrinkage.R
##' @include invGamma.R
##' @include modelData.R
##' @include calculateModel.R
{}
##' Get posterior samples for a specific model configuration
##'
##' @param config the model configuration vector
##' @param nSamples number of samples to simulate
##' @param modelData the data necessary for model estimation, which
##' is the result from \code{\link{modelData}}
##' @return A list with samples from the shrinkage hyperparameter
##' t = g / (g + 1), the regression variance, and the (linear and spline)
##' coefficients.
##'
##' @example examples/getSamples.R
##'
##' @export
##' @keywords regression
##' @author Daniel Sabanes Bove \email{daniel.sabanesbove@@ifspm.uzh.ch}
getSamples <- function(config,
nSamples,
modelData)
{
## checks and extracts:
stopifnot(is.integer(nSamples),
nSamples > 0L,
identical(length(nSamples),
1L),
identical(length(config), modelData$nCovs))
config <- as.numeric(config)
## calculate the (marginal) model
marginal <- calculateModel(config=config,
modelData=modelData)
## sample the shrinkage hyperparameter t,
## yet only if this is not the null model:
samples.t <-
if(marginal$isNullModel)
numeric(0)
else
rshrinkage(n=nSamples,
marginal=marginal,
modelData=modelData)
## compute corresponding parameters of the inverse gamma posterior
aStar <- (modelData$nObs - 1) / 2
sst <- sum((marginal$y - marginal$meanY)^2)
tTimesR2 <-
if(marginal$isNullModel)
0
else
samples.t * marginal$coefR2
bStar <- (1 - tTimesR2) * sst / 2
## now sample the regression variance
samples.sigma2 <- rinvGamma(n=nSamples,
a=aStar,
b=bStar)
## the intercept
samples.intercept <- rnorm(n=nSamples,
mean=mean(modelData$y),
sd=sqrt(samples.sigma2 / modelData$nObs))
## the coefficients of linear effects
samples.linearCoefs <-
if(marginal$isNullModel)
list()
else
{
## cholesky decomposition of XtX (without intercept)
XtX <- crossprod(marginal$X.lin)
XtXroot <- chol(XtX)
## then sample the coefficients:
## first N(0, I_{dim.lin})
simCoefs <- matrix(data=rnorm(n=marginal$dim.lin * nSamples),
nrow=marginal$dim.lin,
ncol=nSamples)
## then N(0, t * sigma2 * I_{dim.lin})
simCoefs <- sweep(simCoefs,
MARGIN=2,
STATS=sqrt(samples.t * samples.sigma2),
FUN="*")
## then N(0, t * sigma2 * (XtX)^{-1})
simCoefs <- backsolve(r=XtXroot,
x=simCoefs,
k=marginal$dim.lin)
## and finally N(t * betaOLS, t * sigma2 * (XtX)^{-1})
simCoefs + sapply(samples.t,
FUN="*",
marginal$betaOLS)
}
## finally the spline coefficients samples
samples.splineCoefs <-
if(marginal$hasOnlyLinear)
list()
else
{
## number of columns of the spline design matrix:
dim.spline <- modelData$dimSplineBasis * length(marginal$whichSpline)
## diagonal rho matrix entries will come in here:
diagRho <- numeric(dim.spline)
## the whole spline design matrix will come in here:
splineDesign <- matrix(nrow=modelData$nObs,
ncol=dim.spline)
## start filling
nSplineTerms <- length(marginal$whichSpline)
splineNames <- character(nSplineTerms)
for(i in seq_len(nSplineTerms))
{
## which covariate?
s <- marginal$whichSpline[i]
## save the name for the list construction below
splineNames[i] <- names(modelData$Z.list[s])
## get the correct rho
d <- config[s]
thisRho <-
if(identical(d, floor(d)))
{
dNumber <- match(d - 1, modelData$splineDegrees)
modelData$rho.list[[s]][dNumber]
}
else
{
getRhos(lambdas=modelData$lambdas.list[[s]],
degrees=d - 1)
}
## positions of rho and spline basis columns
positions <- (i - 1) * modelData$dimSplineBasis + (1:modelData$dimSplineBasis)
## put rho into the diagonal (dimSplineBasis times)
diagRho[positions] <- thisRho
## put the spline basis into the design matrix
splineDesign[, positions] <- modelData$Z.list[[s]]
}
## so we can construct the precision matrix:
precMatrix <- crossprod(splineDesign)
diag(precMatrix) <- diag(precMatrix) + 1 / diagRho
## make a Cholesky decomposition of it
precRoot <- chol(precMatrix)
## then normal variate generation:
## first N(0, I_{dim.spline})
simCoefs <- matrix(data=rnorm(n=dim.spline * nSamples),
nrow=dim.spline,
ncol=nSamples)
## then N(0, sigma2 * I_{dim.spline})
simCoefs <- sweep(simCoefs,
MARGIN=2,
STATS=sqrt(samples.sigma2),
FUN="*")
## then N(0, sigma2 * precMatrix^{-1})
simCoefs <- backsolve(r=precRoot, x=simCoefs, k=dim.spline)
## compute the mean vector samples:
w <- modelData$y -
tcrossprod(rep(1, modelData$nObs), samples.intercept) -
modelData$X[, marginal$whichLinear, drop=FALSE] %*%
samples.linearCoefs
means <- forwardsolve(l=precRoot,
x=
crossprod(splineDesign,
w),
upper.tri=TRUE,
transpose=TRUE)
means <- backsolve(r=precRoot, x=means)
## and finally we have N(means, sigma2 * precMatrix^{-1})
## random vectors:
simCoefs <- simCoefs + means
## now sort these into a list corresponding to the
## covariates they belong to:
splineCoefs <- list()
for(i in seq_len(nSplineTerms))
{
## positions of spline coefficients
positions <- (i - 1) * modelData$dimSplineBasis + (1:modelData$dimSplineBasis)
## put the spline coefficients into the list
splineCoefs[[splineNames[i]]] <-
simCoefs[positions, , drop=FALSE]
}
## "return" the list
splineCoefs
}
if(! marginal$isNullModel)
{
## convert to list with covariates as names
samples.linearCoefs <-
split(samples.linearCoefs,
colnames(marginal$X.lin)[row(samples.linearCoefs)])
}
## return the list with all samples
return(list(t=samples.t,
sigma2=samples.sigma2,
intercept=samples.intercept,
linearCoefs=samples.linearCoefs,
splineCoefs=samples.splineCoefs))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.