# R/helpingFunctions.R In SimoneHermann/BaPreStoPro: Bayesian Prediction of Stochastic Processes

#'
#'
#' Transformation of event times to NHPP
#'
#' @description Transformation of vector of event times to the corresponding counting process variables.
#' @param times vector of event times
#' @param t times of counting process
#'
#' @export
TimestoN <- function(times, t){
lt <- length(t)
dN <- numeric(lt)
ind <- unlist(sapply(times[times <= max(t)], function(a){which(abs(t-a) == min(abs(t-a)))}))
for(a in 1:length(ind)) dN[ind[a]] <- dN[ind[a]] + 1
cumsum(dN)
}

#' Transformation of NHPP variables to event times
#'
#' @description Vector of Poisson process differences are translated to a vector of event times.
#' @param dN vector of differences of counting process
#' @param t times of counting process
#' @export
dNtoTimes <- function(dN, t){
if(any(dN > 1)){
m <- sum(dN > 0)
res <- NULL
he1 <- dN[dN > 0]
he2 <- t[dN > 0]
for(mi in 1:m){
res <- c(res, rep(he2[mi], he1[mi]))
}
}else{
res <- t[dN > 0]
}
res
}

#' Sampling from lognormal proposal density
#'
#' @description Drawing one sample from the lognormal distribution with mean \code{parOld} and standard deviation \code{propSd}. Used in Metropolis Hastings algorithms.
#' @param parOld the parameter from the last iteration step
#' @param propSd proposal standard deviation
#'
#' @examples
#' plot(replicate(100, proposal(1, 0.1)), type = "l")
#' @export
proposal <- function(parOld, propSd){
if(any(parOld < 1e-150)) parOld[parOld < 1e-150] <- 1e-150  # 1e-320 equal to zero ...
mu <- log(parOld) - log( propSd^2/(parOld^2) + 1)/2
sigma2 <- log( propSd^2/(parOld^2)+1)
rlnorm(length(parOld), mu, sqrt(sigma2))
}

#' Proposal ratio of lognormal proposal density
#'
#' @param parOld the parameter from the last iteration step
#' @param parNew drawn candidate
#' @param propSd proposal standard deviation
#'
#' @examples
#' cand <- proposal(1, 0.01)
#' proposalRatio(1, cand, 0.01)
#' @export
proposalRatio <- function(parOld, parNew, propSd){
muOld <- log(parOld) - log( propSd^2/exp(2*log(parOld)) + 1)/2
sigma2Old <- log( propSd^2/exp(2*log(parOld))+1)
muNew <- log(parNew) - log( propSd^2/exp(2*log(parNew)) + 1)/2
sigma2New <- log( propSd^2/exp(2*log(parNew))+1)

prod(dlnorm(parOld, muNew, sqrt(sigma2New))/dlnorm(parNew, muOld, sqrt(sigma2Old)))
}

#' Inversion Method
#'
#' @description Algorithm to sample from cumulative distribution function, if no inverse function is analytically available.
#' @param Fun cumulative distribution function
#' @param len number of samples
#' @param candArea candidate area
#' @param grid fineness degree
#' @param method vectorial ("vector") or not ("free")
#' @examples
#' test <- InvMethod(function(x) pnorm(x, 5, 1), 1000, candArea = c(0, 10), method = "free")
#' plot(density(test))
#' curve(dnorm(x, 5, 1), col = 2, add = TRUE)
#' @references
#' Devroye, L. (1986). Non-Uniform Random Variate Generation. New York: Springer.
#' @export
InvMethod <- function(Fun, len, candArea, grid = 1e-05, method = c("vector", "free")){
method <- match.arg(method)
if(missing(candArea)){
candArea <- findCandidateArea(Fun)
if(method == "vector") candArea <- seq(candArea[1], candArea[2], by = grid)
}else{
if(method == "vector" & length(candArea) == 2) candArea <- seq(candArea[1], candArea[2], by = grid)
if(method == "free" & length(candArea) > 2) candArea <- c(min(candArea), max(candArea))
}

if(method == "vector"){
diFu <- sapply(candArea, Fun)
U <- runif(len, 0, max(diFu))
res <- sapply(U, function(u) candArea[which(diFu >= u)[1]])
}
if(method == "free"){
res <- numeric(len)
U <- runif(len)
for(i in 1:len){
lower <- candArea[1]
upper <- candArea[2]

diff <- upper - lower
while(diff >= grid){
if(Fun(lower+diff/2) < U[i]){
lower <- lower+diff/2
}else{
upper <- lower+diff/2
}
diff <- upper - lower
}
res[i] <- (lower+upper)/2
}
}
res
}

