Adds a density region to an existing plot of a normally-distributed quantity with continuously-varying mean and standard deviation, such as a time series forecast. Automatically computes a reasonable set of ordinates to evaluate the density at, which span the whole forecast space.

1 2 | ```
## S3 method for class 'normal'
densregion(x, mean, sd, ny=20, ...)
``` |

`x` |
Suppose the continuously-varying quantity varies over a
space S. |

`mean` |
Vector of normal means at each point in |

`sd` |
Vector of standard deviations at each point in |

`ny` |
Minimum number of points to calculate the density at for
each |

`...` |
Further arguments passed to |

The plot is shaded by interpolating the value of the density
between grid points, using the algorithm described by Cleveland (1993)
as implemented in the `filled.contour`

function.

Christopher Jackson <chris.jackson@mrc-bsu.cam.ac.uk>

Jackson, C. H. (2008) Displaying uncertainty with
shading. *The American Statistician*, 62(4):340-347. http://www.mrc-bsu.cam.ac.uk/personal/chris/papers/denstrip.pdf

Cleveland, W. S. (1993) *Visualizing Data*. Hobart Press, Summit,
New Jersey.

`densregion`

, `densregion.survfit`

, `denstrip`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 | ```
## Time series forecasting
(fit <- arima(USAccDeaths, order = c(0,1,1),
seasonal = list(order=c(0,1,1))))
pred <- predict(fit, n.ahead = 36)
plot(USAccDeaths, xlim=c(1973, 1982), ylim=c(5000, 15000))
## Compute normal forecast densities automatically (slow)
densregion.normal(time(pred$pred), pred$pred, pred$se,
pointwise=TRUE, colmax="darkgreen")
lines(pred$pred, lty=2)
lines(pred$pred + qnorm(0.975)*pred$se, lty=3)
lines(pred$pred - qnorm(0.975)*pred$se, lty=3)
## Compute forecast densities by hand (more efficient)
nx <- length(pred$pred)
y <- seq(5000, 15000, by=100)
z <- matrix(nrow=nx, ncol=length(y))
for(i in 1:nx)
z[i,] <- dnorm(y, pred$pred[i], pred$se[i])
plot(USAccDeaths, xlim=c(1973, 1982), ylim=c(5000, 15000))
densregion(time(pred$pred), y, z, colmax="darkgreen", pointwise=TRUE)
lines(pred$pred, lty=2)
lines(pred$pred + qnorm(0.975)*pred$se, lty=3)
lines(pred$pred - qnorm(0.975)*pred$se, lty=3)
densregion(time(pred$pred), y+2000, z, colmax="darkblue", pointwise=TRUE)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.