Calculate the variance covariance k x k matrix at time t, under a PCM model

PCMVarAtTime(t, model, W0 = matrix(0, PCMNumTraits(model),
  PCMNumTraits(model)), regime = 1L, verbose = FALSE)

Arguments

t

positive numeric denoting time

model

a PCM model object

W0

a numeric matrix denoting the initial k x k variance covariance matrix at the root (default is the k x k zero matrix).

regime

an integer or a character denoting the regime in model for which to do the calculation; (Defaults to 1L meaning the first regime in the model)

verbose

a logical indicating if (debug) messages should be written on the console (Defaults to FALSE).

Value

A numeric k x k matrix

Examples

# a Brownian motion model with one regime modelBM <- PCM(model = "BM", k = 2) # print the model modelBM
#> Brownian motion model #> S3 class: BM, GaussianPCM, PCM; k=2; p=8; regimes: 1. Parameters/sub-models: #> X0 (VectorParameter, _Global, numeric; trait values at the root): #> [1] 0 0 #> Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; Choleski factor of the unit-time variance rate): #> , , 1 #> #> [,1] [,2] #> [1,] 0 0 #> [2,] 0 0 #> #> Sigmae_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; Choleski factor of the non-heritable variance or the variance of the measurement error): #> , , 1 #> #> [,1] [,2] #> [1,] 0 0 #> [2,] 0 0 #> #>
# assign the model parameters at random: this will use uniform distribution # with boundaries specified by PCMParamLowerLimit and PCMParamUpperLimit # We do this in two steps: # 1. First we generate a random vector. Note the length of the vector equals PCMParamCount(modelBM) randomParams <- PCMParamRandomVecParams(modelBM, PCMNumTraits(modelBM), PCMNumRegimes(modelBM)) randomParams
#> [1] 7.087678 -5.737014 2.103107 -9.209586 9.447748 2.449280 5.622451 #> [8] 2.882372
# 2. Then we load this random vector into the model. PCMParamLoadOrStore(modelBM, randomParams, 0, PCMNumTraits(modelBM), PCMNumRegimes(modelBM), TRUE)
#> [1] 8
# PCMVarAtTime(1, modelBM) # note that the variance at time 0 is not the 0 matrix because the model has a non-zero # environmental deviation PCMVarAtTime(0, modelBM)
#> [,1] [,2] #> [1,] 37.61093 16.205995 #> [2,] 16.20599 8.308067