R/classes.R

Defines functions plot.depmeasure summary.depmeasure print.depmeasure is.depmeasure depmeasure summary.bayesparams print.bayesparams is.bayesparams bayesparams residual.distributions residual.densities plot.bayesfit summary.bayesfit print.bayesfit is.bayesfit bayesfit plot.stepfit summary.stepfit print.stepfit is.stepfit stepfit

Documented in bayesfit bayesparams depmeasure is.bayesfit is.bayesparams is.depmeasure is.stepfit plot.bayesfit plot.depmeasure plot.stepfit print.bayesfit print.bayesparams print.depmeasure print.stepfit stepfit summary.bayesfit summary.bayesparams summary.depmeasure summary.stepfit

## Copyright (C) 2017 Thomas Lugrin
## Definition of methods and related functions
## (Scope:) single framework for 1- and 2-stage functions
## (Note:) summary is tricky especially when multi-dimensional!
## List of functions: - stepfit
##                    - is.stepfit
##                    - print/summary.stepfit
##                    - plot.stepfit
##                    - bayesfit
##                    - is.bayesfit
##                    - print/summary.bayesfit
##                    - plot.bayesfit
##                    - bayesparams
##                    - is.bayesparams
##                    - print/summary.bayesparams
##                    - depmeasure
##                    - is.depmeasure
##                    - print/summary.depmeasure
##                    - plot.depmeasure
##                    - is.int
##################################################

##################################################
## CLASS "STEPFIT"
## constructor
stepfit <- function(){
  x <- list(a=numeric(0), b=numeric(0), res=matrix(0,nrow=1,ncol=1), pars.se=matrix(0,nrow=1,ncol=2), nlag=0)
  class(x) <- "stepfit"
  return(x)
}

## validator
is.stepfit <- function(x){
  if(length(names(x))){
    conds <- names(stepfit()) %in% names(x)# minimal requirement
    conds <- sum(conds) == length(names(stepfit()))
    if(conds){
      conds <- is.numeric(x$a) && is.numeric(x$b) && is.matrix(x$res) && is.matrix(x$pars.se)
      if(!conds) return(FALSE)
      conds <- dim(x$pars.se)[1]==2 && length(x$a)==length(x$b) && dim(x$pars.se)[2]==length(x$a) && length(x$nlag) == 1
      return(inherits(x, "stepfit") && conds)
    }
  }
  return(FALSE)
}

## viewer
print.stepfit <- function(x, ...){
  summary.stepfit(x, ...)
}

summary.stepfit <- function(object, ...){
  cat(" Parameters:\n")
  for(ind in 1:object$nlag){
    cat(paste("alpha(",ind,")  ",signif(object$a,3)," (",signif(object$pars.se[1,ind],3),")\n", sep=""), sep="")
    cat(paste("beta(",ind,")   ",signif(object$b,3)," (",signif(object$pars.se[2,ind],3),")\n", sep=""), sep="")
  }
  cat(" Residuals:\n")
  print(summary(object$res))
}

## plotter
plot.stepfit <- function(x, ...){
  old.par <- par(no.readonly=TRUE)
  on.exit(par(old.par))
  par(ask=TRUE)
  def <- list(xlab="z", ylab=bquote(h(z)), main="")
  if(length(list(...)) > 0){
    pm <- !pmatch(names(def), names(list(...)), 0)# all non-matches
    def <- def[pm]
  }
  for(j in 1:x$nlag){
    if("ylab" %in% names(def)) def$ylab <- bquote(bquote(h(z[.(j)])))
    if("xlab" %in% names(def)) def$xlab <- bquote(bquote(z[.(j)]))
    den <- density(x$res[,j])
    def$x <- x$res[,j]
    do.call(hist, c(def, prob=TRUE, list(...)))
    lines(den, lty="dashed")
  }
}


##################################################
## CLASS "BAYESFIT"
## constructor
bayesfit <- function(){
  z.mat <- matrix(0, nrow=1, ncol=1)
  z.ar  <- array(0, dim=c(1,1,1))
  x <- list(a=z.mat, b=z.mat, sd=z.ar, mean=z.ar, w=z.mat, prec=0,
              ci=z.mat, noo=z.mat, noc=0, prop.sd=z.ar,
              len=0, nlag=0)
  class(x) <- "bayesfit"
  return(x)
}

