## Spread skill relationship ##
## load and prepare data
data(RainIbk)
## mean and standard deviation of square root transformed ensemble forecasts
RainIbk$sqrtensmean <-
apply(sqrt(RainIbk[,grep('^rainfc',names(RainIbk))]), 1, mean)
RainIbk$sqrtenssd <-
apply(sqrt(RainIbk[,grep('^rainfc',names(RainIbk))]), 1, sd)
## quintiles of sqrtenssd
sdcat <- cut(RainIbk$sqrtenssd, c(-Inf, quantile(RainIbk$sqrtenssd,
seq(0.2,0.8,0.2)), Inf), labels = c(1:5))
## mean forecast errors for each quintile
m <- NULL
for(i in levels(sdcat)) {
m <- c(m, mean((sqrt(RainIbk$rain)[sdcat == i] -
RainIbk$sqrtensmean[sdcat == i])^2, na.rm = TRUE))
}
## plot
boxplot((sqrt(rain) - sqrtensmean)^2~sdcat, RainIbk,
xlab = "Quintile of ensemble standard deviation",
ylab = "mean squared error", main = "Spread skill relationship")
Run the code above in your browser using DataLab