Nothing
test_this("Dynamic scoping", {
set.seed(1)
ntips = 200
k = 2 # No. of trait dimensions
tr = ape::rtree(ntips)
X = matrix(rnorm(k*ntips), k, ntips) # Trait matrix
x0 = rnorm(k) # Root value
repar = get_restricted_ou(H='diag', theta=NULL, Sig=NULL, lossmiss='halt')
my_par = function (par, ...) {
phiwV = repar$par(par[1:7], ...)
if (INFO__$node_id > 200) ## If not tip just return the original
return(phiwV)
Sig_e = diag(par[8:9]) ## Our measurement error matrix
phiwV[7:9] = phiwV[7:9] + Sig_e[lower.tri(Sig_e, diag=T)]
phiwV
}
my_jac = function (par, ...) {
new_jac = matrix(0.0, 9, 9)
new_jac[,1:7] = repar$jac(par[1:7], ...)
if (INFO__$node_id <= 200)
new_jac[7,8] = new_jac[9,9] = 1.0
new_jac
}
my_hess = function (par, ...)
lapply(repar$hess(par[1:7], ...), function (H) {
newH = array(0.0, dim=c(dim(H)[1], 9, 9))
newH[,1:7,1:7] = H[,,] ## Copy the original part
newH ## Other entries are just zero
})
mod = glinv(tr, x0, X,
pardims = 9,
parfns = my_par,
parjacs = my_jac,
parhess = my_hess)
expect_type(lik(mod)(c(1,1,0,0,0,0,0,0.5,0.5)), 'double')
expect_type(grad(mod)(c(1,1,0,0,0,0,0,0.5,0.5)), 'double')
expect_type(hess(mod)(c(1,1,0,0,0,0,0,0.5,0.5)), 'double')
})
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.