Commit 79a4d7e6 authored by Nicolas Salamin's avatar Nicolas Salamin
Browse files

added a new function to calculate the probability with a known ancestor at a node

parent e2047e61
......@@ -40,11 +40,19 @@ condLike.old<-function(sigma=1.0, blen, condvec, is.list=TRUE) {
}
}
fN(0.732627, 0.314222, 1, 1, 1)
fN2(0.732627, 0.314222, 1, 1, 1, 0.523424)
sum(dmnorm(x, mean = rep(1, n), varcov = sigma * sigma * vcv(tree), log = TRUE))
fN<-function(xa, xb, vpa, vpb, sigma) {
#the normal density. Sigma is the standard deviation.
#the normal density. Sigma is the standard deviation. Same as log(dnorm(xa-xb, 0, sqrt(sigma * sigma * (vpa + vpb))))
return((-(xa - xb)^2 / (2 * sigma * sigma * (vpa + vpb)) - log(sqrt(2 * pi * sigma * sigma * (vpa + vpb)))));
#return(log(dnorm(xa-xb, sqrt((vpa + vpb) * sigma * sigma))));
}
fN2<-function(xa, xb, vpa, vpb, sigma, anc) {
#same as fN but with a known ancestor instead of xa-xb, see Joe's book eqn 23.10 (some mistakes in the denominator though)
return(-(((xa - anc)^2 / vpa) + ((xb - anc)^2 / vpb)) / (2 * sigma * sigma) - log(sqrt(2 * pi * sigma * sigma * (vpa + vpb))));
}
condLike<-function(sigma=1.0, condvec, kids, try=FALSE) {
......@@ -129,7 +137,7 @@ lnLike<-function(tree, data, anc, sigma) {
kids<-which(tree$edge[,1]==(ntips+1));
condvec[tnodes, 1:3]<-condLike(sigma=sigma, condvec=condvec, kids)
print(condvec[tnodes,])
print(condvec)
return(condvec[tnodes, 1] + fN(condvec[tnodes, 2], anc, condvec[tnodes, 3], 0, sigma));
}
......@@ -142,7 +150,7 @@ fitBM <- function(tree, x, anc, sigma, pruning=TRUE, optim=TRUE) {
lik <- function(par, tree, x, n) lnLike(tree, x, par[1], par[2]);
}
else {
lik <- function(par, tree, x, n) sum(dmnorm(x, mean = rep(par[1], n), varcov = par[2] * vcv(tree), log = TRUE));
lik <- function(par, tree, x, n) sum(dmnorm(x, mean = rep(par[1], n), varcov = par[2] * 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));
......@@ -163,40 +171,27 @@ fitBM <- function(tree, x, anc, sigma, pruning=TRUE, optim=TRUE) {
require(phytools);
tree<-pbtree(n=5, scale=1);
data<-data.frame(char=fastBM(tree, sig2=5, a=0), row.names=tree$tip.label);
tree<-pbtree(n=2, scale=1);
data<-data.frame(char=fastBM(tree, sig2=1, a=0), row.names=tree$tip.label);
fitBM(tree, data[,1], 0, 2.5, pruning=FALSE, optim=TRUE)
fitBM(tree, data[,1], 0, 2.5, pruning=TRUE, optim=TRUE)
fitBM(tree, data[,1], 0.646066, 0.816942, pruning=TRUE, optim=FALSE)
tree$edge
xxx<-condvecInit(tree, data[,1], 5, 9)
condLike(sigma=1, xxx, c(4,5), try=TRUE)
xxx[c(4,5),]
dnorm(-0.111984+1.561902, 0, 0.125411*2)
dnorm(-0.111984, -0.8369430, 0.125411)*2
lnLike(tree, data, 0, 2)
system.time(lnLike(tree, data, 0, 2.5))
s <- seq(0.9, 10, by=0.1)
s <- seq(0.5, 2.5, by=0.1)
res1 <- numeric(length(s))
res2 <- numeric(length(s))
for(i in 1:length(s)) {
res1[i] <- fitBM(tree, data[,1], -0.5, s[i], pruning=TRUE, optim=FALSE)
res2[i] <- fitBM(tree, data[,1], -0.5, s[i], pruning=FALSE, optim=FALSE)
res1[i] <- fitBM(tree, data[,1], 0.942176, s[i], pruning=TRUE, optim=FALSE)
res2[i] <- fitBM(tree, data[,1], 0.942176, s[i], pruning=FALSE, optim=FALSE)
}
plot(s, res2, t="l", col="red")
lines(s, res1)
which(res2 == max(res2))
optim(c(1), lnLike, tree=tree, data=data, method="L-BFGS-B", lower=c(1e-12), upper=c(Inf), control=list(fnscale=-1))
lnLike(tree, data, 0.7378779)
x<-seq(0.1, 5, 0.01)
plot(x, sapply(x, lnLike, tree=tree, data=data), type="l")
condLike(sigma=1, 0.01, c(0.0, 10, 0), is.list=FALSE)
condvecInit(tree, data[,1], 3, 5)
\ 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