## validator
is.bayesfit <- function(x){
  if(length(names(x))){
    conds <- names(bayesfit()) %in% names(x)
    conds <- sum(conds) == length(names(bayesfit()))
    if(all(conds)){
      conds <- is.matrix(x$a) && is.matrix(x$b) &&
        is.array(x$sd) && is.array(x$mean) &&
        is.matrix(x$w) && is.numeric(x$prec) && is.array(x$prop.sd) &&
        is.matrix(x$ci) && is.matrix(x$noo) && is.numeric(x$noc) && is.numeric(x$len) && is.numeric(x$nlag) &&
        is.int(x$len) && length(x$len)==1 && is.int(x$nlag) && length(x$nlag)==1 &&
        is.int(x$noc) && is.int(x$noo) && is.int(x$ci)
      return(inherits(x, "bayesfit") && conds)
    }
  }
  return(FALSE)
}

## viewer
print.bayesfit <- function(x, ...){
  summary.bayesfit(x, ...)
}

summary.bayesfit <- function(object, ...){
  cat(" Posterior medians:\n")
  for(j in 1:object$nlag){
    cat("alpha(",j,")  ",median(object$a[,j]),"\n", sep="")
    cat("beta(",j,")   ",median(object$b[,j]),"\n", sep="")
    cat("mean(",j,",",1,") ",median(object$mean[,1,j]),"\n", sep="")
    cat("sd(",j,",",1,")   ",median(object$sd[,1,j]),"\n", sep="")
  }
}

plot.bayesfit <- function(x, which=1:3, ...){
  old.par <- par(no.readonly=TRUE)
  on.exit(par(old.par))
  stopifnot(is.bayesfit(x))
  stopifnot(is.int(which))
  stopifnot(min(which) >= 1 && max(which) <= 3)
  which <- sort(unique(which))
  ugm <- all(x$a==0) && all(x$b==0)
  if(ugm) which <- intersect(which,1:2)
  if(x$nlag > 1) par(ask=TRUE)
  if(length(which) == 2) par(mfrow=c(1,2))
  else if(length(which) > 2) par(mfrow=c(2,2))
  def <- list(xlab="x", ylab="y", main="")
  if(length(list(...)) > 0){
    pm <- !pmatch(names(def), names(list(...)), 0)# all non-matches
    def <- def[pm]
  }
  for(j in 1:x$nlag){
    if(1 %in% which){
      if("xlab" %in% names(def)) def$xlab <- bquote(bquote(z[.(j)]))
      if("ylab" %in% names(def)) def$ylab <- bquote(bquote(h(z[.(j)])))
      if("main" %in% names(def)) def$main <- "Sample of residual densities"
        grid <- seq(-100,100,2)
        de <- residual.densities(x,j,grid)
        grid <- seq(min(de$x),max(de$x),length.out=151)
        de <- residual.densities(x,j,grid)
        xy <- list(x=de$x, y=de$y[,1])
        if(!("ylim" %in% names(list(...)))) xy$ylim <- c(0,max(de$y))
        do.call(plot, c(xy, type="l",def,list(...)))
        for(it in 1:dim(de$y)[2])
          lines(de$x, de$y[,it])
    }
    if(2 %in% which){
      if("xlab" %in% names(def)) def$xlab <- bquote(bquote(z[.(j)]))
      if("ylab" %in% names(def)) def$ylab <- bquote(bquote(H(z[.(j)])))
      if("main" %in% names(def)) def$main <- "Residual distribution"
      grid <- seq(-100,100,2)
      ds <- residual.distributions(x,j,grid)
      grid <- seq(min(ds$x),max(ds$x),length.out=151)
      ds <- residual.distributions(x,j,grid)
      xy <- list(x=ds$x, y=ds$y[,1])
      do.call(plot, c(xy, type="l",def,list(...)))
      for(q in 2:3)
        lines(ds$x, ds$y[,q], lty="dashed")
    }
    if(3 %in% which){
      if("xlab" %in% names(def)) def$xlab <- bquote(bquote(alpha[.(j)]))
      if("ylab" %in% names(def)) def$ylab <- bquote(bquote(beta[.(j)]))
      if("main" %in% names(def)) def$main <- "Joint posterior of H-T par."
      lims <- c(max(min(x$a[,j]),-1),min(max(x$a[,j]),1),max(min(x$b[,j]),0),min(max(x$b[,j]),1))
      kd <- kde2d(x$a[,j], x$b[,j], lims=lims)
      do.call(contour, c(kd,nlevels=5,def,list(...)))
    }
  }
}

