Add some OU-process utilities

DESCRIPTION

@@ -13,4 +13,4 @@ Imports:
     beanplot
 License: MIT
 LazyData: True
-RoxygenNote: 6.0.1
+RoxygenNote: 6.1.1

NAMESPACE

@@ -37,6 +37,10 @@ export(plotCrossSection)
 export(plotHPDs)
 export(plotJellyBeanHPDs)
 export(plotJellyBeans)
+export(plotOUDensity)
+export(plotOUDensityEmpirical)
+export(plotOUProcessFlexibility)
+export(plotOUProcessFlexibilityEmpirical)
 export(plotOUProcessHPD)
 export(plotOUProcessHPDEmpirical)
 export(plotOUProcessPriors)

R/OUProcessPlot.R

@@ -1,19 +1,21 @@
-#' Functions for plotting OU processes
-#'
-
+#######################################
+# Functions for plotting OU processes #
+#######################################
+# Louis du Plessis, 2018              #
+#######################################
 
 
 #' Plot expected 95% quantiles (0.025 and 0.975) and median of OU-process, 
 #' conditioning on values of x0, mu, sigma, nu
 #' 
 #' Can also plot a gradient of quantiles (by default plots percentiles - n=101)
 #' 
+#' TODO: Name is misleading, quantiles are NOT HPDs.
 #' 
 #' @export
 plotOUProcessHPD <- function(x0, t, mu, sigma, nu, 
                              nGrad=101, plotGrad=TRUE, ylim=NULL, ...) {
 
-
   # Get theoretical 95% HPD and median
   limits <- quantilesOU(c(0.025, 0.5, 0.975), x0, t, mu, sigma, nu)
 
@@ -37,7 +39,7 @@ plotOUProcessHPD <- function(x0, t, mu, sigma, nu,
 }
 
 
-#' Plot empirical 95% quantiles (0.025 and 0.975) and median of OU-process, 
+#' Plot empirical 95% HPD interval and median of OU-process, 
 #' from simulating nTraj replicates, which can also be plotted
 #' 
 #' Instead of specifying constant parameter values, prior functions can also
@@ -49,7 +51,7 @@ plotOUProcessHPDEmpirical <- function(x0_prior, t, mu_prior, sigma_prior, nu_pri
 
   # Setup
   X <-  matrix(0, nrow=nTraj, ncol=length(t))
-  plot(1,type="n",xlim=range(t), ylim=ylim,axes=TRUE,xlab=NA,ylab=NA, ...)    
+  plot(1,type="n",xlim=range(t), ylim=ylim,axes=TRUE, ...)    
 
   # Do simulations and plot trajectories
   for (i in 1:nTraj) {
@@ -72,4 +74,177 @@ plotOUProcessHPDEmpirical <- function(x0_prior, t, mu_prior, sigma_prior, nu_pri
 
   return(X)
 
-}
+}
+
+
+
+
+#' Plot the density of an OU-process with parameters x0, mu, sigma, nu after time-interval
+#' dt, at the range of values given by xrange.
+#' 
+#' Density is only plotted between the 0.025 and 0.975 quantiles.
+#' 
+#' @export
+plotOUDensity <- function(xrange, dt, x0, mu, sigma, nu, col=pal.dark(cblue), col.alpha=0.25, 
+                          new=TRUE, plotMedian=FALSE, plotX0=TRUE, plotGrid=TRUE, ...) { 
+
+  OUdensity <- densityOU(xrange,dt,x0,mu,sigma,nu)
+  OUlimits  <- quantilesOU(c(0.025, 0.5, 0.975), x0, dt, mu, sigma, nu)
+  CIlimits  <- findInterval(OUlimits, xrange)
+
+  if (new) {
+    plot(c(xrange[1],xrange), c(0,OUdensity), type='n', ...)
+    if (plotGrid) {
+        grid(col='black')
+    }
+  }
+
+  lines(xrange[CIlimits[1]:CIlimits[3]], OUdensity[CIlimits[1]:CIlimits[3]], lwd=2, col=col, xpd=TRUE)
+  polygon(c(xrange[max(1,CIlimits[1])],xrange[CIlimits[1]:CIlimits[3]],xrange[CIlimits[3]]), c(0,OUdensity[CIlimits[1]:CIlimits[3]],0), col=scales::alpha(col,col.alpha), border=NA)
+
+  if (plotX0) {
+    points(x0,0,pch=16,cex=2,col=col,xpd=TRUE)
+  }
+
+  if (plotMedian) {
+    abline(v=OUlimits[2], col=col, lwd=2, lty=2)
+  }
+
+}
+
+
+#' Plot the density of ntraj OU-processes simulated to time dt, using priors
+#' x0_prior, mu_prior, sigma_prior and nu_prior on the OU-process parameters
+#' (fixed values are also supported).
+#' 
+#' Density is only plotted in the 95% HPD interval, at the points in xrange.
+#'
+#' @export
+plotOUDensityEmpirical <- function(xrange, dt, x0_prior, mu_prior, sigma_prior, nu_prior, nTraj=1000, nSteps=100, col=pal.dark(cblue), col.alpha=0.25, 
+                                   bw='sj', new=TRUE, plotMedian=FALSE, plotX0=TRUE, plotGrid=TRUE, ...) { 
+
+
+  # Setup
+  t <- seq(0,dt,length.out=nSteps)
+  X <- matrix(0, nrow=nTraj, ncol=length(t))
+
+  # Do simulations and plot trajectories
+  for (i in 1:nTraj) {
+    x0    <- eval(x0_prior)
+    mu    <- eval(mu_prior)
+    sigma <- eval(sigma_prior)
+    nu    <- eval(nu_prior)
+
+    # Simulate trajectory
+    X[i,] <- simulateOU(x0, t, mu, sigma, nu)
+  }
+  # Empirical 95% HPD and median
+  OUlimits  <- getHPD(X[,ncol(X)])
+  OUdensity <- density(X[,ncol(X)], bw=bw, from=OUlimits[1], to=OUlimits[3])
+
+  if (new) {
+    plot(c(xrange[1],OUdensity$x,xrange[length(xrange)]), c(0,OUdensity$y,0), type='n', ...)
+    if (plotGrid) {
+        grid(col='black')
+    }
+  }
+
+  lines(OUdensity$x, OUdensity$y, lwd=2, col=col)
+  polygon(c(OUdensity$x[1], OUdensity$x, OUdensity$x[length(OUdensity$x)]), c(0, OUdensity$y, 0), col=scales::alpha(col,col.alpha), border=NA)
+
+  if (plotX0) {
+    points(x0,0,pch=16,cex=2,col=col,xpd=TRUE)
+  }
+
+  if (plotMedian) {
+    abline(v=OUlimits[2], col=col, lwd=2, lty=2)
+  }
+
+}
+
+
+#' Plot a summary of the flexibility of an OU-process with parameters mu, sigma and nu, starting 
+#' from different x0 values, over a time period of dt
+#' 
+#' Answers the questions: 
+#' - How big is the 95% CI after one time interval? 
+#  - How much can the mean revert after one time interval?
+#' 
+#' @export
+plotOUProcessFlexibility <- function(xlim, dt, x0, mu, sigma, nu, nSteps=100) { 
+
+  t    <- seq(0,dt, by=dt/nSteps)
+  x    <- seq(xlim[1],xlim[2], length.out=1000)
+  n    <- length(x0)
+  cols <- RColorBrewer::brewer.pal(max(3,n),"Set2")
+
+  layout(matrix(c(1:n,rep(n+1,n)),byrow=TRUE, nrow=2))
+  for (i in 1:n) {
+    plotOUProcessHPD(x0[i], t, mu, sigma, nu, ylim=range(x), xlab="t",ylab="x",main=paste("x0 =",x0[i]))      
+  }
+
+  for (i in 1:n) {
+    plotOUDensity(x, dt, x0[i], mu, sigma, nu, col=cols[i], new=ifelse(i==1,TRUE,FALSE),xaxs='i', yaxs='i', xlab='x', ylab="", bty='n')
+  } 
+  abline(v=mu,col='red',lty=3, lwd=2)
+  legend('topright',col=cols, pch=16, bty='n', legend=paste("x0 =",x0), cex=1.25)
+}
+
+
+
+
+#' Plot a summary of the flexibility of an OU-process with parameters sampled from priors 
+#' mu_prior, sigma_prior and nu_prior, starting from different x0 values, over a time period 
+#' of dt 
+#' 
+#' Answers the questions: 
+#' - How big is the 95% CI after one time interval? 
+#  - How much can the mean revert after one time interval?
+#' 
+#' @export
+plotOUProcessFlexibilityEmpirical <- function(xlim, dt, x0, mu_p, sigma_p, nu_p, nSteps=100) { 
+
+  t    <- seq(0,dt, by=dt/nSteps)
+  x    <- seq(xlim[1],xlim[2], length.out=1000)
+  n    <- length(x0)
+  cols <- RColorBrewer::brewer.pal(max(3,n),"Set2")
+
+  layout(matrix(c(1:n,rep(n+1,n)),byrow=TRUE, nrow=2))
+  for (i in 1:n) {
+    plotOUProcessHPDEmpirical(x0[i], t, mu_p, sigma_p, nu_p, nTraj=1000, plotTraj=TRUE, ylim=range(x), xlab="t",ylab="x", main=paste("x0 =",x0[i]))
+  }
+
+  for (i in 1:n) {
+    plotOUDensityEmpirical(x, dt, x0[i], mu_p, sigma_p, nu_p, col=cols[i], new=ifelse(i==1,TRUE,FALSE),xaxs='i', yaxs='i', xlab='x', ylab="", bty='n')
+  } 
+
+  mu <- c()
+  for (i in 1:5000) {
+    mu <- c(mu, eval(mu_p))
+  }
+  #print(length(mu))
+  abline(v=median(mu),col='red',lty=3, lwd=2)
+  legend('topright',col=cols, pch=16, bty='n', legend=paste("x0 =",x0), cex=1.25)
+}
+
+
+
+###############################################################################
+
+testFlexibility <- function() {
+
+  x0   <- c(1,1.5,2,5)
+  dt   <- 0.5
+  xlim <- c(0,8)
+
+  nu    <- 1
+  mu    <- 1
+  sigma <- 1
+  plotOUProcessFlexibility(xlim, dt, x0, mu, sigma, nu) 
+
+  mu_p    <- getPrior(rlnorm, 1, meanlog=0, sdlog=1.25)
+  sigma_p <- getPrior(rlnorm, 1, meanlog=0, sdlog=0.25)
+  nu_p    <- getPrior(rexp, 1, rate=1)
+  plotOUProcessFlexibilityEmpirical(xlim, dt, x0, mu_p, sigma_p, nu_p) 
+}
+

examples/OUProcessLkTests.R

@@ -74,6 +74,13 @@ nu    <- 2
 l     <- logLikOU(x,t,mu,sigma,nu, removeconstant = TRUE)
 print(paste0("{",paste(x,collapse=","),"}  {",paste(t,collapse=","),"}  ",l))
 
+x     <- c(10,9,8,7,6,5,5,5,5,6,5)
+t     <- seq(from=0, to=10, by=1)
+x0    <- 10
+mu    <- 5
+sigma <- 1
+nu    <- 2
+l     <- logLikOU(x,t,mu,sigma,nu, removeconstant = TRUE)
+print(paste0("{",paste(x,collapse=","),"}  {",paste(t,collapse=","),"}  ",l))
 
 
-# 4) Log-likelihood of a trajectory with priors