diff --git a/DESCRIPTION b/DESCRIPTION
index 42f4b0c..5c8fb68 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -13,4 +13,4 @@ Imports:
     beanplot
 License: MIT
 LazyData: True
-RoxygenNote: 6.0.1
+RoxygenNote: 6.1.1
diff --git a/NAMESPACE b/NAMESPACE
index 2d8ddf6..f171123 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/OUProcessPlot.R b/R/OUProcessPlot.R
index 5b22ace..f0dfb39 100644
--- a/R/OUProcessPlot.R
+++ b/R/OUProcessPlot.R
@@ -1,6 +1,8 @@
-#' 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, 
@@ -8,12 +10,12 @@
 #' 
 #' 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)
   
-}
\ No newline at end of file
+}
+
+
+
+
+#' 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) 
+}
+
diff --git a/examples/OUProcessLkTests.R b/examples/OUProcessLkTests.R
index 26b401c..de9d08f 100644
--- a/examples/OUProcessLkTests.R
+++ b/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
\ No newline at end of file
diff --git a/man/plotCrossSection.Rd b/man/plotCrossSection.Rd
index 637274e..cea38ed 100644
--- a/man/plotCrossSection.Rd
+++ b/man/plotCrossSection.Rd
@@ -8,7 +8,8 @@ Logfile is a BEAST2 log file}
 plotCrossSection(par1, par2, lk, lk.upper = NULL, lk.lower = NULL,
   par1.truth = NA, par2.truth = NA, col = pal.dark(cblue),
   fill = pal.dark(cblue, 0.5), xlims = NULL, ylims = NULL,
-  par1.label = NA, par2.label = NA, lk.label = NA, sections = -1, ...)
+  par1.label = NA, par2.label = NA, lk.label = NA, sections = -1,
+  ...)
 }
 \description{
 Plot cross section of the likelihood surface of par1, for given values of par2
diff --git a/man/plotHPDs.Rd b/man/plotHPDs.Rd
index c2944e6..837d34a 100644
--- a/man/plotHPDs.Rd
+++ b/man/plotHPDs.Rd
@@ -7,12 +7,12 @@
 plotHPDs(data, type = "beans", hpdlines = c(0.05, 0.5),
   linewidths = c(0.75, 1), maxwidth = 0.5, hpdonly = TRUE,
   showlimits = FALSE, plotmedian = TRUE, plothalf = FALSE,
-  col = pal.dark(cblack), fill = pal.dark(cgray, 0.25), col.axis = NA,
-  lwd = 1, add = FALSE, new = TRUE, xaxis = TRUE, yaxis = TRUE,
-  xlab = "", ylab = "", xline = 1, yline = 1, xticks = NULL,
-  yticks = NULL, xticklabels = NULL, ygrid = FALSE, xgrid = FALSE,
-  gridcol = pal.light(cgray), ylims = NULL, cex.label = 1.4,
-  cex.axis = 1, axispadding = 0, side = 2, ...)
+  col = pal.dark(cblack), fill = pal.dark(cgray, 0.25),
+  col.axis = NA, lwd = 1, add = FALSE, new = TRUE, xaxis = TRUE,
+  yaxis = TRUE, xlab = "", ylab = "", xline = 1, yline = 1,
+  xticks = NULL, yticks = NULL, xticklabels = NULL, ygrid = FALSE,
+  xgrid = FALSE, gridcol = pal.light(cgray), ylims = NULL,
+  cex.label = 1.4, cex.axis = 1, axispadding = 0, side = 2, ...)
 }
 \arguments{
 \item{data}{List with each element being a vector that represents a posterior sample
diff --git a/man/plotOUDensity.Rd b/man/plotOUDensity.Rd
new file mode 100644
index 0000000..3f33156
--- /dev/null
+++ b/man/plotOUDensity.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/OUProcessPlot.R
+\name{plotOUDensity}
+\alias{plotOUDensity}
+\title{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.}
+\usage{
+plotOUDensity(xrange, dt, x0, mu, sigma, nu, col = pal.dark(cblue),
+  col.alpha = 0.25, new = TRUE, plotMedian = FALSE, plotX0 = TRUE,
+  plotGrid = TRUE, ...)
+}
+\description{
+Density is only plotted between the 0.025 and 0.975 quantiles.
+}
diff --git a/man/plotOUDensityEmpirical.Rd b/man/plotOUDensityEmpirical.Rd
new file mode 100644
index 0000000..7abdbb6
--- /dev/null
+++ b/man/plotOUDensityEmpirical.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/OUProcessPlot.R
+\name{plotOUDensityEmpirical}
+\alias{plotOUDensityEmpirical}
+\title{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).}
+\usage{
+plotOUDensityEmpirical(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, ...)
+}
+\description{
+Density is only plotted in the 95% HPD interval, at the points in xrange.
+}
diff --git a/man/plotOUProcessFlexibility.Rd b/man/plotOUProcessFlexibility.Rd
new file mode 100644
index 0000000..189b5db
--- /dev/null
+++ b/man/plotOUProcessFlexibility.Rd
@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/OUProcessPlot.R
+\name{plotOUProcessFlexibility}
+\alias{plotOUProcessFlexibility}
+\title{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}
+\usage{
+plotOUProcessFlexibility(xlim, dt, x0, mu, sigma, nu, nSteps = 100)
+}
+\description{
+Answers the questions: 
+- How big is the 95% CI after one time interval?
+}
diff --git a/man/plotOUProcessFlexibilityEmpirical.Rd b/man/plotOUProcessFlexibilityEmpirical.Rd
new file mode 100644
index 0000000..094eeea
--- /dev/null
+++ b/man/plotOUProcessFlexibilityEmpirical.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/OUProcessPlot.R
+\name{plotOUProcessFlexibilityEmpirical}
+\alias{plotOUProcessFlexibilityEmpirical}
+\title{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}
+\usage{
+plotOUProcessFlexibilityEmpirical(xlim, dt, x0, mu_p, sigma_p, nu_p,
+  nSteps = 100)
+}
+\description{
+Answers the questions: 
+- How big is the 95% CI after one time interval?
+}
diff --git a/man/plotOUProcessHPD.Rd b/man/plotOUProcessHPD.Rd
index 1a9c5eb..20d5c46 100644
--- a/man/plotOUProcessHPD.Rd
+++ b/man/plotOUProcessHPD.Rd
@@ -2,15 +2,15 @@
 % Please edit documentation in R/OUProcessPlot.R
 \name{plotOUProcessHPD}
 \alias{plotOUProcessHPD}
-\title{Functions for plotting OU processes}
+\title{Plot expected 95% quantiles (0.025 and 0.975) and median of OU-process, 
+conditioning on values of x0, mu, sigma, nu}
 \usage{
 plotOUProcessHPD(x0, t, mu, sigma, nu, nGrad = 101, plotGrad = TRUE,
   ylim = NULL, ...)
 }
 \description{
-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)
 }
 \details{
-Can also plot a gradient of quantiles (by default plots percentiles - n=101)
+TODO: Name is misleading, quantiles are NOT HPDs.
 }
