-
Notifications
You must be signed in to change notification settings - Fork 3
/
GuessInitVecParams.Rd
144 lines (117 loc) · 4.78 KB
/
GuessInitVecParams.Rd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/GuessInitVecParams.R
\name{GuessInitVecParams}
\alias{GuessInitVecParams}
\title{A heuristic-based guess of the optimal parameters of a model}
\usage{
GuessInitVecParams(o, regimes = as.character(PCMRegimes(o)),
oParent = o, accessExpr = "", n = 1L,
argsPCMParamLowerLimit = NULL, argsPCMParamUpperLimit = NULL,
X = NULL, tree = NULL, SE = NULL, tableAnc = NULL,
varyParams = FALSE, returnWithinBoundsOnly = TRUE,
res = PCMParamRandomVecParams(o = oParent, k = PCMNumTraits(oParent), R
= PCMNumRegimes(oParent), n = n, argsPCMParamLowerLimit =
argsPCMParamLowerLimit, argsPCMParamUpperLimit = argsPCMParamUpperLimit),
...)
}
\arguments{
\item{o}{a PCM model object.}
\item{n}{integer denoting the number of parameter vectors to generate.}
\item{argsPCMParamLowerLimit, argsPCMParamUpperLimit}{lists of parameters
passed to \code{PCMParamLowerLimit} and \code{PCMParamUpperLimit}
respectively. (see the argument \code{returnWithinBoundsOnly}).}
\item{tableAnc}{an ancestor table for tree. Default \code{NULL}.}
\item{varyParams}{logical indicating if fixed (FALSE) or varying (TRUE)
parameter values should be returned for each of the \code{n} vectors. The
default value is TRUE.}
\item{returnWithinBoundsOnly}{logical indication if the returned parameter
vectors should be within the lower and upper bound of the model.}
\item{res}{an n x p matrix where p is the number of parameters in o.
Internal use only.}
\item{...}{additional arguments passed to implementations. Currently an
N x N matrix argument named treeVCVMat can be passed which is equal to the
phylogenetic covariance matrix for tree.}
\item{k}{integer denoting the number of traits in o.
Default: \code{PCMNumTraits(o)}.}
}
\value{
an n x p matrix where p is the number of parameters in o.
}
\description{
This is an S3 generic function that returns a
\code{n}x\code{PCMParamCount(o)} matrix of parameter vectors.
Implementations try to deduce parameter values based on the passed tree
and data, that should be suitable for starting points for likelihood
optimization or MCMC inference.
}
\examples{
library(PCMBase)
library(PCMFit)
modelBM.ab <- modelBM.ab.Guess <- PCM(
PCMDefaultModelTypes()["B"], k = 2, regimes = c("a", "b"))
modelBM.ab$X0[] <- c(5, 2)
# in regime 'a' the traits evolve according to two independent BM
# processes (starting from the global vecto X0).
modelBM.ab$Sigma_x[,, "a"] <- rbind(c(1.6, 0.01),
c(0, 2.4))
# in regime 'b' there is a correlation between the traits
modelBM.ab$Sigma_x[,, "b"] <- rbind(c(2.4, 0.08),
c(0.0, 0.8))
modelOU.ab.Guess <- PCM(
PCMDefaultModelTypes()["F"], k = 2, regimes = c("a", "b"))
modelMG.ab.Guess <- MixedGaussian(
k = 2, modelTypes = MGPMDefaultModelTypes(),
mapping = c(
a = unname(MGPMDefaultModelTypes()["B"]),
b = unname(MGPMDefaultModelTypes()["F"])))
param <- double(PCMParamCount(modelBM.ab))
# load the current model parameters into param
PCMParamLoadOrStore(modelBM.ab, param, offset=0, load=FALSE)
print(param)
# make results reproducible
set.seed(2, kind = "Mersenne-Twister", normal.kind = "Inversion")
# number of regimes
R <- 2
# number of extant tips
N <- 100
tree.a <- PCMTree(ape::rtree(n=N))
PCMTreeSetLabels(tree.a)
PCMTreeSetPartRegimes(
tree.a, part.regime = c(`101` = "a"), setPartition = TRUE)
lstDesc <- PCMTreeListDescendants(tree.a)
splitNode <- names(lstDesc)[which(sapply(lstDesc, length) > N/2 &
sapply(lstDesc, length) < 2*N/3)][1]
tree.ab <- PCMTreeInsertSingletons(
tree.a, nodes = as.integer(splitNode),
positions = PCMTreeGetBranchLength(tree.a, as.integer(splitNode))/2)
PCMTreeSetPartRegimes(
tree.ab,
part.regime = structure(
c("a", "b"), names = as.character(c(N+1, splitNode))),
setPartition = TRUE)
traits <- PCMSim(tree.ab, modelBM.ab, modelBM.ab$X0)
likFunBM <- PCMCreateLikelihood(traits, tree.ab, modelBM.ab.Guess)
likFunOU <- PCMCreateLikelihood(traits, tree.ab, modelOU.ab.Guess)
likFunMG <- PCMCreateLikelihood(traits, tree.ab, modelMG.ab.Guess)
vecGuessBM <-
GuessInitVecParams(modelBM.ab.Guess, X = traits, tree = tree.ab)
vecGuessOU <-
GuessInitVecParams(modelOU.ab.Guess, X = traits, tree = tree.ab)
vecGuessMG <-
GuessInitVecParams(modelMG.ab.Guess, X = traits, tree = tree.ab)
# likelihood at true parameters used to generate the data
print(param)
likFunBM(param)
# likelihood at guessed parameter values for the BM model
print(vecGuessBM)
likFunBM(vecGuessBM)
# likelihood at guessed parameter values for the OU model
print(vecGuessOU)
likFunOU(vecGuessOU)
# likelihood at guessed parameter values for the BM model
print(vecGuessMG)
likFunMG(vecGuessMG)
# likelihood at completely random parameter values
vecRand <- PCMParamRandomVecParams(modelBM.ab.Guess, k = 2, R = 2)
likFunBM(vecRand)
}