## extract.popsize.R (2004-07-4)
## Extract table with population size in dependence of time
## from mcmc output generated by mcmc.popsize
## Copyright 2004 Rainer Opgen-Rhein and Korbinian Strimmer
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
extract.popsize<-function(mcmc.out, credible.interval=0.95, time.points=200, thinning=1, burn.in=0)
{
# construct a matrix with the positions of the jumps
b<-burn.in+1
i<-1
k<-array(dim=ceiling((length(mcmc.out$pos)-burn.in)/thinning))
while(i<=length(k)) {
k[i]<-length(mcmc.out$pos[[b]]);
(i<-i+1);
b<-b+thinning
}
o<-max(k)
b<-burn.in+1
i<-1
pos.m<-matrix(nrow=length(k), ncol=o)
while(i<=length(k)) {
pos.m[i,]<-c(mcmc.out$pos[[b]], array(dim=o-length(mcmc.out$pos[[b]])));
i<-i+1;
b<-b+thinning
}
# construct a matrix with the heights of the jumps
b<-burn.in+1
i<-1
h.m<-matrix(nrow=length(k), ncol=o)
while(i<=length(k)) {
h.m[i,]<-c(mcmc.out$h[[b]], array(dim=o-length(mcmc.out$h[[b]])));
i<-i+1;
b<-b+thinning
}
prep<-list("pos"=pos.m, "h"=h.m)
####################
step <- (max(prep$pos, na.rm=TRUE)-min(prep$pos, na.rm=TRUE))/(time.points-1)
nr <- time.points
p<-min(prep$pos, na.rm=TRUE)
i<-1
me<-matrix(nrow=nr, ncol=5)
prep.l<-prep
prep.l$pos<-cbind(prep$pos,prep$pos[,length(prep$pos[1,])])
prep.l$h<-cbind(prep$h,prep$h[,length(prep$h[1,])])
while (p<=max(prep$pos, na.rm=TRUE))
{
#Vector with position of heights
l.prep<-prep$pos<=p
l.prep[is.na(l.prep)]<-FALSE
pos.of.h<-l.prep%*% array(data=1, dim=dim(prep$pos)[2])
#Vector with heights
z<-array(data=(1:dim(prep$pos)[1]), dim=dim(prep$pos)[1])
index.left<-cbind(z,pos.of.h)
index.right<-cbind(z, pos.of.h+1)
mixed.heights<-((((p-prep$pos[index.left])/(prep$pos[index.right]-prep$pos[index.left]))*
(prep$h[index.right]-prep$h[index.left]))+prep$h[index.left])
me[i,2]<-mean(mixed.heights)
#library(MASS)
#me[i,2]<-huber(mixed.heights)$mu
me[i,3]<-median(mixed.heights)
me[i,4]<-quantile(mixed.heights, probs=(1-credible.interval)/2, na.rm=TRUE)
me[i,5]<-quantile(mixed.heights, probs=(1+credible.interval)/2, na.rm=TRUE)
me[i,1]<-p
p<-p+step
i<-i+1
}
#av.jumps<-round((length(prep$pos)-sum(is.na(prep$pos)))/length(prep$pos[,1])-2,2)
#print("average jumps")
#print((length(prep$pos)-sum(is.na(prep$pos)))/length(prep$pos[,1])-2)
colnames(me) <- c("time", "mean", "median", "lower CI", "upper CI")
class(me) <- "popsize"
return(me)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.