diff --git a/man/plotOUProcessHPDEmpirical.Rd b/man/plotOUProcessHPDEmpirical.Rd
index e5db8b1..3b1c6a0 100644
--- a/man/plotOUProcessHPDEmpirical.Rd
+++ b/man/plotOUProcessHPDEmpirical.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/OUProcessPlot.R
 \name{plotOUProcessHPDEmpirical}
 \alias{plotOUProcessHPDEmpirical}
-\title{Plot empirical 95% quantiles (0.025 and 0.975) and median of OU-process, 
+\title{Plot empirical 95% HPD interval and median of OU-process, 
 from simulating nTraj replicates, which can also be plotted}
 \usage{
 plotOUProcessHPDEmpirical(x0_prior, t, mu_prior, sigma_prior, nu_prior,
diff --git a/man/plotPrior.Rd b/man/plotPrior.Rd
index 7ead521..6fe2e41 100644
--- a/man/plotPrior.Rd
+++ b/man/plotPrior.Rd
@@ -4,11 +4,11 @@
 \alias{plotPrior}
 \title{Plot a parametric distribution where the distribution is already available in R}
 \usage{
-plotPrior(priorfun = "norm", prior_args = list(), col = pal.dark(cblue),
-  fill = pal.dark(cblue, 0.25), ylab = "Density", xlab = "x",
-  plotquantile = TRUE, grid = TRUE, positive = TRUE, invert = FALSE,
-  meaninrealspace = FALSE, scaling = 1, add = FALSE, xlims = NULL,
-  ylims = NULL, ...)
+plotPrior(priorfun = "norm", prior_args = list(),
+  col = pal.dark(cblue), fill = pal.dark(cblue, 0.25),
+  ylab = "Density", xlab = "x", plotquantile = TRUE, grid = TRUE,
+  positive = TRUE, invert = FALSE, meaninrealspace = FALSE,
+  scaling = 1, add = FALSE, xlims = NULL, ylims = NULL, ...)
 }
 \arguments{
 \item{priorfun}{Family of distributions e.g. norm or beta}
diff --git a/man/plotSkyline.Rd b/man/plotSkyline.Rd
index 86c2430..b33886a 100644
--- a/man/plotSkyline.Rd
+++ b/man/plotSkyline.Rd
@@ -6,7 +6,8 @@
 \usage{
 plotSkyline(times, skyline_mat, type = "smooth", traces = 1000,
   col = pal.dark(cblack), fill = pal.dark(cgray, 0.25), lwd = 1,
-  lty = 1, new = TRUE, add = FALSE, xlims = NULL, ylims = NULL, ...)
+  lty = 1, new = TRUE, add = FALSE, xlims = NULL, ylims = NULL,
+  ...)
 }
 \arguments{
 \item{times}{The time points to draw the skyline at. length(times) should be equal to ncol(skyline_mat) + 1.
diff --git a/man/plotSkylinePretty.Rd b/man/plotSkylinePretty.Rd
index ca31620..912f1de 100644
--- a/man/plotSkylinePretty.Rd
+++ b/man/plotSkylinePretty.Rd
@@ -5,12 +5,12 @@
 \title{Plot a pretty skyline. Contains default options for drawing axes to the plot}
 \usage{
 plotSkylinePretty(times, skyline_mat, type = "smooth", traces = 1000,
-  col = pal.dark(cblack), fill = pal.dark(cgray, 0.25), col.axis = NA,
-  lwd = 1, lty = 1, xaxis = TRUE, yaxis = TRUE, xlab = "",
-  ylab = "", xline = 1, yline = 1, xticks = NULL, yticks = NULL,
-  xticklabels = NULL, xlims = NULL, ylims = NULL, xgrid = FALSE,
-  ygrid = FALSE, gridcol = pal.light(cgray), cex.label = 1.4,
-  cex.axis = 1, axispadding = 0, side = 2, ...)
+  col = pal.dark(cblack), fill = pal.dark(cgray, 0.25),
+  col.axis = NA, lwd = 1, lty = 1, xaxis = TRUE, yaxis = TRUE,
+  xlab = "", ylab = "", xline = 1, yline = 1, xticks = NULL,
+  yticks = NULL, xticklabels = NULL, xlims = NULL, ylims = NULL,
+  xgrid = FALSE, ygrid = FALSE, gridcol = pal.light(cgray),
+  cex.label = 1.4, cex.axis = 1, axispadding = 0, side = 2, ...)
 }
 \arguments{
 \item{xline}{Line to draw x-axis label on}