R/extract.popsize.R

   1 ## extract.popsize.R (2004-07-4)
   2
   3 ##  Extract table with population size in dependence of time
   4 ##      from mcmc output generated by mcmc.popsize
   5
   6 ## Copyright 2004 Rainer Opgen-Rhein and Korbinian Strimmer
   7
   8 ## This file is part of the R-package `ape'.
   9 ## See the file ../COPYING for licensing issues.
  10
  11 extract.popsize<-function(mcmc.out, credible.interval=0.95, time.points=200, thinning=1, burn.in=0)
  12 {
  13
  14   # construct a matrix with the positions of the jumps
  15     b<-burn.in+1
  16     i<-1
  17     k<-array(dim=ceiling((length(mcmc.out$pos)-burn.in)/thinning))
  18     while(i<=length(k)) {
  19        k[i]<-length(mcmc.out$pos[[b]]);
  20        (i<-i+1);
  21        b<-b+thinning
  22     }
  23     o<-max(k)
  24
  25     b<-burn.in+1
  26     i<-1
  27     pos.m<-matrix(nrow=length(k), ncol=o)
  28     while(i<=length(k)) {
  29         pos.m[i,]<-c(mcmc.out$pos[[b]], array(dim=o-length(mcmc.out$pos[[b]])));
  30         i<-i+1;
  31         b<-b+thinning
  32     }
  33
  34   # construct a matrix with the heights of the jumps
  35     b<-burn.in+1
  36     i<-1
  37     h.m<-matrix(nrow=length(k), ncol=o)
  38     while(i<=length(k)) {
  39         h.m[i,]<-c(mcmc.out$h[[b]], array(dim=o-length(mcmc.out$h[[b]])));
  40         i<-i+1;
  41         b<-b+thinning
  42      }
  43   prep<-list("pos"=pos.m, "h"=h.m)
  44
  45 ####################
  46
  47   step <- (max(prep$pos, na.rm=TRUE)-min(prep$pos, na.rm=TRUE))/(time.points-1)
  48   nr <- time.points
  49
  50   p<-min(prep$pos, na.rm=TRUE)
  51   i<-1
  52   me<-matrix(nrow=nr, ncol=5)
  53
  54   prep.l<-prep
  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,])])
  57
  58   while (p<=max(prep$pos, na.rm=TRUE))
  59   {
  60     #Vector with position of heights
  61     l.prep<-prep$pos<=p
  62     l.prep[is.na(l.prep)]<-FALSE
  63     pos.of.h<-l.prep%*% array(data=1, dim=dim(prep$pos)[2])
  64
  65     #Vector with heights
  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)
  69
  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])
  72
  73     me[i,2]<-mean(mixed.heights)
  74
  75     #library(MASS)
  76     #me[i,2]<-huber(mixed.heights)$mu
  77
  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)
  81     me[i,1]<-p
  82     p<-p+step
  83     i<-i+1
  84   }
  85
  86   #av.jumps<-round((length(prep$pos)-sum(is.na(prep$pos)))/length(prep$pos[,1])-2,2)
  87   #print("average jumps")
  88
  89   #print((length(prep$pos)-sum(is.na(prep$pos)))/length(prep$pos[,1])-2)
  90
  91   colnames(me) <- c("time", "mean", "median", "lower CI", "upper CI")
  92   class(me) <- "popsize"
  93
  94   return(me)
  95 }