densregion.normal {denstrip} | R Documentation |
Density regions based on normal distributions
Description
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.
Usage
## S3 method for class 'normal'
densregion(x, mean, sd, ny=20, ...)
Arguments
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 |
Details
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.
Author(s)
Christopher Jackson <chris.jackson@mrc-bsu.cam.ac.uk>
References
Jackson, C. H. (2008) Displaying uncertainty with shading. The American Statistician, 62(4):340-347.
Cleveland, W. S. (1993) Visualizing Data. Hobart Press, Summit, New Jersey.
See Also
densregion
, densregion.survfit
, denstrip
Examples
## 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)
## Not run:
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)
## End(Not run)
## 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)