Nothing
#' exDQLM Diagnostics
#'
#' The function computes the following for the model(s) provided: the posterior predictive loss criterion based off the check loss, the one-step-ahead distribution sequence and its KL divergence from normality. The function also plots the following: the qq-plot and ACF plot corresponding to the one-step-ahead distribution sequence, and a time series plot of the MAP standard forecast errors.
#'
#' @inheritParams exdqlmPlot
#' @param m2 An optional additional object of class "`exdqlm`" to compare with `m1`.
#' @param plot If `TRUE`, the following will be plotted for `m1` and `m2` (if provided): a qq-plot and ACF plot of the MAP one-step-ahead distribution sequence, and a time series plot of the standardized forecast errors.
#' @param cols Color(s) used to plot diagnostics.
#' @param ref Reference sample of size `length(y)` from a standard normal distribution used to compute the KL divergence.
#'
#' @return A list containing the following is returned:
#' \itemize{
#' \item `m1.uts` - The one-step-ahead distribution sequence of `m1`.
#' \item `m1.KL` - The KL divergence of `m1.uts` and a standard normal.
#' \item `m1.pplc` - The posterior predictive loss criterion of `m1` based off the check loss function.
#' \item `m1.qq` - The ordered pairs of the qq-plot comparing `m1.uts` with a standard normal distribution.
#' \item `m1.acf` - The autocorrelations of `m1.uts` by lag.
#' }
#' If `m2` is provided, analogous results for `m2` are also included in the list.
#' @export
#'
#' @examples
#' \donttest{
#' y = scIVTmag[1:100]
#' model = polytrendMod(1,mean(y),10)
#' M0 = exdqlmISVB(y,p0=0.85,model,df=c(0.95),dim.df = c(1),
#' gam.init=-3.5,sig.init=15)
#' check.out = exdqlmChecks(y,M0,plot=FALSE)
#' check.out$m1.KL
#' check.out$m1
#' }
#'
exdqlmChecks <- function(y,m1,m2=NULL,plot=TRUE,cols=c("grey","grey"),ref=NULL){
# check inputs
check_ts(y)
TT = length(y)
if(!is.exdqlm(m1)){
stop("m1 must be an output from 'exdqlmISVB()' or 'exdqlmMCMC()'")
}
if(dim(m1$samp.theta)[2] != TT){
stop("length of dynamic quantile in m1 does not match data")
}
cols = c(matrix(cols,2,1))
### m1
# m1 seq.uts
m1.uts = stats::pnorm(m1$map.standard.forecast.errors)
# m1 KL divergence
TT = length(m1$map.standard.forecast.errors)
if(is.null(ref)){
ref = stats::rnorm(TT)
}else{
ref = c(ref)
if(length(ref)!=TT){
stop("ref should be a sample of size 'length(y)' from a standard normal distribution")
}
}
m1.KL = mean(FNN::KL.divergence(ref,m1$map.standard.forecast.errors))
# m1 pplc
m1.loss = matrix(NA,TT,dim(m1$samp.post.pred)[2])
for(t in 1:TT){m1.loss[t,] = CheckLossFn(m1$p0,y[t]-m1$samp.post.pred[t,])}
m1.pplc = sum(rowMeans(m1.loss))
# m1 qqplot
m1.qq = stats::qqnorm(m1$map.standard.forecast.errors,plot=FALSE)
# m1 acf
m1.acf = stats::acf(m1.uts,plot=FALSE)
#
retlist = list(m1.uts=m1.uts,m1.KL=m1.KL,m1.pplc=m1.pplc,m1.qq=m1.qq,m1.acf=m1.acf)
### m2
if(!is.null(m2)){
# check inputs
if(!is.exdqlm(m2)){
stop("m2 must be an output from 'exdqlmISVB()' or 'exdqlmMCMC()'")
}
if(dim(m1$samp.theta)[2] != TT){
stop("length of dynamic quantile in m2 does not match data")
}
if(TT != length(m2$map.standard.forecast.errors)){
stop("length of m1 quantile does not match length of m2 quantile")
}
if(m1$p0 != m2$p0){
stop("quantiles estimated in m1 and m2 do not match")
}
# m2 seq.uts
m2.uts = stats::pnorm(m2$map.standard.forecast.errors)
retlist[["m2.uts"]] = m2.uts
# m2 KL divergence
retlist[["m2.KL"]] = mean(FNN::KL.divergence(ref,m2$map.standard.forecast.errors))
# m2 pplc
m2.loss = matrix(NA,TT,dim(m2$samp.post.pred)[2])
for(t in 1:TT){m2.loss[t,] = CheckLossFn(m2$p0,y[t]-m2$samp.post.pred[t,])}
retlist[["m2.pplc"]] = sum(rowMeans(m2.loss))
# m2 qqplot
retlist[["m2.qq"]] = stats::qqnorm(m2$map.standard.forecast.errors,plot=FALSE)
# m2 acf
retlist[["m2.acf"]] = stats::acf(m2.uts,plot=FALSE)
}
if(plot){
if(is.null(m2)){
qq.x.range = range(m1.qq$x)
qq.y.range = range(m1.qq$y)
acf.y.range = range(m1.acf$acf)
fe.y.range = range(m1$map.standard.forecast.errors)
}else{
qq.x.range = range(c(m1.qq$x,retlist$m2.qq$x))
qq.y.range = range(c(m1.qq$y,retlist$m2.qq$y))
acf.y.range = range(c(m1.acf$acf,retlist$m2.acf$acf))
fe.y.range = range(c(m1$map.standard.forecast.errors,m2$map.standard.forecast.errors))
}
# m1 qqplot
plot(m1.qq,main="",col=cols[1],pch=20,xlab="Theoretical Quantiles",ylab="M1 Sample Quantiles",xlim=qq.x.range,ylim=qq.y.range)
graphics::abline(a=0,b=1)
# m1 acf
plot(m1.acf,ylab="M1 ACF",col=cols[1],main="",ylim=acf.y.range)
# m1 forecast errors
ts.xy = grDevices::xy.coords(y)
graphics::plot(ts.xy$x,m1$map.standard.forecast.errors,ylab="M1 standard forecast errors",xlab="time",col=cols[1],pch=20,type="l",ylim=fe.y.range)
graphics::abline(h=0,lty=2)
### m2
if(!is.null(m2)){
# m2 qqplot
plot(retlist[["m2.qq"]],main="",col=cols[2],pch=20,xlab="Theoretical Quantiles",ylab="M2 Sample Quantiles",xlim=qq.x.range,ylim=qq.y.range)
graphics::abline(a=0,b=1)
# m2 acf
plot(retlist[["m2.acf"]],ylab="M2 ACF",col=cols[2],main="",ylim=acf.y.range)
# m2 forecast errors
graphics::plot(ts.xy$x,m2$map.standard.forecast.errors,ylab="M2 standard forecast errors", xlab="time",col=cols[2],pch=20,type="l",ylim=fe.y.range)
graphics::abline(h=0,lty=2)
}
}
# return model checks
return(invisible(retlist))
}
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.