ordinal brier scoring

ordinal_brier.R

###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
### calculates Brier Score for a single forecast on an Ordered IFP
### calculation based on definition in Appendix B of ACE_TE_Plan_OY1_v3.
### procedure reduces penalty to probability assigned to bins adjacent
### to the realized outcome
###~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ScoreOrdered <- function(fcast.values, 
                         correct.values,
                         scoring.fn){
  # Step 1 – Take the original categories (A-B-C-D) and break them up 
  # into a set of binary categories (A-BCD; AB-CD; ABC-D)
  n.pairs <- length(fcast.values) - 1
  
  pairs <- array(NA, dim=c(n.pairs,2))
  for(i in 1:n.pairs){
    forecast.sub1 <- sum(fcast.values[1:i])
    forecast.sub2 <- sum(fcast.values[-(1:i)])
    pairs[i,] <- c(forecast.sub1,forecast.sub2)
  }
  
  # Step 1a define ref.value (correct option) for each pair
  correct.arr <- array(NA, dim=dim(pairs))
  for(i in 1:n.pairs){
    ref.sub1 <- sum(correct.values[1:i])
    ref.sub2 <- sum(correct.values[-(1:i)])
    correct.arr[i,] <- c(ref.sub1,ref.sub2)
  } 
  if(max(correct.arr) > 1){
    stop(paste(Sys.time(),"BrierScoreOrdered: ref.vals calculation is flawed, values > 1 [",
               paste(correct.arr,collapse=","),"]"))
  }
  
  # Step 2 – Apply the scoring rule to each binary pair
  n.pairs <- dim(pairs)[1]
  score.pairs <- array(dim=n.pairs)
  for(i in 1:n.pairs){
    score.pairs[i] <- scoring.fn(fcast.values   = pairs[i,], 
                                 correct.values = correct.arr[i,])
  }
  
  # Step 3 - Take an average across the binary pair scores.
  score <- mean(score.pairs)   
  return(score)
}