Commit 48820d87 authored by Nicolas Salamin's avatar Nicolas Salamin
Browse files

final version that works and is up-to-date

parent 1724e908
......@@ -59,9 +59,9 @@ condLike<-function(sigma2=1.0, condvec, kids) {
return(c(lnL, anc, S));
}
#set the condL for the tips based on the data and compute the distance between each node and the root
#set the condL for the tips based on the data and the internal nodes
condvecInit<-function(tree, data, ntips, tnodes) {
condvec<-matrix(nrow=tnodes, ncol=4); #three columns, the first for the logLn, the second for the state, the third for extra BM variance and the fourth is the brlen subtending the node
condvec<-matrix(nrow=tnodes, ncol=4); #four columns: the first for the logLn, the second for the state, the third for extra BM variance and the fourth is the brlen subtending the node
for(i in 1:tnodes) {
row<-which(tree$edge[,2] == i);
......@@ -116,30 +116,31 @@ lnLike<-function(tree, data, sigma2, getLnl=TRUE) {
}
}
library(mnormt)
fitBM <- function(tree, x, anc, sigma2, pruning=TRUE, optim=TRUE, getLnl=TRUE) {
n <- length(x);
if(optim) {
if(pruning) {
lik <- function(par, tree, x) lnLike(tree, x, par, getLnl=TRUE); #if you want to optimize, the function MUST return the likelihood, whatever the wish of the user...
fit <- optim(c(1), lik, tree = tree, x = x, method = "L-BFGS-B", lower = c(1e-12), upper = c(Inf), control = list(fnscale = -1));
return(list(sig2 = fit$par, logL = fit$value, convergence = fit$convergence));
}
else {
lik <- function(par, tree, x, n) sum(dmnorm(x, mean = rep(par[1], n), varcov = par[2] * vcv(tree), log = TRUE));
fit <- optim(c(0, 1), lik, tree = tree, x = x, n = n, method = "L-BFGS-B", lower = c(-Inf, 1e-12), upper = c(Inf, Inf), control = list(fnscale = -1));
return(list(x0 = fit$par[1], sig2 = fit$par[2], logL = fit$value, convergence = fit$convergence));
}
n <- length(x);
if(optim) {
if(pruning) {
lik <- function(par, tree, x) lnLike(tree, x, par, getLnl=TRUE); #if you want to optimize, the function MUST return the likelihood, whatever the wish of the user...
fit <- optim(c(1), lik, tree = tree, x = x, method = "L-BFGS-B", lower = c(1e-12), upper = c(Inf), control = list(fnscale = -1));
return(list(sig2 = fit$par, logL = fit$value, convergence = fit$convergence));
}
else {
if(pruning) {
return(lnLike(tree, x, sigma2, getLnl=getLnl));
}
else {
return(sum(dmnorm(x, mean = rep(anc, n), varcov = sigma2 * vcv(tree), log = TRUE)));
}
require(mnormt);
lik <- function(par, tree, x, n) sum(dmnorm(x, mean = rep(par[1], n), varcov = par[2] * vcv(tree), log = TRUE));
fit <- optim(c(0, 1), lik, tree = tree, x = x, n = n, method = "L-BFGS-B", lower = c(-Inf, 1e-12), upper = c(Inf, Inf), control = list(fnscale = -1));
return(list(x0 = fit$par[1], sig2 = fit$par[2], logL = fit$value, convergence = fit$convergence));
}
}
else {
if(pruning) {
return(lnLike(tree, x, sigma2, getLnl=getLnl));
}
else {
require(mnormt);
return(sum(dmnorm(x, mean = rep(anc, n), varcov = sigma2 * vcv(tree), log = TRUE)));
}
}
}
###########################
......@@ -189,4 +190,17 @@ plot(res, ylim=c(0,2*trueS))
abline(h=trueS, col="red")
plot(anc, ylim=c(0, 3*trueA))
abline(h=trueA, col="red")
\ No newline at end of file
abline(h=trueA, col="red")
########################
### Running the code ###
########################
tree<-pbtree(n=20, scale=1);
data<-data.frame(char=fastBM(tree, sig2=2, a=0), row.names=tree$tip.label);
#step 1: estimate the sigma
s2<-fitBM(tree, data[,1], pruning=TRUE, optim=TRUE, getLnl=TRUE)$sig2
#step 2: getting the values at each node:
fitBM(tree, data[,1], sigma2=s2, pruning=TRUE, optim=FALSE, getLnl=FALSE)
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment