Nothing
######################################################################
# Function to perform nowcast at a specific day "now" using a procedure
# which takes truncation of the available observations into
# account. The full documentation is available in the nowcast.Rd file.
#
# Author: Michael Hoehle <http://www.math.su.se/~hoehle>
#
# Parameters:
# now - a Date object representing today
# when - a vector of Date objects representing the days to do the forecast for.
# A requirement is that for all elements in when are smaller or equal
# than "now".
# data - the Database containing columns dEventCol and dReportCol, which
# contain the date of the event and of when the report arrives in
# the database.
# dEventCol - name of column in data containing time of event occurrence
# dReportCol - name of column in data containing time of reprt arrival
# method - which method to use
# D - maximum delay to consider
# m - moving window for delay estimation
# control - a list containing the following arguments
# * gd.prior.kappa - prior for delay is symmetric Dirichlet
# with concentration parameter gd.prior.kappa
#
# Note: As predictions are done simultaneously the entire vector of observations
# is casted. Then the subset specified in "when" is returned.
#
# Returns:
# stsNC object with reporting triangle, delay estimate and prediction interval in the appropriate slots.
#
# Todo:
# * yt.support to N.tInf support in nowcast??
# * bayes.notrunc and bayes.notrunc.bnb could become one code segment
# * Enable user to provide reporting triangle directly.
# * Function should work for weekly and monthly data as well
######################################################################
nowcast <- function(now, when, data, dEventCol="dHospital", dReportCol="dReport",
method=c("bayes.notrunc","bayes.notrunc.bnb","lawless","bayes.trunc","unif","bayes.trunc.ddcp"),
aggregate.by="1 day",
D=15, m=NULL, m.interpretation=c("hoehle_anderheiden2014", "lawless1994"),
control=list(
dRange=NULL,alpha=0.05,nSamples=1e3,
N.tInf.prior=c("poisgamma","pois","unif"),
N.tInf.max=300, gd.prior.kappa=0.1,
ddcp=list(ddChangepoint=NULL,
cp_order=c("zero","one"),
Wextra=NULL,
logLambda=c("iidLogGa","tps","rw1","rw2"),
responseDistr=c("poisson", "negbin"),
mcmc=c(burnin=2500,sample=10000,thin=1, adapt=1000, store.samples=FALSE)),
score=FALSE,predPMF=FALSE)) {
#Check if the runjags package is available (required for bayes.trunc.ddcp to work!
if ("bayes.trunc.ddcp" %in% method) {
if (!requireNamespace("runjags",quietly=TRUE)) {
stop("The \"bayes.trunc.ddcp\" method requires the runjags package to be installed, which is available from CRAN.")
}
}
if ((!inherits(now,"Date")) | (length(now)>1)) {
stop("The parameter 'now' has to be a single Date.")
}
#Check if all when_i<= now
if (!all(when<=now)) {
stop("Assertion when <= now failed.")
}
#Check that specified methods are all valid
method <- match.arg(method,c("bayes.notrunc","bayes.notrunc.bnb","lawless","bayes.trunc","unif","bayes.trunc.ddcp"),several.ok=TRUE)
######################################################################
# Time aggregation. Make sure it's a valid aggregational level and
# move all dates to the "first" of this level.
# @hoehle: Should work for day, weeks and month. Quarter and year not atm.
######################################################################
aggregate.by <- match.arg(aggregate.by,c("1 day","1 week", "1 month"),several.ok=FALSE)
epochInPeriodStr <- switch(aggregate.by, "1 day"="1","1 week"="%u", "1 month"="%d")
if (aggregate.by != "1 day") {
warning("Moving dates to first of each epoch.")
#Move dates back to first of each epoch unit
for (colName in c(dEventCol, dReportCol)) {
data[,colName] <- data[,colName] - as.numeric(format(data[,colName],epochInPeriodStr)) + 1
}
#Check now and when
if (!all( format( c(now,when),epochInPeriodStr) == 1)) {
stop("The variables 'now' and 'when' needs to be at the first of each epoch")
}
}
#Choose the correct difference function
if (aggregate.by == "1 day") {
timeDelay <- function(d1,d2) {as.numeric(d2-d1)}
}
if (aggregate.by == "1 week") {
timeDelay <- function(d1,d2) { floor(as.numeric(difftime(d2,d1,units="weeks"))) } #Count the number of full weeks
}
if (aggregate.by == "1 month") {
timeDelay <- function(d1,d2) {
#Helper function from http://stackoverflow.com/questions/1995933/number-of-months-between-two-dates
monnb <- function(d) {
lt <- as.POSIXlt(as.Date(d, origin="1900-01-01"))
lt$year*12 + lt$mon
}
monnb(d2) - monnb(d1) #count the number of full months
}
}
## Check the value of the m interpretation
m.interpretation <- match.arg(m.interpretation, c("hoehle_anderheiden2014", "lawless1994"))
if (!is.null(m) & (method == "lawless") & (m.interpretation != "lawless1994")) {
warning("Selected method is Lawless (1994), but the interpretation of m is a horizontal cut in the reporting triangle (as in Hoehle and an der Heiden (2014)) and not as in Lawless (1994).")
}
if (!is.null(m) & (method != "lawless") & (m.interpretation == "lawless1994")) {
stop("The selected nowcasting method only works with m.interpretation = 'hoehle_anderheiden2014'")
}
######################################################################
#If there is a specification of dateRange set dMin and dMax accordingly
#Otherwise use as limits the range of the data
######################################################################
if (is.null(control[["dRange",exact=TRUE]])) {
dMin <- min(data[,dEventCol],na.rm=TRUE)
dMax <- max(data[,dEventCol],na.rm=TRUE)
} else {
dMin <- control$dRange[1]
dMax <- control$dRange[length(control$dRange)]
}
#@hoehle - check that dRange is proper
if (!all( format( c(dMin, dMax), epochInPeriodStr) == 1)) {
stop("The variables in dRange needs to be at the first of each epoch.")
}
dateRange <- seq(dMin, dMax, by=aggregate.by)
######################################################################
# Additional manipulation of the control arguments
######################################################################
#Check if alpha is specified
if (is.null(control[["alpha",exact=TRUE]])) {
control$alpha <- 0.05
}
if (is.null(control[["N.tInf.prior",exact=TRUE]])) {
control$N.tInf.prior <- "unif"
}
if (is.null(control[["N.tInf.max",exact=TRUE]])) {
control$N.tInf.max <- 300
}
if (is.null(control[["gd.prior.kappa",exact=TRUE]])) {
control$gd.prior.kappa <- 0.1
}
if (is.null(control[["nSamples",exact=TRUE]])) {
control$nSamples <- 1e3
}
if (is.null(control[["score",exact=TRUE]])) {
control$score <- FALSE
}
#Checking for the bayes.trunc.ddcp procedure. If so make sure params are set up.
if ("bayes.trunc.ddcp" %in% method) {
#If no parameters at all set to defaults.
if (is.null(control[["ddcp",exact=TRUE]])) {
control$ddcp <- list(ddChangepoint=NULL,
cp_order="zero",
Wextra=NULL,
logLambda=c("iidLogGa","tps","rw1","rw2"),
tau.gamma=1,
response.distr=c("poisson"),
mcmc=c(burnin=2500, sample=10000, thin=1, adapt=1000, store.samples=FALSE))
}
#Check form of logLambda
if (is.null(control[["ddcp",exact=TRUE]][["logLambda",exact=TRUE]])) {
control[["ddcp"]] <- modifyList(control[["ddcp",exact=TRUE]], list(logLambda="iidLogGa"))
} else {
control[["ddcp"]]$logLambda <- match.arg(control[["ddcp"]][["logLambda"]],c("iidLogGa","tps","rw1","rw2"))
}
#Check breakpoint to use in case of bayes.trunc.ddcp (delay distribution with breakpoint)
if (is.null(control[["ddcp"]][["ddChangepoint"]]) ||
!inherits(control[["ddcp"]][["ddChangepoint"]], "Date")) {
stop("Please specify a Date object as changepoint in control$ddChangepoint.")
}
if (any(control[["ddcp"]][["ddChangepoint"]] > now)) {
warning("Some of the elements in ddChangepoint are beyond 'now'. This might be problematic!")
}
##Check cp_order variable
if (!is.null(control[["ddcp",exact=TRUE]][["cp_order",exact=TRUE]])) {
control[["ddcp"]]$cp_order <- match.arg(control[["ddcp"]][["cp_order"]],c("zero","one"))
} else {
control[["ddcp"]]$cp_order <- "zero"
}
##Check Wextra argument
if (!is.null(control[["ddcp",exact=TRUE]][["Wextra",exact=TRUE]])) {
if (!is.array(control[["ddcp",exact=TRUE]][["Wextra",exact=TRUE]])) {
stop("Wextra is not an array.")
}
}
#Make this an accessible variable
ddChangepoint <- control$ddcp$ddChangepoint
Wextra <- control$ddcp$Wextra
ncol_Wextra <- if (is.null(Wextra)) 0 else ncol(Wextra)
colnames_Wextra <- if (is.null(Wextra)) NULL else colnames(Wextra)
cp_order <- control$ddcp$cp_order
#Response distribution in the model
if (is.null(control[["ddcp",exact=TRUE]][["response.distr",exact=TRUE]])) {
control[["ddcp"]]$response.distr <- "poisson"
} else {
stopifnot(control[["ddcp",exact=TRUE]][["response.distr",exact=TRUE]] %in% c("poisson", "negbin"))
}
#Precision parameter for gamma coefficients for hazard delay distribution
if (is.null(control[["ddcp",exact=TRUE]][["tau.gamma",exact=TRUE]])) {
control[["ddcp"]]$tau.gamma <- 1
}
#Prior for eta ~ [eta.mu, eta.prec]
if (is.null(control[["ddcp",exact=TRUE]][["eta.mu",exact=TRUE]])) {
control[["ddcp"]]$eta.mu <- rep(0,length(ddChangepoint) + ncol_Wextra)
} else {
if (length(control[["ddcp"]]$eta.mu) != length(ddChangepoint) + ncol_Wextra) {
stop("length of eta.mu is different from the number of change points in 'ddChangepoint'.")
}
}
if (is.null(control[["ddcp",exact=TRUE]][["eta.prec",exact=TRUE]])) {
if (length(ddChangepoint) == 1) {
control[["ddcp"]]$eta.prec <- 1
} else {
control[["ddcp"]]$eta.prec <- diag(rep(1, length(ddChangepoint)))
}
} else {
#Custom option
if (length(ddChangepoint) == 1) {
if (length(control[["ddcp"]]$eta.prec) != 1) {
stop("length of eta.prec is different from the number of change points in 'ddChangepoint'.")
} else {
if (!( (nrow(control[["ddcp"]]$eta.prec) == (ncol(control[["ddcp"]]$eta.prec))) & (nrow(control[["ddcp"]]$eta.prec) == length(ddChangepoint)))) {
stop(paste0("dimension ", dim(control[["ddcp"]]$eta.prec), " of eta.prec is different from the number of change points in 'ddChangepoint' (",length(ddChangepoint),"."))
}
}
}
}
#Check MCMC options
if (is.null(control[["ddcp",exact=TRUE]][["responseDistr",exact=TRUE]])) {
control[["ddcp"]][["responseDistr"]] <- "poisson"
} else {
# Check that it's a valid response distribution
stopifnot(control[["ddcp"]][["responseDistr"]] %in% c("negbin","poisson"))
}
#Check MCMC options
if (is.null(control[["ddcp",exact=TRUE]][["mcmc",exact=TRUE]])) {
control[["ddcp"]][["mcmc"]] <- c(burnin=2500,sample=10000,thin=1, adapt=1000, store.samples=FALSE)
} else {
if (!all(c("burnin","sample","thin", "adapt", "store.samples") %in% names(control[["ddcp",exact=TRUE]][["mcmc",exact=TRUE]]))) {
stop("mcmc option list needs names 'burnin', 'sample', 'thin', 'adapt' and 'store.samples'.")
}
}
#done with options for bayes.trunc.ddcp
}
######################################################################
# Do preprocessing of the data
######################################################################
hasNADates <- is.na(data[,dEventCol]) | is.na(data[,dReportCol])
data <- data[!hasNADates,]
message(paste0("Removed ",sum(hasNADates), " records due to NA dates."))
#Create a column containing the reporting delay using the timeDelay
#function
data$delay <- timeDelay(data[,dEventCol],data[,dReportCol])
#Handle delays longer than D.
#@hoehle - handle that the unit might not just be days
#notThereButDThere <- (data[,dReportCol] > now) & ((data[,dEventCol]) + D <= now)
notThereButDThere <- (timeDelay(data[,dReportCol],now) < 0) & (timeDelay(data[,dEventCol],now) >= D)
if (sum(notThereButDThere,na.rm=TRUE)) {
warning(paste(sum(notThereButDThere,na.rm=TRUE), " observations > \"now\" due to a delay >D. If delay cut to D they would be there."),sep="")
}
#Which observations are available at time s
#@hoehle: data.sub <- data[ na2FALSE(data[,dReportCol] <= now),]
data.sub <- data[ na2FALSE(timeDelay(data[,dReportCol],now) >= 0),]
if (nrow(data.sub)==0) {
stop(paste("No data available at now=",now,"\n"))
}
#Create an sts object containing the observed number of counts until s
sts <- linelist2sts(data.sub,dEventCol,aggregate.by=aggregate.by,dRange=dateRange)
sts <- as(sts,"stsNC")
#Create an extra object containing the "truth" based on data
sts.truth <- linelist2sts(data, dEventCol, aggregate.by=aggregate.by, dRange=dateRange)
#List of scores to calculate. Can become an argument later on
scores <- c("logS","RPS","dist.median","outside.ci")
#Initialize scoring rule results - to be saved in control slot -- dirty
SR <- array(0,dim=c(nrow(sts),length(method),length(scores)))
#List for storing the predictive PMFs.
if (is.null(control[["predPMF",exact=TRUE]])) {
control$predPMF <- FALSE
}
#Prepare a list of different estimated of the delay CDF
delayCDF <- list()
######################################################################
# Done manipulating the control list with default arguments
######################################################################
sts@control <- control
#Save truth
sts@truth <- sts.truth
#Reserve space for returning the predictive PMFs
sts@predPMF <- list()
######################################################################
# Consistency checks
######################################################################
#Check if support of N.tInf is large enough
if (2*control$N.tInf.max < max(observed(sts),na.rm=TRUE)) {
warning("N.tInf.max appears too small. Largest observed value is more than 50% of N.tInf.max, which -- in case this number is extrapolated -- might cause problems.\n")
}
#Create a vector representing the support of N.tInf
N.tInf.support <- 0:control$N.tInf.max
#======================================================================
#======================================================================
# Build reporting triangle and derived parameters for delay
#======================================================================
#======================================================================
cat("Building reporting triangle...\n")
#Time origin t_0
t0 <- min(dateRange)
#Sequence from time origin until now (per day??)
t02s <- seq(t0, now,by=aggregate.by)
#Maximum time index
T <- length(t02s)-1
#Check if the maximum delay is longer than the available time series
if (D>T) {
stop("D>T. Cannot estimate the long delays.")
}
#How many observations to take for estimating the delay distribution
if (is.null(m)) {
m <- T
}
if (m<1) { stop("Assertion m>=1 not fulfilled.") }
#Define the observation triangle
n <- matrix(NA,nrow=T+1, ncol=T+1, dimnames=list(as.character(t02s), NULL))
#Loop over time points. (more efficient than delay and then t)
for (t in 0:T) {
#Extract all reports happening at time (index) t.
#@hoehle: data.att <- data.sub[na2FALSE(data.sub[,dEventCol] == t02s[t+1]), ]
data.att <- data.sub[na2FALSE(timeDelay(data.sub[,dEventCol], t02s[t+1])) == 0, ]
#Loop over all delays
for (x in 0:(T-t)) {
#Count number with specific delay
n[t+1,x+1] <- sum(data.att[,"delay"] == x)
}
}
cat("No. cases: ",sum(n,na.rm=TRUE),"\n")
#Handle delays longer than D
#@hoehle: Not done! Just fix them to have delay D.
# browser()
nLongDelay <- apply(n[,(D+1)+seq_len(T-D), drop=FALSE],1,sum,na.rm=TRUE)
if (any(nLongDelay>0)) {
warning(paste(sum(nLongDelay)," cases with a delay longer than D=",D," days forced to have a delay of D days.\n",sep=""))
n <- n[,1:(D+1)]
n[,(D+1)] <- n[,(D+1)] + nLongDelay
} else {
#No problems. Just extract up to D+1
n <- n[,1:(D+1)]
}
#Calculate n.x and N.x as in (2.7) and (2.8) and Fig.2 of Lawless (1994)
#Note the different moving window definition as in the Lawless article.
n.x <- rep(0,times=D+1)
N.x <- rep(0,times=D+1)
##Compute n.x and N.x
for (x in 0:D) {
##Define time of occurrence sliding window index set (see documentation)
if (m.interpretation == "hoehle_anderheiden2014") {
toc_index_set <- max(0,T-m):(T-x)
} else { #hoehle: Lawless definition is max(0,T-m-x)
toc_index_set <- max(0,T-m-x):(T-x)
}
## Count
for (t in toc_index_set) {
#cat("x=",x,"\tt=",t,":\n")
n.x[x+1] <- n.x[x+1] + n[t+1,x+1]
for (y in 0:x) {
#cat("x=",x,"\tt=",t,"\ty=",y,":\n")
N.x[x+1] <- N.x[x+1] + n[t+1,y+1]
}
}
}
cat("No. cases within moving window: ",sum(n.x,na.rm=TRUE),"\n")
#Available observations at time T, definition of N(t;T) on p.17.
N.tT <- sapply(0:T, function(t) sum(n[t+1, 0:min(D+1,(T-t)+1)]))
#Truth - already in another object. Delete??
N.tInf <- table( factor(as.character(data[,dEventCol]),levels=as.character(t02s)))
#Store results of the reporting triangle in the control slot together with additional
#attributes for fast access of, e.g., summaries or defining variables.
reportingTriangle <- n
attr(reportingTriangle, "n.x") <- n.x
attr(reportingTriangle, "N.x") <- N.x
attr(reportingTriangle, "N.tT") <- N.tT
attr(reportingTriangle, "N.tInf") <- N.tInf
attr(reportingTriangle, "T") <- T
attr(reportingTriangle, "D") <- D
attr(reportingTriangle, "t02s") <- t02s
sts@reportingTriangle <- reportingTriangle
#======================================================================
# Calculations are jointly for all t values.
#======================================================================
#List of casts each containing a table 0..N.tInf.max with the PMF
Ps <- list()
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
# Lawless (1994) method without adjustment for overdispersion
#
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ("lawless" %in% method) {
#Hazard function estimates, i.e. g-function estimate as in (2.9)
#of Lawless (1994). NAs are set to zero (consequences??)
g.hat <- ifelse( !is.na(n.x/N.x), n.x/N.x, 0)
#Force g.hat(0)=1 as stated just below (2.1)
g.hat[1] <- 1
#Check how the estimated CDF looks
#F <- NULL ; for (d in 0:D) { i <- d+seq_len(D-d) ; F[d+1] <- prod(1-g.hat[i+1]) }
#plot(0:D,F)
#Compute weights Wt.hat as in eqn. (2.13) of Lawless (1994). Use T1=Inf.
#Note: Wt.hat estimates F_t(T-t).
T1 <- Inf
What.t <- sapply(0:T, function(t) {
if (t<T-D) {
1
} else {
x = T-t + seq_len(min(T1-t,D)-(T-t)) #(2.13) with modifications. knifflig.
prod(1-g.hat[x+1])
}
})
#Store result of delay distribution estimation for support 0:T
delayCDF[["lawless"]] <- rev(What.t)
#Do the prediction as in (2.12)
Nhat.tT1 <- N.tT / What.t
#V.Wt as in (2.15). Note: only depends implicitly on m through N.x.
Vhat.Wt <- sapply(0:T, function(t) {
if (t<T-D) {
0
} else {
x = T-t + seq_len(min(T1-t,D)-(T-t))
What.t[t+1]^2 * sum( g.hat[x+1]/(N.x[x+1]*(1-g.hat[x+1])),na.rm=TRUE)
}
})
#(2.16).
Vhat.Zt <- sapply(0:T, function(t) {
if (t<T-D) {
0
} else {
(N.tT[t+1]*(N.tT[t+1]+1))/What.t[t+1]^4*Vhat.Wt[t+1]
+ N.tT[t+1]*(1-What.t[t+1])/What.t[t+1]^2
}
})
#Upper and lower 95% prediction limits based on the asymptotic normal
#conf.level <- 0.95
#U <- Nhat.tT1 + qnorm( (1+conf.level)/2)* sqrt(Vhat.Zt)
#L <- pmax(N.tT, Nhat.tT1 - qnorm( (1+conf.level)/2)* sqrt(Vhat.Zt))
#Discretize result: all mass below actual observation is put to observation
PMFs <- matrix(NA, nrow=control$N.tInf.max+1,ncol=T+1)
#CDF of a left truncated normal, truncated at "at"
ltruncpnorm <- function(x, mean, sd, at) {
ifelse( x < at, 0, pnorm(x, mean, sd) / (1-pnorm(at, mean,sd)))
}
#Deduce PMF for each time point in the past
for (i in 1:length(Nhat.tT1)) {
#Safeguard against NAs
if (is.na(Vhat.Zt[i])) { warning("Vhat.Zt[i] is NA") ; Vhat.Zt[i] <- -1e99 }
#If Vhat.Zt = 0 then Nhat.tT1 \equiv 0! Special care needed here
if (Vhat.Zt[i] > 0) {
CDF <- c(0,ltruncpnorm(N.tInf.support, mean=Nhat.tT1[i], sd=sqrt(Vhat.Zt[i]),at=N.tT[i]))
PMFs[,i] <- diff(CDF)
} else {
PMFs[,i] <- (N.tInf.support == Nhat.tT1[i])*1
}
}
Ps[["lawless"]] <- PMFs
} #end lawless procedure
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
# Bayesian method (simple model, clever sampling -> no MCMC)
#
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#Part jointly for both bayes and bayes.notrunc
if (("bayes.trunc" %in% method) | ("bayes.notrunc" %in% method)) {
cat("bayes prep...\n")
######################################################################
# Prior of N(t,\infty)
######################################################################
N.tInf.prior <- control$N.tInf.prior
#Extract prior parameters from prior choice
if (N.tInf.prior == "pois") {
lambda <- attr(N.tInf.prior,"lambda",exact=TRUE)
} else {
if (N.tInf.prior == "poisgamma") {
#Find size parameters such that mean variance is as target.
var.prior <- function(size.prior) { mean.prior + mean.prior^2/size.prior }
#If mean & var specified
if (all(c("mean.lambda","var.lambda") %in% names(attributes(N.tInf.prior)))) {
mean.prior <- attr(N.tInf.prior,"mean.lambda",exact=TRUE)
var.prior.target <- attr(N.tInf.prior,"var.lambda",exact=TRUE)
size.prior <- uniroot( function(size.prior) { var.prior(size.prior) - var.prior.target},interval=c(1e-12,50))$root
#Check result
cat("(E,V) of prior for lambda = (",paste(c(mean.prior,var.prior(size.prior)),collapse=","),")\n")
} else {
stop("mean.lambda and var.lambda not part of prior specification")
}
} else {
if (N.tInf.prior == "unif") {
N.tInf.prior.max <- attr(N.tInf.prior,"N.tInf.prior.max",exact=TRUE)
} else {
#No option applied
stop("Not a valid prior!")
}
}
}
######################################################################
# Define function to generate PMF for max(0,T-D),..,T by sampling.
#
# Parameters:
# alpha.star, beta.star - vector containing the posterior GD params
######################################################################
pmfBySampling <- function(alpha.star, beta.star) {
#Sample from posterior distribution, i.e. sample from the reverse distribution
#and reverse result
p.sample <- rgd(control$nSamples,alpha.star,beta.star)[,(length(alpha.star)+1):1]
#All the time points where extrapolation is to be done
tSet <- max(0,(T-D)):T
######################################################################
# Procedure to generate nowcasts of all time points up to T-D,...,T.
# This is based on the posterior samples available in p.sample.
# Current code adds up the PMF tables instead of a pure sample based
# procedure and also prevents PMF=0 better than tabulating the samples.
######################################################################
N.tT1.pred <- array(NA, dim=c(dim(p.sample)[1],control$N.tInf.max+1,dim(p.sample)[2]),dimnames=list(NULL,seq_len(control$N.tInf.max+1)-1L,tSet))
for (j in 1:control$nSamples) {
#Extract delay PMF from sample
p <- p.sample[j,]
#Proportion reported up to x, x=0,..,T
F <- c(rep(1,T-D),rev(cumsum(p)))
#Guard against numerical instability: ensure that not larger than 1.
F <- ifelse(F>1,1,F)
#Loop over all time points to nowcast
for (i in 1:length(tSet)) {
t <- tSet[i]
N.tT1.pred[j,,i] <- switch(N.tInf.prior,
"poisgamma"=dpost.bnb(N.tT[t+1],sumpd=F[t+1],mu=mean.prior,size=size.prior,N.tInf.max=control$N.tInf.max))
}
}
#Average the PMFs as in Step (2) of the algorithm
PMF <- apply(N.tT1.pred,MARGIN=c(2,3),mean)
#Add part, where no prediction needs to be done
if (T-D>0) {
#Empty PMFs
determined <- matrix(0,nrow=control$N.tInf.max+1,ncol=T-D-1+1)
#Add "1" entry at the observed
for (t in 0:(T-D-1)) {
determined[N.tT[t+1]+1,t+1] <- 1
}
PMF <- cbind(determined,PMF)
}
return(PMF)
} #done definition of pmfBySampling
}
if ("bayes.trunc" %in% method) {
cat("bayes.trunc...\n")
######################################################################
#Prior of reporting delay as parameters of generalized Dirichlet prior
######################################################################
#Define symmetric dirichlet as prior, just as in the other case
alpha.prior <- rep(control$gd.prior.kappa, D)
beta.prior <- rep(0,D)
beta.prior[D] <- control$gd.prior.kappa
for (i in (D-1):1) {
beta.prior[i] <- alpha.prior[i+1] + beta.prior[i+1]
}
######################################################################
# Posterior section
######################################################################
#Deduce posterior distribution of delay distribution, i.e. it is again
#a generalized Dirichlet
alpha <- beta <- rep(NA,D)
for (d in 0:(D-1)) {
alpha[d+1] <- n.x[D-d+1] ##Note: +1 coz index 1 is delay 0.
beta[d+1] <- N.x[D-d+1] - n.x[D-d+1]
}
#Check if there are any points without data and warn about it.
if (any(alpha + beta == 0)) {
warning("The delays ",paste(D-which(alpha+beta==0)-1,collapse=",")," have no observations. Results might be instable and depend all on prior.")
}
#Add up. Note: Delay zero (i.e. element D+1) is ignored as this is
#not modelled explicitly by the GD distribution (sum to 1 constraints)
alpha.star <- alpha.prior + alpha
beta.star <- beta.prior + beta
#Compute the expectation of the GD distribution and store this as the delay
delayCDF[["bayes.trunc"]] <- cumsum(rev(Egd(alpha.star,beta.star)))
#Save result
Ps[["bayes.trunc"]] <- pmfBySampling(alpha.star, beta.star)
} # end "bayes.trunc" %in% method
#======================================================================
# Bayesian version which ignores truncation
#======================================================================
if ("bayes.notrunc" %in% method) {
cat("bayes.notrunc...\n")
######################################################################
# Prior section
######################################################################
alpha.prior <- rep(control$gd.prior.kappa, D) #symmetric dirichlet
beta.prior <- rep(0,D)
beta.prior[D] <- control$gd.prior.kappa
for (i in (D-1):1) {
beta.prior[i] <- alpha.prior[i+1] + beta.prior[i+1]
}
######################################################################
# Posterior section
######################################################################
#Deduce posterior distribution of delay distribution, i.e. it is again
#a generalized Dirichlet
alpha <- beta <- rep(NA,D)
for (d in 0:(D-1)) {
alpha[d+1] <- n.x[D-d+1]
beta[d+1] <- sum(n.x[D - (d+1):D + 1])
}
#Check if there are any points without data and warn about it.
if (any(alpha + beta == 0)) {
warning("The delays ",paste(D-which(alpha+beta==0)-1,collapse=",")," have no observations. Results might be instable and depend all on prior.")
}
#Posterior parameters.
alpha.star <- alpha.prior + alpha
beta.star <- beta.prior + beta
#Check that its a ordinary Dirichlet
for (i in (D-1):1) {
if (!all.equal(beta.star[i], alpha.star[i+1] + beta.star[i+1])) {
warning("Posterior at i=",i," is not an ordinary Dirichlet as it's supposed to be.")
}
}
#Save resulting delay distribution
delayCDF[["bayes.notrunc"]] <- cumsum(rev(Egd(alpha.star,beta.star)))
Ps[["bayes.notrunc"]] <- pmfBySampling(alpha.star,beta.star)
} # end bayes.notrunc
######################################################################
# Unadjusted procedure using beta negative binomial. ToDo:
# integrate code with other Bayesian procedures
######################################################################
if ("bayes.notrunc.bnb" %in% method) {
cat("bayes.notrunc.bnb...\n")
######################################################################
# Prior section (same as for all methods)
######################################################################
alpha.prior <- rep(control$gd.prior.kappa, D) #symmetric dirichlet
beta.prior <- rep(0,D)
beta.prior[D] <- control$gd.prior.kappa
for (i in (D-1):1) {
beta.prior[i] <- alpha.prior[i+1] + beta.prior[i+1]
}
######################################################################
# Posterior section
######################################################################
#Deduce posterior distribution of delay distribution, i.e. it is again
#an ordinary Dirichlet
alpha <- beta <- rep(NA,D)
for (d in 0:(D-1)) {
alpha[d+1] <- n.x[D-d+1]
beta[d+1] <- sum(n.x[D - (d+1):D + 1])
}
#Check if there are any points without data and warn about it.
if (any(alpha + beta == 0)) {
warning("The delays ",paste(D-which(alpha+beta==0)-1,collapse=",")," have no observations. Results might be instable and depend all on prior.")
}
#Posterior parameters.
alpha.star <- alpha.prior + alpha
beta.star <- beta.prior + beta
#Check that its a ordinary Dirichlet
for (i in (D-1):1) {
if (!all.equal(beta.star[i], alpha.star[i+1] + beta.star[i+1])) {
warning("Posterior at i=",i," is not an ordinary Dirichlet as it's supposed to be.")
}
}
#Save resulting delay distribution (i.e. no truncation adjustment)
delayCDF[["bayes.notrunc"]] <- cumsum(rev(Egd(alpha.star,beta.star)))
#Allocate PMF to return
PMF <- matrix(0,nrow=control$N.tInf.max+1,ncol=length(max(0,(T-D)):T))
#Concentration parameter vector of the ordinary Dirichlet distribution
#Note. alpha.star vector is reversed (shortest delay last).
alpha <- rev(c(alpha.star,beta.star[length(beta.star)]))
#consistency check
if (!all.equal(rev(Egd(alpha.star,beta.star)),alpha/sum(alpha))) {
stop("Problem. GD and Dirichlet do not correspond...")
}
tSet <- max(0,(T-D)):T
for (i in 1:length(tSet)) {
t <- tSet[i]
alpha.i <- cumsum(alpha)[T-t+1]
beta.i <- sum(alpha) - alpha.i
if (T-t==D) {
PMF[,i] <- ifelse( N.tInf.support == N.tT[t+1], 1, 0)
} else {
#Calculate PMF knowing the q ~ Beta( , ) by the aggregation
#property.
#Note: Vector N.tT starts at time zero, i.e. time T corresponds to T+1
PMF[,i] <- dbnb( N.tInf.support - N.tT[t+1],n=N.tT[t+1]+1, alpha=alpha.i, beta=beta.i)
}
} #done looping over all time points
#Add part, where no prediction needs to be done
if (T-D>0) {
#Empty PMFs
determined <- matrix(0,nrow=control$N.tInf.max+1,ncol=T-D-1+1)
#Add "1" entry at the observed
for (t in 0:(T-D-1)) {
determined[N.tT[t+1]+1,t+1] <- 1
}
PMF <- cbind(determined,PMF)
}
Ps[["bayes.notrunc.bnb"]] <- PMF
} # end bayes.notrunc.bnb
######################################################################
# Fully Bayes version with MCMC
######################################################################
if ("bayes.trunc.ddcp" %in% method) {
#Allocate result
PMF <- matrix( 0,ncol=(T+1),nrow=control$N.tInf.max+1)
#Fix seed value of JAGS RNG for each chain
n.chains <- 3
init <- lapply(1:n.chains,function(i) {
list(.RNG.name="base::Mersenne-Twister",.RNG.seed=i*10)
})
#Make design matrix for a quadratic TPS spline in time
makeTPSDesign <- function(T,degree=2) {
nbeta=degree + 1
X <- matrix(NA,ncol=nbeta, nrow=T+1)
for (t in 0:T) {
#Form a centered time covariate
t.centered <- t - T/2
for(pow in 0:degree) {
X[t+1,pow+1]<- t.centered^(pow)
}
}
#Make the knot points evenly spaced between 0,T not including these points
knots <- seq(0,T,length=min(round(T/6)+2,22))
knots <- knots[-c(1,length(knots))]
#Remove knots which are beyond T-maxDelay/2
knots <- knots[knots <= T-D/2]
knots <- knots - T/2
nknots <- length(knots)
#Penalty as REs - setup design matrix
Z <- matrix(NA,nrow=T+1,ncol=length(knots))
for (t in 0:T){
t.center <- t - T/2
for(k in 1:nknots){
Z[t+1,k]<- pmax((t.center-knots[k]),0)^degree
}
}
return(list(X=X,Z=Z,knots=knots,nknots=nknots,nbeta=nbeta))
}
tps <- makeTPSDesign(T=T,degree=2)
##browser()
#Design matrix for logistic discrete time hazard model containing
#changepoints. Could be extended s.t. the user provides W.
W <- array(NA, dim=c(T+1, length(ddChangepoint) + ncol_Wextra, D+1), dimnames=list(as.character(t02s), c(as.character(ddChangepoint),colnames_Wextra),paste("delay",0:D,sep="")))
for (t in 0:T) {
for (i in 1:length(ddChangepoint)) {
# Shape design matrix for change-points
if (cp_order == "zero") {
W[t+1, i, ] <- as.numeric( (t02s[t+1] + 0:D) >= ddChangepoint[i])
} else if (cp_order == "one") {
W[t+1, i, ] <- pmax(0, as.numeric( (t02s[t+1] + 0:D) - ddChangepoint[i]))
}
}
# Add additional effects as part of the design matrix
if (ncol_Wextra > 0) {
W[t + 1, length(ddChangepoint) + 1:ncol(Wextra), ] <- Wextra[t+1,, ]
}
}
#Priors. Uniform on the delays
D.prime <- round( D/2-0.4)+1
p.prior <- rep(1/(D.prime+1), D.prime+1)
mu.gamma <- qlogis( p.prior[1])
for (d in 1:(D.prime-1)) {
mu.gamma <- c(mu.gamma, qlogis( p.prior[d+1] / (1-sum(p.prior[1:d]))))
}
tau.gamma <- rep(control$ddcp$tau.gamma,times=D.prime)
#Prepare data for JAGS
jagsData <- list(#Data
rT=n,T=T+1,m=m+1,maxDelay=D,
#Time dependent logistic discrete hazard model
W=W, eta.mu=control$ddcp$eta.mu, eta.prec=control$ddcp$eta.prec,
mu.gamma=mu.gamma, tau.gamma=tau.gamma,
#Epidemic curve
alpha.lambda=2500/3000,beta.lambda=50/3000,
#Spline related stuff
X=tps$X,Z=tps$Z,nknots=tps$nknots,beta.mu=rep(0,tps$nbeta),beta.prec=1e-6*diag(tps$nbeta)
)
#Select appropriate model (one of: "tps","rw2","iid", "rw" as specified in the options)
logLambda.method <- control$ddcp$logLambda
responseDistr<- control$ddcp$responseDistr
### browser()
#Load the BUGS specification of the Bayesian hierarchical Poisson model
bugsModel <- readLines(file.path(path.package('surveillance'), 'jags',"bhpm.bugs"))
# Load local file
#bugsModel <- readLines(file.path("bhpm.bugs"))
bugsModel <- gsub(paste("#<",logLambda.method,">",sep=""),"",bugsModel)
bugsModel <- gsub(paste("#<",responseDistr,">",sep=""),"",bugsModel)
##browser()
#Problem when eta is scalar (TODO: Improve the solution.)
if (length(ddChangepoint) == 1) {
#Make eta ~ dnorm( , ) instead of eta ~ dmnorm
bugsModel <- gsub("(^[ ]*eta ~ )(dmnorm)","\\1dnorm",bugsModel)
#Use eta[1] instead of eta for matrix multiplication
bugsModel <- gsub("(eta)(.*%\\*%)","eta\\[1\\]\\2",bugsModel)
}
#cat(paste(bugsModel,collapse="\n"))
bugsFile <- tempfile(pattern = "nowcast-")
writeLines(bugsModel, bugsFile)
##browser()
## if (FALSE) {
## #Try to compile the model with ordinary rjags to see if there are any problems
## #before doing 3 chains parallelized using runjags.
## model <- jags.model(bugsFile,
## data = jagsData,
## init=init, #Fix seed value of JAGS as well
## n.chains = n.chains, n.adapt = 100)
## list.samplers(model)
## coda.samples(model,variable.names='logLambda',n.iter=100)
## }
######################################################################
# runjags way -- ToDo: parametrize using control options!
######################################################################
runjagsMethod <- 'rjparallel' #'rjags'
monitor <- c('gamma','eta','logLambda','NtInf', ifelse(control$ddcp$responseDistr == "negbin", "r", NA))
samples.rj <- runjags::run.jags(bugsFile,#bugsModel,
monitor = monitor, data=jagsData, n.chains=3,
inits = init,
burnin = control$ddcp$mcmc["burnin"],
sample = control$ddcp$mcmc["sample"],
thin = control$ddcp$mcmc["thin"],
adapt = control$ddcp$mcmc["adapt"],
summarise = FALSE,
method=runjagsMethod)
#Extract posterior median of discrete survival time delay distribution model parameters
dt.surv.samples <- coda::as.mcmc.list(samples.rj, vars = c('gamma','^eta'))
post.median <- dt.surv.pm <- apply( as.matrix(dt.surv.samples), 2, median)
#Posterior median of the lambda's
lambda.post <- exp(apply( as.matrix(coda::as.mcmc.list(samples.rj, vars = c('^logLambda'))), 2,
quantile, prob=c(0.025,0.5,0.975)))
#Extract posterior median of model parameters
gamma.red <- post.median[grep("gamma",names(post.median))]
eta <- matrix(post.median[grep("^eta",names(post.median))])
rownames(eta) <- colnames(W)
#Map from reduced set to full set
gamma <- gamma.red[round( (0:(D-1))/2 - 0.4) + 1]
#browser()
#Compute the resulting PMF from the model. Possibly put this in a separate function.
pmf <- matrix(NA, nrow=nrow(W), ncol=D+1, dimnames=list(as.character(t02s), paste("delay", 0:D, sep="")))
#Determine PMF
for (t in 1:length(t02s)) {
if (as.character(t02s[t]) %in% rownames(W)) {
lin.pred <- ( gamma + t(eta) %*% W[t,,0:D])
pmf[t,] <- haz2pmf(c(plogis(lin.pred),1))
}
}
#Store result as CDF
delayCDF[["bayes.trunc.ddcp"]] <- t(apply(pmf, 1, cumsum))
# Convert to coda compatible output.
samples <- coda::as.mcmc.list(samples.rj)
# Store model as attribute
if(control$ddcp$logLambda != "tps") tps <- NULL
# Configure list with model output and store is as an attribute.
list_return <- list(post.median=dt.surv.pm,W=W,lambda.post=lambda.post,tps=tps, gamma=gamma, eta=eta)
if (control[["ddcp",exact=TRUE]][["mcmc",exact=TRUE]][["store.samples", exact=TRUE]]) {
list_return <- modifyList(list_return, list(mcmc_samples = samples))
}
attr(delayCDF[["bayes.trunc.ddcp"]],"model") <- list_return
# Extract PMFs
for (t in 0:T) {
#Extract samples related to this time point
vals <- as.matrix(samples[,paste("NtInf[",t+1,"]",sep="")])
#PMF
PMF[,t+1] <- prop.table(table(factor(vals,levels=0:control$N.tInf.max)))
}
Ps[["bayes.trunc.ddcp"]] <- PMF
}
#======================================================================
#A really bad forecast -- the uniform
#======================================================================
if ("unif" %in% method) {
#Allocate result
PMF <- matrix( 0,ncol=(T+1),nrow=control$N.tInf.max+1)
#Loop over all time points to nowcast and put U(N.tT[t],Nmax)
for (t in 0:T) {
#How many values are there in N.tT .. Nmax
noVals <- max(0,control$N.tInf.max - N.tT[t+1]) + 1
#PMF at t is 0,...0 (N.tT-1 times), 1/noVals,...,1/noVals
PMF[,t+1] <- c(rep(0,N.tT[t+1]),rep(1/noVals,times=noVals))
}
Ps[["unif"]] <- PMF
}
######################################################################
#Loop over all time points in the vector "when". Only these are
#returned.
######################################################################
idxt <- which(dateRange %in% when)
for (i in idxt) {
#Save PMFs if that is requested
if (control$predPMF) {
res <- list()
for (j in 1:length(method)) {
res[[method[j]]] <- Ps[[method[j]]][,i]
}
sts@predPMF[[as.character(dateRange[i])]] <- res
}
#Evaluate scoring rules, if requested
if (control$score) {
#Infer the true value
ytinf <- observed(sts.truth)[i,]
#Evaluate all scores for all predictive distributions
for (i.P in 1:length(method)) {
for (i.score in 1:length(scores)) {
#cat("i=",i," i.P=",i.P," (",method[i.P],") i.score=",i.score,"\n")
SR[i,i.P,i.score] <- do.call(scores[i.score],args=list(P=Ps[[method[i.P]]][,i],y=ytinf,alpha=control$alpha))
}
}
} #end if control$score
#Add first nowcast & ci to stsNC slots
sts@upperbound[i,] <- median(N.tInf.support[which.max( cumsum(Ps[[method[1]]][,i])>0.5)])
sts@pi[i,,] <- N.tInf.support[c(which.max(cumsum(Ps[[method[1]]][,i]) > control$alpha/2),which.max(cumsum(Ps[[method[1]]][,i]) > 1-control$alpha/2))]
dimnames(sts@pi) <- list(as.character(dateRange),NULL,paste( c(control$alpha/2*100,(1-control$alpha/2)*100),"%",sep=""))
} #end of loop over time points
#Add scoring rule to output
if (control$score) {
dimnames(SR) <- list(as.character(dateRange),method,scores)
sts@SR <- SR
}
######################################################################
#Other arguments to save in control object
######################################################################
sts@control$N.tInf.support <- N.tInf.support
sts@control$method <- sts@control$name <- method
#Store variables relevant for the nowcast
sts@control$D <- D
sts@control$m <- m
sts@control$now <- now
sts@control$when <- when
sts@control$timeDelay <- timeDelay
#Store delayCDF object
sts@delayCDF <- delayCDF
#For backwards compatibility -- change this in the future TODODODODODO!
sts@control$yt.support <- sts@control$N.tInf.support
sts@control$y.prior.max <- sts@control$N.tInf.max
##Store the call options
theCall <- list(now=now,when=when,data=data,dEventCol=dEventCol,dReportCol=dReportCol,method=method,aggregate.by=aggregate.by,D=D, m=m)
sts@control$call <- theCall
##Done
return(sts)
}
######################################################################
# Helper functions
######################################################################
#Helper function
na2FALSE <- function(x) {x[is.na(x)] <- FALSE ; return(x) }
######################################################################
# Logarithmic score
#
# Parameters:
# P - predictive distribution, given as a vector containing the PMF
# with support 0,...,N.prior.max
# y - the actual observation. Can be a vector.
#
# Returns:
# -log P(y). If y outside 0,..,N.prior.max then -Inf.
######################################################################
logS <- function(P, y, ...) {
return(ifelse( y>=0 & y<=length(P)-1, -log(P[y+1]), -Inf))
}
######################################################################
# Ranked probability score
#
# Parameters:
# P - predictive distribution, given as a vector containing the PMF
# with support 0,...,N.prior.max
# y - the actual observation. Can be a vector.
#
# Returns:
# -log P(y). If y outside 0,..,N.prior.max then -Inf.
######################################################################
RPS <- function(P,y, ...) {
N.support <- 0:(length(P)-1)
sum( (cumsum(P) - (y <= N.support))^2)
}
#Some other scoring rules which are not proper.
dist.median <- function(P,y, ...) {
point.estimate <- which.max(cumsum(P)>=0.5) - 1
return(abs(point.estimate - y))
}
#0/1 indicator of observed value outside equal tailed (1-alpha/2) CI
outside.ci <- function(P,y,alpha) {
N.support <- 0:(length(P)-1)
ci <- N.support[c(which.max(cumsum(P) > alpha/2),which.max(cumsum(P) >
1-alpha/2))]
ifelse( y>=ci[1] & y<=ci[2], 0, 1)
}
######################################################################
# Helper functions for sampling the predictive distribution
######################################################################
#Unnormalized in Binomial-Negative-Binomial Hierarchy. Should work for vectors of N.tInf!
#Only kernel parts for N.tInf need to be taken into account
dpost.bnb.unorm <- function(N.tInf, N.tT, sumpd, mu, size) {
dbinom(N.tT, size=N.tInf, prob=sumpd)*dnbinom(N.tInf, mu=mu,size=size)
#Direct implementation - appears to be less stable...
#ifelse(N.tInf >= N.tT,
# exp(lgamma(N.tInf+size)-lgamma(N.tInf-N.tT+1) + N.tInf*log( (1-sumpd)*(mu/(mu+size)))),0)
#Compare the 2
## foo.a <- dbinom(N.tT, size=N.tInf, prob=sumpd)*dnbinom(N.tInf, mu=mu,size=size)
## foo.b <- ifelse(N.tInf >= N.tT, #& N.tInf <= size,
## exp(lgamma(N.tInf+size)-lgamma(N.tInf-N.tT+1) + N.tInf*log( (1-sumpd)*(mu/(mu+size)))),0)
## plot(foo.a/sum(foo.a))
## points(foo.b/sum(foo.b),col="red")
}
#Sample in binomial-negative-binomial hierarchy
rpost.bnb <- function(n=1, N.tT, sumpd, mu,size, N.tInf.max=1e4) {
p <- dpost.bnb.unorm(0:N.tInf.max,N.tT=N.tT,sumpd=sumpd, mu=mu,size=size)
#Set NA values to zero (why would they be NA?)
#if (is.na(sum(p))) { warning("rpost.bnb: sum is NA") ; browser(p)}
#Normalize the distribution - safe this for time reasons
#p <- p/sum(p)
#Sample
sample(0:N.tInf.max, size=n, replace=TRUE, prob=p)
}
#PMF for the predictive distribution in binomial-negative-binomial hierarchy.
#Returns entire vector for 0:N.tInf.max
dpost.bnb <- function(N.tT, sumpd, mu,size, N.tInf.max=1e4) {
p <- dpost.bnb.unorm(0:N.tInf.max,N.tT=N.tT,sumpd=sumpd, mu=mu,size=size)
#Set NA values to zero (why would they be NA?)
#if (is.na(sum(p))) { warning("rpost.bnb: sum is NA") ; browser(p)}
#Normalize the distribution - safe this for time reasons
return(p/sum(p))
}
######################################################################
# PMF of the beta-negatative binomial distribution
# See Teerapabolarn (2008)
#
# Parameters:
# k - where to evaluate. can be a vector.
#
# Returns:
# PMF.
######################################################################
dbnb <- function(k,n,alpha,beta) {
#Check if k's outside the support are requested.
neg <- k<0
k[neg] <- 0
#Calculate the density of the beta-negbin. See Teerapabolarn (2008)
num <- lgamma(n+alpha)+lgamma(k+beta)+lgamma(n+k)+lgamma(alpha+beta)
den <- lgamma(n+k+alpha+beta)+lgamma(n)+lgamma(k+1)+lgamma(alpha)+lgamma(beta)
res <- exp(num-den)
res[neg] <- 0
return( res)
}
######################################################################
# Convert discrete time hazard function on 0,...,Dmax to a probability
# mass function.
#
# Parameters:
# haz - vector with entries for (0,...,Dmax)
# Returns:
# vector with PMF on 0,...,Dmax.
######################################################################
haz2pmf <- function(haz) {
PMF <- 0*haz
for (i in 0:(length(haz)-1)) {
PMF[i+1] <- haz[i+1] * (1-sum(PMF[seq(i)]))
}
return(PMF)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.