densregion.normal {denstrip} | R Documentation |
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.
## 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.
Cleveland, W. S. (1993) Visualizing Data. Hobart Press, Summit, New Jersey.
densregion
, densregion.survfit
, denstrip
## 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)