fanplot (version 3.1)

fan0: Fan Plot of Distributions Percentiles Over Time.

Description

This function is no longer recommended. Was intended for use with original pn function. Replaced by fan, which avoids the requirement of a pn object and offers a greater range of styles. This function is called within fan, is psims is a pn.ts object. Plots percentiles of sequential distributions provided by pn object. Typically called after the plot area has been set.

Usage

fan0(psims, fan.col = heat.colors(floor(dim(psims)[2]/2)), 
      ln.col = fan.col[1], ln = seq(10,90,10), txt = ln, ...)

Arguments

psims

pn object created by the pn function.

fan.col

Palette of colours used in the fan plot

ln.col

Line colour to be imposed on top of the fan

ln

Vector of number to plot contour lines on top the fan. Must correspond to calculated percentiles in psims. By default takes every decile contained within the pn object. If the

txt

Vector of number to plot text at the end of corresponding percentiles. Must be in calculated percentiles of psims.

...

Additional arguments passed to lines.

Value

• See Details

Details

Returns a plot of the percentiles given in the psims argument. Additional lines and text are added to illustrate major contours on the probability distribution. By default these run on every decile. Text is automatically added where lines are drawn. Hence, if text is only required where there are contour lines, there is no need to specify anything in the txt argument. Text can also be suppressed and added later using fan.txt. Lines can be suppressed by adding type = NULL. Colours are by default taken from the heat.colors palette. Alternatives can be specified using fan.col. The joining of a forecast fan to data is controlled be the anchor argument in the pn function.

See Also

pn, fan.txt, fan.fill

Examples

# calculate percentiles across time
th.pn <- pn(sims = th.mcmc)

##
##Default
##
# empty plot
plot(NULL, type = "n", xlim = c(1, 945),  ylim = range(th.pn), ylab = "Theta")

# add fan
fan0(th.pn)

##
##Text at both ends and lines also at 1st and 99th percentiles
##
# empty plot
plot(NULL, type = "l", xlim = c(1, 945), xlab = "Time", xaxt = "n",
     ylim = range(th.pn), ylab = "Theta")

# add fan
fan0(th.pn, txt = NULL, ln = c(1,10,30,50,70,90,99))

# add text labels for percentiles
fan.txt(th.pn, pn.r = c(1,10,50,90,99))
fan.txt(th.pn, pn.l = c(1,30,70,99))


##
##Different colours
##
# calculate percentiles across time
sigma.pn <- pn(sims = sqrt(exp(th.mcmc)))

# define colours
pal <- colorRampPalette(c("royalblue", "grey", "white"))
fancol <- pal(50)

# empty plot
plot(NULL, type = "l", xlim = c(1, 945), xlab = "Time", xaxt = "n",
     ylim = range(sigma.pn), ylab = "Standard Deviation")
     
# add fan
fan0(sigma.pn, fan.col = fancol, ln = c(1, 10, 50, 90, 99))

##
##Forecast fans 1) default 2) coarser colours, different contour line colour and anchoring fan
##
# create BUGS model
library("tsbugs")
r <- ts(ew[2:167]/ew[1:166]-1, start=1841)
y <- diff(r) 
pop.bug <- sv.bugs(y, k=25, sim=TRUE,
                  sv.mean.prior2 = "dgamma(0.000001,0.000001)", 
                  sv.ar.prior2 = "dunif(-0.999, 0.999)")

# run in OpenBUGS
library("R2OpenBUGS")
writeLines(pop.bug$bug, "pop.txt")
pop <- bugs(data = pop.bug$data,
            inits = list(list(psi0.star=exp(12), psi1=0.5, itau2=0.5)),
            param = c("psi0", "psi1", "tau", "y.new", "y.sim"),
            model = "pop.txt", 
            n.iter = 11000, n.burnin = 1000, n.chains = 1)

# derive posterior predictive distributions of the population growth rate and population total
ynew.mcmc <- pop$sims.list$y.new
rnew.mcmc <- apply(ynew.mcmc, 1, diffinv, xi = tail(r,1))
rnew.mcmc <- t(rnew.mcmc[-1,])

pnew.mcmc <- apply(1+rnew.mcmc, 1, cumprod) * tail(ew,1)
pnew.mcmc <- t(pnew.mcmc)

# calculate percentiles across time
r0 <- tsp(r)[2]
rnew.pn <- pn(sims = rnew.mcmc, start = r0 + 1)
pnew.pn <- pn(sims = pnew.mcmc/1e+06, start = r0 + 1, anchor = tail(ew/1e+06, 1),) 

# plot ew
par(mfrow = c(1 ,2))
plot(ew/1e+06, ylim = c(40, 80), xlim = c(1940, 2040), ylab = "Population (m)", lwd = 2)
     
# add fan    
fan0(pnew.pn)

# plot ew
plot(ew/1e+06, ylim = c(40, 80), xlim = c(1940, 2040), ylab = "Population (m)", lwd = 2)
     
# add fan
fan0(pnew.pn2, ln.col = "black")