Verification and SpecsVerification

require(SpecsVerification)

## Loading required package: SpecsVerification

require(verification) ##NCAR - daje jako lepe rezultate 

## Loading required package: verification

## Loading required package: fields

## Loading required package: spam

## Loading required package: grid

## Spam version 1.4-0 (2016-08-29) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

## 
## Attaching package: 'spam'

## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve

## Loading required package: maps

## Loading required package: boot

## Loading required package: CircStats

## Loading required package: MASS

## Loading required package: dtw

## Loading required package: proxy

## 
## Attaching package: 'proxy'

## The following object is masked from 'package:spam':
## 
##     as.matrix

## The following objects are masked from 'package:stats':
## 
##     as.dist, dist

## The following object is masked from 'package:base':
## 
##     as.matrix

## Loaded dtw v1.18-1. See ?dtw for help, citation("dtw") for use in publication.

require(ggplot2)

## Loading required package: ggplot2

rm(list=ls())

First we need an ensamble to work with

There are two options basically - generate an ensemble around a set of obervations or generate the obs and the ens together

The ClimEns function a climatological ENS around a given set of observations. The funciton can be constructed so that the output does not contain the observation

obs <- rnorm(10) ## time series of 10 observed values 
ens<-ClimEns(obs,leave.one.out = FALSE) ##calculates a 10 member ens around each value 

The GenerateToyData function generates a synthetic ensemble and observation data, using a latent variable model. See Details.

The function takes the following arguments

args(GenerateToyData)

## function (N = 20, mu.y = 0, s.s = 7, s.eps = 6, mu.x = 0, beta = 0.2, 
##     s.eta = 8, K = 10, mu.x.ref = NA, beta.ref = NA, s.eta.ref = NA, 
##     K.ref = NA) 
## NULL

ens_toy<-GenerateToyData()
ens<-ens_toy$ens
obs<-ens_toy$obs

1. Rank histogram

PlotRankhist(Rankhist(ens,obs))

rank.hist <- Rankhist(ens,obs) 
PlotRankhist(rank.hist, mode="prob.paper")

2. Brier score

I guess it makes sense to use the Brier even when it’s not a binary event like rain/did not rain. Take the median or the average for the observations as a threshold.

EnsBrier(ens, obs, tau) - SpecsVerification

ens - N*K matrix. N cases of a K-member ensemble forecast. obs - vector of length N. The corresponding verifying observations.

tau -The threshold value whose exceedance defines the binary event. If ensemble members and observations are coded as 0???s and 1???s already, setting tau=0.5 (the default) will produce the desired result. Value

mean(EnsBrier(ens,obs,tau=0)) ## This is the Brier score 

## [1] 0.2735

brier(obs, pred, baseline, thresholds = seq(0,1,0.1), bins = TRUE, … ) - verification

obs - vector of binary observation pred - vector of probabilistic predictins [0,1]

obs<-rbinom(n=100,size=1,prob=0.5) ## size sets the range 
pred<-rbeta(n=100,shape1=1,shape2=3)
brier(obs,pred)

## $baseline.tf
## [1] FALSE
## 
## $bs
## [1] 0.3369
## 
## $bs.baseline
## [1] 0.2499
## 
## $ss
## [1] -0.3481393
## 
## $bs.reliability
## [1] 0.1049585
## 
## $bs.resol
## [1] 0.01795848
## 
## $bs.uncert
## [1] 0.2499
## 
## $y.i
##  [1] 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95
## 
## $obar.i
##  [1] 0.3666667 0.6842105 0.5000000 0.5000000 0.4444444 0.4000000 0.6000000
##  [8] 1.0000000 0.0000000        NA
## 
## $prob.y
##  [1] 0.30 0.19 0.12 0.18 0.09 0.05 0.05 0.01 0.01   NA
## 
## $obar
## [1] 0.49
## 
## $thres
##  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
## 
## $check
## [1] 0.3369
## 
## $bins
## [1] TRUE

EnsBrierSs(ens, ens.ref, obs, tau)

Calculate Brier skill score or continuous ranked probability skill score between two competing ensemble forecasts for the same observation. The sampling standard deviation is estimated by propagation of uncertainty. The higher the skill score, the higher the improvement of ???ens??? over ???ens.ref???.

dat.1<-GenerateToyData()
ens.1<-dat.1$ens
obs.1<-dat.1$obs

dat.ref<-GenerateToyData() ##usually the reference is the climatology or in our case kvasi climatology 
ens.ref<-dat.ref$ens
obs.ref<-dat.ref$obs

EnsBrierSs(ens.1,ens.ref,obs.ref,1) 

## $bss
## [1] -0.2440087
## 
## $bss.sigma
## [1] 0.1820807

BrierScoreDecomposition(p, y, calibration=list(method=“bin”, bins=10), probs=NA, n.boot=0)

Decomposition of the empirical Brier Score into its Reliability, Resolution, and Uncertainty com- ponent

Najvece pitanje je kako da iz ensambla mozemo u R/u da izracunamo verovatnoce

p <- rbeta(n=100, shape1=1, shape2=3) 
y <- rbinom(n=100, size=1, prob=p^0.9) 
BrierScoreDecomposition(p, y, calibration=list(method="bin", bins=5))

##         REL         RES         UNC 
## 0.009639413 0.029888312 0.182400000

================================================ SIDEBAR ================================================

The beta distribution is basically a standard distribution for probabilistic values [0-1]

dbeta - density pbeta - pdf qbeta - quantile rbeta - randombeta

Shape1 je alfa, shape2 je beta URL wiki

and the r- prefix is always a random number generator for some distribution or other rbinom rnorm; d - for a density function dbinom dnorm, p for a distribution function (PDF) pbinom pnorm