R/DressCrps.R

In SpecsVerification: Forecast Verification Routines for Ensemble Forecasts of Weather and Climate

#' Calculate the Continuous Ranked Probability Score (CRPS) for a mixture of Normal distributions, for example generated by ensemble dressing
#'
#' @param dressed.ens a list with elements `ens`, a N*R matrix representing N time instances of kernel centers, and `ker.wd`, a N*R matrix with corresponding kernel standard deviations. See function `DressEnsemble`
#' @param obs a numeric vector of length N with real-valued observations
#' @return numeric vector of length N with the CRPS values
#' @examples
#' data(eurotempforecast)
#' dressed.ens <- DressEnsemble(ens)
#' mean(DressCrps(dressed.ens, obs))
#' @seealso EnsCrps, ScoreDiff, SkillScore
#' @references Grimit et al (2006): The continuous ranked probability score for circular variables and its application to mesoscale forecast ensemble verification. Q.J.R. Meteorol. Soc. \doi{10.1256/qj.05.235}
#' @export
DressCrps <- function(dressed.ens, obs) {

  stopifnot(all(dim(dressed.ens[["ens"]]) == dim(dressed.ens[["ker.wd"]])))

  N <- nrow(dressed.ens[["ens"]])
  K.vec <- rowSums(!is.na(dressed.ens[["ens"]]))

  crps <- sapply(
            seq_len(nrow(dressed.ens[['ens']])), 
            function(ii) {
              dresscrps_cpp(dressed.ens[['ens']][ii,], 
                            dressed.ens[['ker.wd']][ii,], 
                            obs[ii])
            })

  return(crps)
}