residual.densities <- function(x, lag, grid){
  nsamp<- min(25,x$len)# nbr of samples
  noc  <- dim(x$mean)[2]
  nobs <- sum(x$noo[1,])
  ind <- sample.int(x$len, nsamp, replace=FALSE)
  dens <- vapply(ind, function(it,j){
    de <- vapply(1:noc, function(c,j,it){
        x$noo[it,c]*dnorm(grid, x$mean[it,c,j], x$sd[it,c,j])
      }, grid, j=j, it=it)
      rowSums(de)# sample density
    }, grid, j=lag)
  dens <- dens/nobs
  maxs <- vapply(1:nsamp, function(j) max(dens[,j]), 0)
  wmin <- which.min(maxs)# expected to be the flattest
  minimax <- min(maxs)
  fact <- 1e3# trim tails which are at least [fact] times smaller than minimax
  lo <- 1; up <- length(grid)
  i <- 1
  while(i < up && dens[i,wmin] < minimax/fact)
    i <- i+1
  lo <- max(lo,i-3); i <- up
  while(i > lo && dens[i,wmin] < minimax/fact)
    i <- i-1
  up <- min(up,i+3)
  return(list(x=grid[lo:up], y=dens[lo:up,]))
}

residual.distributions <- function(x, lag, grid){
  nsamp<- min(500,x$len)# nbr of samples
  noc  <- dim(x$mean)[2]
  nobs <- sum(x$noo[1,])
  ind <- sample.int(x$len, nsamp, replace=FALSE)
  dstr <- vapply(ind, function(it, j){
    ds <- vapply(1:noc, function(c,j,it){
      x$noo[it,c]*pnorm(grid, x$mean[it,c,j], x$sd[it,c,j])
    }, grid, j=j, it=it)
    rowSums(ds)# sample distribution
  }, grid, j=lag)
  dstr <- dstr/nobs
  dstr <- vapply(seq_along(grid), function(gr){
    quantile(dstr[gr,], probs=c(.5,.05,.95))
  }, numeric(3))
  dstr <- t(dstr)
  # trim tails of distribution
  lo <- 1; up <- length(grid)
  if(grid[2]-grid[1] > 1){
    eps.prob <- 1e-3
    i <- 1
    while(i < up && dstr[i,1] < eps.prob)
      i <- i+1
    lo <- max(lo,i-1); i <- up
    while(i > lo && dstr[i,1] > 1-eps.prob)
      i <- i-1
    up <- min(up,i+1)
  }
  return(list(x=grid[lo:up], y=dstr[lo:up,]))
}


##################################################
## CLASS "BAYESPARAMS"
## constructor
bayesparams <- function(prop.a=0.02, prop.b=0.02,
                        prior.mu=c(0,10), prior.nu=c(2,1/2), prior.eta=c(2,2),
                        trunc=100, comp.saved=15,
                        maxit=30000, burn=5000, thin=1,
                        adapt=5000, batch.size=125, mode=1){
  x <- list(prop.a=prop.a, prop.b=prop.b,
              prior.mu=prior.mu, prior.nu=prior.nu, prior.eta=prior.eta,
              trunc=trunc, comp.saved=comp.saved,
              maxit=maxit, burn=burn, thin=thin,
              adapt=adapt, batch.size=batch.size, mode=mode)
  class(x) <- "bayesparams"
  if(is.bayesparams(x)) return(x)
  else stop("Wrong type of argument. See help for details.")
}

