# Authors: Forrest W. Crawford and Marc A. Suchard ################################### library(TreeSim) source("yulerewards.r") source("treeutils.r") ################################### # simulate topologies and brownian motions # to estimate a known variance. sim_estimate_variance = function(ntopologies=100) { cat("this might take a while...\n") quartz(width=8,height=6) par(mfrow=c(4,1)) par(mar=c(4,4,0,0)) num_species = 3:12 sigma2 = 1 nbranchingtimes_simtype = c(1,5,10,100) nsim_simtype = c(1,5,10,100) for(simtype in 1:4) { plot(c(0,0), type="n", ylim=c(0,sigma2+4), xlim=c(min(num_species)-0.25,max(num_species)+0.25), main="", xlab="", ylab="", frame.plot=FALSE, axes=FALSE) if(simtype==4) axis(1, at=num_species) axis(2, at=0:(sigma2+2)) lines(c(min(num_species), max(num_species)), c(sigma2,sigma2), lwd=2, lty="dashed", col="gray") nbranchingtimes = nbranchingtimes_simtype[simtype] nsim = nsim_simtype[simtype] #nbranchingtimes = 1 #nsim = 1 t = 1 lam = 2 estimates = array(NA, dim=c(length(num_species), ntopologies)) for(i in 1:length(num_species)) { n = num_species[i] cat("n =", n, "\n") for(j in 1:ntopologies) { reftree = sim.bd.taxa.age(n, 1, lambda=lam,mu=0,age=t,mrca=TRUE)[[1]] a = get_rewards(reftree) md = rep(0,nbranchingtimes) for(k in 1:nbranchingtimes) { newtree = resample_yule_branch_lengths(reftree,lam) d = rep(0,nsim) for(l in 1:nsim) { x = rTraitCont(newtree, sigma=sqrt(sigma2)) #d[j] = var(x) d[l] = sum((x-mean(x))^2) / n } md[k] = mean(d) } f = function(v) {getmloglik_rewards(md,v,1,n,t,lam,a)} vmin = min(md)/get_maxreward(2,n,t,a) vmax = max(md)/get_minreward(2,n,t,a) #cat(" estimate in interval (",vmin,",",vmax,")\n") vs = seq(vmin,vmax,length.out=10)[c(-1,-10)] vstart = vs[which.min(sapply(vs, f))] #cat("vstart =", vstart, "\n") obj = nlm(f, vstart, stepmax=0.3) #cat(" code =", obj$code, "iterations =", obj$iterations, "\n") estimates[i,j] = ifelse(obj$iterations==0,NA,obj$estimate) cat("topology", j, "/", ntopologies, ": vhat =", obj$estimate, "iters =", obj$iterations, "code =", obj$code, "\n") } points(jitter(rep(n,ntopologies),amount=0.05), estimates[i,], pch=20) } } } sim_estimate_variance(10)