FitAkdParameters(ens, obs)
p(y|ens) = 1 / K * sum {dnorm(y, z.i, s)} where s = (4/3/K)^0.4 * (s1 + s2 * a^2 * var(ens)) and z.i = r1 + r2 * mean(ens) + a * ens
The parameters r1, r2, a, s1, s2 are fitted by minimizing the continuously ranked probability score (CRPS) of the ensemble forecasts `ens` with respect to the observations `obs`. The optimization is carried out using the R function `optim(...)`.
The dressing variance is a linear function of the ensemble variance, and can thus become negative for certain parameter settings. In order to ensure nonnegativity, a logarithmic barrier is introduced during the optimization. That is, the method does not minimize the mean CRPS, but the function:
mean(CRPS(ens, obs)) * (1 + 0.01 * barrier) where barrier = max(0, -log( minimal dressing variance ))
Since the evaluation of the CRPS is numerically expensive, the optimization can take a long time. Speed can be increased by optimizing the parameters only for a part of the forecast instances.
# Example:
ens <- matrix(rnorm(500),100,5)
obs <- rnorm(100)
FitAkdParameters(ens, obs)
Run the code above in your browser using DataLab