moving: Moving Window Mean De-biasing

In jonasbhend/biascorrection: Bias Correct Ensemble Forecasts

Moving Window Mean De-biasing

Description

Computes moving window mean debiasing. A short-hand for multiplicative (moving_mul) mean debiasing is also provided. Furthermore, a mean debiasing approach in which the algorithm moves on in jumps rather than sequantially along lead times is provided via fastmoving and fastmoving_mul, respectively.

Usage

moving(
  fcst,
  obs,
  fcst.out = fcst,
  window = min(nrow(fcst), 31),
  jump = 1,
  multiplicative = FALSE,
  lower.bound = NULL,
  ...
)

moving_mul(..., multiplicative = TRUE)

fastmoving(fcst, window = min(nrow(fcst), 31), jump = min(window, 11), ...)

fastmoving_mul(
  fcst,
  window = min(nrow(fcst), 31),
  jump = min(window, 11),
  multiplicative = TRUE,
  ...
)

Arguments

fcst	n x m x k array of n lead times, m forecasts, of k ensemble members
obs	n x m matrix of veryfing observations
fcst.out	array of forecast values to which bias correction should be applied (defaults to fcst)
window	width of window of lead times to be used
jump	minimum number of lead times the moving window approach jumps
multiplicative	logical, is mean debiasing to be applied multiplicatively?
lower.bound	is used to truncate output if set (e.g. to zero for precipitation)
...	additional arguments for compatibility with other bias correction methods

Details

The mean debiasing is lead time dependent, parameter window is used to select the number of lead times around the target lead time to be used for bias correction. For the begining and end of the series, the lead-time interval is kept constant, so that to estimate the correction for the first lead time, the first window lead times are used.

Examples

## initialise forcast observation pairs
nens <- 51
signal <- outer(sin(seq(0,4,length=215)), sort(rnorm(30, sd=0.2)), '+') + 2
fcst <- list(raw=array(rgamma(length(signal)*nens, shape=2, scale=2), c(dim(signal), nens)) *
  c(signal) * (runif(length(signal)*nens) > 0.1))
obs <- array(rgamma(length(signal), shape=3, scale=1), dim(signal)) *
  signal * (runif(length(signal)) > 0.3)
fcst$moving <- biascorrection:::moving(fcst$raw[,1:20,], 
  obs[,1:20], fcst.out=fcst$raw, lower.bound=0)
fcst$moving_mul <- biascorrection:::moving_mul(fcst$raw[,1:20,], 
  obs[,1:20], fcst.out=fcst$raw, lower.bound=0)
oprob <- (seq(obs[,21:30]) - 1/3) / (length(obs[,21:30]) + 1/3)
oldpar <- par(no.readonly=TRUE)
par(mfrow=c(1,2))
plot(density(obs[,21:30], from=0, to=80, bw=1), type='n',
     main='Distribution in validation period')
for (i in 1:length(fcst)) lines(density(fcst[[i]][,21:30,1], from=0, to=80, bw=1), lwd=2, lty=i)
lines(density(obs[,21:30], from=0, to=80, bw=1), lwd=2, col=2)
legend('topright', c('Observations', 'No bias correction', names(fcst)[-1]), 
       lwd=2, col=c(2,rep(1, length(fcst))), lty=c(1,seq(fcst)), inset=0.05)
plot(quantile(obs[,21:30], type=8, oprob), 
  quantile(fcst[[1]][,21:30,], type=8, oprob),
  type='l', lwd=2, xlab='Observed quantiles',
  ylab='Forecast quantiles',
  main='Out-of-sample validation for moving')
abline(c(0,1), lwd=2, col=2)
for (i in 2:length(fcst)) lines(quantile(obs[,21:30], type=8, oprob),
  quantile(fcst[[i]][,21:30,], type=8, oprob), lwd=2, lty=i)
legend('topleft', c('No bias correction', names(fcst)[-1]), lwd=2, lty=seq(along=fcst), inset=0.05)
par(oldpar)

jonasbhend/biascorrection documentation built on Nov. 11, 2023, 1:16 a.m.