R/PCM.R
PCMVarAtTime.Rd
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)
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). |
A numeric k x k matrix
# 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