1 ## extract.popsize.R (2004-07-4)
3 ## Extract table with population size in dependence of time
4 ## from mcmc output generated by mcmc.popsize
6 ## Copyright 2004 Rainer Opgen-Rhein and Korbinian Strimmer
8 ## This file is part of the R-package `ape'.
9 ## See the file ../COPYING for licensing issues.
11 extract.popsize<-function(mcmc.out, credible.interval=0.95, time.points=200, thinning=1, burn.in=0)
14 # construct a matrix with the positions of the jumps
17 k<-array(dim=ceiling((length(mcmc.out$pos)-burn.in)/thinning))
19 k[i]<-length(mcmc.out$pos[[b]]);
27 pos.m<-matrix(nrow=length(k), ncol=o)
29 pos.m[i,]<-c(mcmc.out$pos[[b]], array(dim=o-length(mcmc.out$pos[[b]])));
34 # construct a matrix with the heights of the jumps
37 h.m<-matrix(nrow=length(k), ncol=o)
39 h.m[i,]<-c(mcmc.out$h[[b]], array(dim=o-length(mcmc.out$h[[b]])));
43 prep<-list("pos"=pos.m, "h"=h.m)
47 step <- (max(prep$pos, na.rm=TRUE)-min(prep$pos, na.rm=TRUE))/(time.points-1)
50 p<-min(prep$pos, na.rm=TRUE)
52 me<-matrix(nrow=nr, ncol=5)
55 prep.l$pos<-cbind(prep$pos,prep$pos[,length(prep$pos[1,])])
56 prep.l$h<-cbind(prep$h,prep$h[,length(prep$h[1,])])
58 while (p<=max(prep$pos, na.rm=TRUE))
60 #Vector with position of heights
62 l.prep[is.na(l.prep)]<-FALSE
63 pos.of.h<-l.prep%*% array(data=1, dim=dim(prep$pos)[2])
66 z<-array(data=(1:dim(prep$pos)[1]), dim=dim(prep$pos)[1])
67 index.left<-cbind(z,pos.of.h)
68 index.right<-cbind(z, pos.of.h+1)
70 mixed.heights<-((((p-prep$pos[index.left])/(prep$pos[index.right]-prep$pos[index.left]))*
71 (prep$h[index.right]-prep$h[index.left]))+prep$h[index.left])
73 me[i,2]<-mean(mixed.heights)
76 #me[i,2]<-huber(mixed.heights)$mu
78 me[i,3]<-median(mixed.heights)
79 me[i,4]<-quantile(mixed.heights, probs=(1-credible.interval)/2, na.rm=TRUE)
80 me[i,5]<-quantile(mixed.heights, probs=(1+credible.interval)/2, na.rm=TRUE)
86 #av.jumps<-round((length(prep$pos)-sum(is.na(prep$pos)))/length(prep$pos[,1])-2,2)
87 #print("average jumps")
89 #print((length(prep$pos)-sum(is.na(prep$pos)))/length(prep$pos[,1])-2)
91 colnames(me) <- c("time", "mean", "median", "lower CI", "upper CI")
92 class(me) <- "popsize"