## validator
is.bayesparams <- function(x){
  if(length(names(x))){
    names.bp <- c("prop.a","prop.b","prior.mu","prior.nu","prior.eta",
                  "trunc","comp.saved","maxit","burn","thin",
                  "adapt","batch.size","mode")
    conds <- names.bp %in% names(x)
    if(all(conds)){
      conds <- is.numeric(x$prop.a) && is.numeric(x$prop.b) &&
        length(x$prop.a)==1 && length(x$prop.b)==1 && x$prop.a>0 && x$prop.b>0 &&
        is.numeric(x$prior.mu) && is.numeric(x$prior.nu) && is.numeric(x$prior.eta) &&
        length(x$prior.mu)==2 && length(x$prior.nu)==2 && length(x$prior.eta)==2 &&
        x$prior.mu[2] > 0 && all(x$prior.nu>0) && all(x$prior.eta>0) &&
        length(x$trunc)==1 && is.int(x$trunc) && length(x$comp.saved)==1 && is.int(x$comp.saved) &&
        length(x$maxit)==1 && is.int(x$maxit) && length(x$burn)==1 && is.int(x$burn) &&
        length(x$thin)==1 && is.int(x$burn) && length(x$adapt)==1 && is.int(x$adapt) &&
        length(x$batch.size)==1 && is.int(x$batch.size) && length(mode)==1 && x$mode %in% c(0,1,2)
      return(inherits(x, "bayesparams") && conds)
    }     
  }
  return(FALSE)
}

## viewer
print.bayesparams <- function(x, ...){
  summary.bayesparams(x, ...)
}

summary.bayesparams <- function(object, ...){
  if(!is.bayesparams(object)){
    summary.default(object)
  }else{
    cat("proposal for alpha            ",object$prop.a,"\n", sep="")
    cat("proposal for beta             ",object$prop.b,"\n", sep="")
    cat("prior par. for means          (",object$prior.mu[1],", ",object$prior.mu[2],")\n", sep="")
    cat("prior par. for sd             (",object$prior.nu[1],", ",object$prior.nu[2],")\n", sep="")
    cat("prior par. for precision par. (",object$prior.eta[1],", ",object$prior.eta[2],")\n", sep="")
    cat("DP truncation                 ",object$trunc,"\n", sep="")
    cat("nbr of components saved       ",object$comp.saved,"\n", sep="")
    cat("max. nbr of MCMC iterations   ",object$maxit,"\n", sep="")
    cat("length of burn-in             ",object$burn,"\n", sep="")
    cat("thinning par.                 ",object$thin,"\n", sep="")
    cat("length of adaption            ",object$adapt,"\n", sep="")
    cat("size of adaption batches      ",object$batch.size,"\n", sep="")
    cat("debug mode                    ",object$mode,"\n", sep="")
  }
  invisible()
}


##################################################
## CLASS "DEPMEASURE"
## constructor
depmeasure <- function(type=c("theta","chi","steptheta","runs")){
  z.mat <- matrix(0, nrow=1, ncol=1)
  x <- list(probs=0, levels=0, nlag=1)
  if(type[1]=="theta" || type[1]=="chi"){
    x$fit   <- bayesfit()
    x$distr <- z.mat
    if(type[1]=="theta")
      x$theta <- z.mat
    else if(type[1]=="chi")
      x$chi <- z.mat
  }else if(type[1]=="steptheta"){
    x$fit <- stepfit()
    x$theta <- z.mat
  }else if(type[1]=="runs"){
    x$theta <- z.mat
    x$nbr.exc <- 0
  }else{
    stop("type not recognised. Must be either 'theta' or 'chi'.")
  }
  class(x) <- "depmeasure"
  return(x)
}