#' Rejection Sampling Algorithm
#'
#' @description Algorithm to sample from an arbitrary density function.
#' @param Fun cumulative distribution function
#' @param dens density
#' @param len number of samples
#' @param cand candidate area
#' @param grid fineness degree
#' @param method vectorial ("vector") or not ("free")
#' @references
#' Devroye, L. (1986). Non-Uniform Random Variate Generation. New York: Springer.
#' @examples
#' plot(density(RejSampling(dens = function(x) dnorm(x, 5, 1),
#'    len = 500, cand = seq(2, 9, by = 0.001), method = "free")))
#' lines(density(RejSampling(dens = function(x) dnorm(x, 5, 1), len = 500,
#'       cand = seq(2, 9, by = 0.001), method = "vector")), col=2)
#' curve(dnorm(x, 5, 1), from = 2, to = 8, add = TRUE, col = 3)
#' @export
RejSampling <- function(Fun, dens, len, cand, grid = 1e-03, method = c("vector", "free")){  # for negative support?!?
method <- match.arg(method)
if(method == "free"){
res <- numeric(len)
for(i in 1:len){
if(missing(cand)){
ca <- findCandidateArea(function(t) Fun(t))
}else{
ca <- range(cand)
}
mp <- optimize(f = function(t) dens(t), ca, maximum = TRUE)\$objective
resi <- NULL
while(is.null(resi)){
u <- runif(1,0,mp)
candi <- runif(1, ca[1], ca[2])
prob <- dens(candi)
if(u <= prob){
resi <- candi
}
}
res[i] <- resi
}
}
if(method == "vector"){
res <- numeric(len)
if(missing(cand)){
ca <- findCandidateArea(Fun)
cand <- seq(ca[1], ca[2], by = grid)
}
prob <- vapply(cand, function(v) dens(v), FUN.VALUE = numeric(1))
cand <- cand[prob != 0]
prob <- prob[prob != 0]
mp <- max(prob)
count <- 1
while(count <= len){
u <- runif(1, 0, mp)
ind <- sample(1:length(cand), 1)
if(u <= prob[ind]){
res[count] <- cand[ind]
count <- count + 1
}
}
}
return(res)
}

#' Adaptation of proposal standard deviation
#'
#' @description Adaptive MCMC: if acceptance rate of the chain is smaller than \code{lower} or larger than \code{upper},
#' the proposal standard deviation \code{propSd}=exp(l) is adapted with respect to function \code{delta.n},
#' that means, the new proposal standard deviation
#' is equal to exp(l-\code{delta.n(batch)}), respectively exp(l+\code{delta.n(batch)}).
#'
#' @param chain Markov chain
#' @param propSd current proposal standard deviation
#' @param batch number of batch (of chain)
#' @param lower lower bound
#' @param upper upper bound
#' @param delta.n function of batch number
#'
#' @return adapted proposal standard deviation
#' @references Rosenthal, J. S. (2011). Optimal Proposal Distributions and Adaptive MCMC. In: Handbook of Markov Chain Monte Carlo, pp. 93-112.
#' @export
ad.propSd <- function(chain, propSd, batch, lower = 0.3, upper = 0.6, delta.n = function(n) min(0.05, 1/sqrt(n))){
ar <- length(unique(chain))/length(chain)
lsi <- log(propSd)

lsi[ar < lower] <- lsi - delta.n(batch)
lsi[ar > upper] <- lsi + delta.n(batch)
exp(lsi)
}

findCandidateArea <- function(VFun, start = 1, pos.support = TRUE, quasi.null = 1e-05){
if(pos.support){
cand <- start
while(1 - VFun(cand) > quasi.null){
cand <- cand*2
}
upper <- cand

cand <- start
while(VFun(cand) > 0){
cand <- cand/2
}
lower <- cand
} else{
cand <- start
while(1 - VFun(cand) > quasi.null){
cand <- cand*2
}
upper <- cand

cand <- -start
while(VFun(cand) > quasi.null){
cand <- cand*2
}
lower <- cand
}

c(lower, upper)
}

#' Calculation of a proposal for burn-in phase and thinning rate
#'
#' @description The proposed burn-in is calculated by dividing the Markov chains into \code{m} blocks and calculate the 95\% credibility intervals and the respective mean.
#' Starting in the first one, the block is taken as burn-in as long as the mean of the current block is not in the credibility interval of the following block or vice versa.
#' The thinning rate is proposed by the first lag which leads to a chain autocorrelation less than \code{dependence}.
#' It is not easy to automate these choices, so it is highly recommended to verify the chains manually.
#' @param chain vector of Markov chain samples
#' @param dependence allowed dependence for the chain
#' @param m number of blocks
#' @return vector of burn-in and thinning
#' @export
diagnostic <- function(chain, dependence = 0.8, m = 10) {
lc <- length(chain)
K <- floor(lc/m)
thinning <- min(which(acf(chain[-(1:floor(lc/5))], plot=F)\$acf <= dependence)[1], floor(K/10), na.rm = TRUE)
he1 <- sapply(1:m, function(i) chain[((i - 1) * K + 1):(i * K)][seq(1, K, by = thinning)])

he2 <- apply(he1, 2, quantile, c(0.025, 0.975))
he.mean <- apply(he1, 2, mean)
is.in <- (he.mean[-m] >= he2[1, -1] & he.mean[-m] <= he2[2, -1]) | (he.mean[-1] >= he2[1, -m] & he.mean[-1] <= he2[2, -m])
#    burnIn <- 0
burnIn <- K
for (i in 1:(m-1)) {
if (sum(is.in) < length(is.in)) {
is.in <- is.in[-1]
burnIn <- burnIn + K
}
}
burnIn <- min((m-1)*K, burnIn)
return(c(burnIn = burnIn, thinning = thinning))
}
SimoneHermann/BaPreStoPro documentation built on May 10, 2017, 1:42 p.m.