## validator
is.depmeasure <- function(x){
  if(length(names(x))){
    conds <- names(depmeasure("theta")) %in% names(x)
    if(all(conds)){
      conds <- is.bayesfit(x$fit) && is.matrix(x$theta)
      return(inherits(x, "depmeasure") && conds)
    }
    conds <- names(depmeasure("chi")) %in% names(x)
    if(all(conds)){
      conds <- is.bayesfit(x$fit) && is.matrix(x$chi)
      return(inherits(x, "depmeasure") && conds) 
    }
    conds <- names(depmeasure("steptheta")) %in% names(x)
    if(all(conds)){
      conds <- is.stepfit(x$fit) && is.matrix(x$theta)
      return(inherits(x, "depmeasure") && conds) 
    }
    conds <- names(depmeasure("runs")) %in% names(x)
    if(all(conds)){
      conds <- is.matrix(x$theta)
      return(inherits(x, "depmeasure") && conds)
    }
  }
  return(FALSE)
}

## viewer
print.depmeasure <- function(x, ...){
  summary.depmeasure(x, ...)
}

summary.depmeasure <- function(object, ...){
  cat("Bayes fit\n")
  summary(object$fit)
  cat(paste(names(object)[2],"posterior estimates\n"))
  print(summary(object[[2]]))
  invisible()
}

## plotter
plot.depmeasure <- function(x, ...){
  stopifnot(is.depmeasure(x))
  if("nbr.exc" %in% names(x)){
    ord <- order(colnames(x$theta))
    def <- list(xlab="x", ylab=bquote(bquote(theta(x,.(x$nlag)))), ylim=range(x$theta))
    if(length(list(...)) > 0){
      pm <- !pmatch(names(def), names(list(...)), 0)# all non-matches
      def <- def[pm]
    }
    def$x <- x$levels; def$y <- x$theta[,"estimate"]
    do.call(plot, c(type="l", def, list(...)))
    ord <- ord[-length(ord)]
    # add credible/confidence lines
    if(length(ord) > 1){
      for(i in 0:(floor(length(ord)/2)-1)){
        lines(x$levels, x$theta[,ord[i*2+1]], lty=i+2)
        lines(x$levels, x$theta[,ord[i*2+2]], lty=i+2)
      }
    }
  }else if("chi" %in% names(x)){
    ord <- order(colnames(x$chi))
    def <- list(xlab="x", ylab=bquote(bquote(chi[.(x$nlag)](x))), ylim=range(x$chi))
    if(length(list(...)) > 0){
      pm <- !pmatch(names(def), names(list(...)), 0)
      def <- def[pm]
    }
    def$x <- x$levels; def$y <- x$chi[,"median"]
    do.call(plot, c(def,type="l",list(...)))
    ord <- ord[-c(length(ord)-1,length(ord))]
    # add credible/confidence lines
    if(length(ord) > 1){
      for(i in 0:(floor(length(ord)/2)-1)){
        lines(x$levels, x$chi[,ord[i*2+1]], lty=i+2)
        lines(x$levels, x$chi[,ord[i*2+2]], lty=i+2)
      }
    }
  }else{
    ord <- order(colnames(x$theta))
    def <- list(xlab="x", ylab=bquote(bquote(theta(x,.(x$nlag)))), ylim=range(x$theta))
    if(length(list(...)) > 0){
      pm <- !pmatch(names(def), names(list(...)), 0)
      def <- def[pm]
    }
    def$x <- x$levels; ifelse(is.stepfit(x$fit), def$y <- x$theta[,"estimate"], def$y <- x$theta[,"median"])
    do.call(plot, c(def,type="l",list(...)))
    ord <- ord[-length(ord)]
    if(is.bayesfit(x$fit)) ord <- ord[-length(ord)]
    # add credible/confidence lines
    if(length(ord) > 1){
      for(i in 0:(floor(length(ord)/2)-1)){
        lines(x$levels, x$theta[,ord[i*2+1]], lty=i+2)
        lines(x$levels, x$theta[,ord[i*2+2]], lty=i+2)
      }
    }
  }
}


##################################################
## ADDITIONAL VALIDATORS
is.int <- function(x){
  if(all(x%/%1==x)) return(TRUE)
  else return(FALSE)
}

Try the tsxtreme package in your browser

Any scripts or data that you put into this service are public.

tsxtreme documentation built on April 24, 2021, 1:07 a.m.