scoringRules (version 1.1.3)

scores_sample_multiv_weighted: Weighted Multivariate Scoring Rules for Simulated Forecast Distributions

Description

Compute weighted versions of multivariate scores \(S(y, dat)\), where \(S\) is a proper scoring rule, \(y\) is a d-dimensional realization vector and \(dat\) is a simulated sample of multivariate forecasts. The weighted scores allow particular outcomes of interest to be emphasised during forecast evaluation. Threshold-weighted and outcome-weighted versions of three multivariate scores are available: the energy score, a score based on a Gaussian kernel (mmds_sample, see details below) and the variogram score of order \(p\).

Usage

twes_sample(y, dat, a = -Inf, b = Inf, chain_func = NULL, w = NULL)
owes_sample(y, dat, a = -Inf, b = Inf, weight_func = NULL, w = NULL)
twmmds_sample(y, dat, a = -Inf, b = Inf, chain_func = NULL, w = NULL)
owmmds_sample(y, dat, a = -Inf, b = Inf, weight_func = NULL, w = NULL)
twvs_sample(
  y,
  dat,
  a = -Inf,
  b = Inf,
  chain_func = NULL,
  w = NULL,
  w_vs = NULL,
  p = 0.5
)
owvs_sample(
  y,
  dat,
  a = -Inf,
  b = Inf,
  weight_func = NULL,
  w = NULL,
  w_vs = NULL,
  p = 0.5
)

Value

Value of the score. A lower score indicates a better forecast.

Arguments

y

realized values (numeric vector of length d).

dat

numeric matrix of data (columns are simulation draws from multivariate forecast distribution).

a

numeric vector of of length d containing lower bounds for the indicator weight function w(z) = 1{a[1] < z[1] < b[1], ..., a[d] < z[d] < b[d]}.

b

numeric vector of of length d containing upper bounds for the indicator weight function w(z) = 1{a[1] < z[1] < b[1], ..., a[d] < z[d] < b[d]}.

chain_func

function used to target particular outcomes in the threshold-weighted scores; the default corresponds to the weight function above.

w

numeric vector of weights for forecast draws (length equal to number of columns of dat)

weight_func

function used to target particular outcomes in the outcome-weighted scores; the default corresponds to the weight function above.

w_vs

numeric matrix of weights for dat used in the variogram score. This matrix must be square and symmetric, with all elements being non-negative. If no weights are specified, constant weights (with all elements of w_vs equal to one) are used.

p

order of variogram score. Standard choices include \(p = 1\) and \(p = 0.5\).

Author

Sam Allen

Details

In the input matrix dat each column is expected to represent a sample from the multivariate forecast distribution, the number of rows of dat thus has to match the length of the observation vector y, and the number of columns of dat is the number of simulated samples.

The threshold-weighted scores (twes_sample, twmmds_sample, twvs_sample) transform y and dat using the chaining function chain_func and then call the relevant unweighted score function (es_sample, mmds_sample, vs_sample). The outcome-weighted scores (owes_sample, owmmds_sample, owvs_sample) weight y and dat using the weight function weight_func and then call the relevant unweighted score function (es_sample, mmds_sample, vs_sample). See the documentation for e.g. es_sample for further details.

The default weight function used in the weighted scores is w(z) = 1{a[1] < z[1] < b[1], ..., a[d] < z[d] < b[d]}, which is equal to one if z is in the 'box' defined by the vectors a and b, and is equal to zero otherwise. This weight function emphasises outcomes between the vectors a and b, and is commonly used in practical applications when interest is on values above a threshold along multiple dimensions.

Alternative weight functions can also be employed using the chain_func and weight_func arguments. Computation of the threshold-weighted scores for samples from a predictive distribution requires a chaining function rather than a weight function. This is why a chaining function is an input for twes_sample, twmmds_sample, and twvs_sample, whereas a weight function is an input for owes_sample, owmmds_sample, and owvs_sample.

The chain_func and weight_func arguments are functions that will be applied to the elements in y and dat. weight_func must input a numeric vector of length d, and output a single numeric value. An error will be returned if weight_func returns negative values. chain_func must input a numeric vector of length d, and return a numeric vector of length d.

If no custom argument is given for a, b, chain_func or weight_func, then all weighted scores are equivalent to the standard unweighted scores es_sample, mmds_sample, and vs_sample.

The w argument is also present in the unweighted scores. w is used to weight the draws from the predictive distribution, and does not weight particular outcomes within the weighted scoring rules. This should not be confused with the weight_func argument.

References

Allen, S. (2024): `Weighted scoringRules: Emphasising Particular Outcomes when Evaluating Probabilistic Forecasts', Journal of Statistical Software. tools:::Rd_expr_doi("10.18637/jss.v110.i08")

Threshold-weighted scores

Allen, S., Ginsbourger, D. and J. Ziegel (2023): `Evaluating forecasts for high-impact events using transformed kernel scores', SIAM/ASA Journal on Uncertainty Quantification 11, 906-940. tools:::Rd_expr_doi("10.1137/22M1532184")

Outcome-weighted scores:

Holzmann, H. and B. Klar (2017): `Focusing on regions of interest in forecast evaluation', Annals of Applied Statistics 11, 2404-2431. tools:::Rd_expr_doi("10.1214/17-AOAS1088")

See Also

scores_sample_multiv for standard (unweighted) scores based on simulated multivariate forecast distributions. scores_sample_univ_weighted for weighted scores based on simulated univariate forecast distributions

Examples

if (FALSE) {
d <- 3  # number of dimensions
m <- 10  # number of samples from multivariate forecast distribution

# parameters for multivariate normal example
mu0 <- rep(0, d)
mu <- rep(1, d)
S0 <- S <- diag(d)
S0[S0==0] <- 0.2
S[S==0] <- 0.1

# generate samples from multivariate normal distributions
obs <- drop(mu0 + rnorm(d) %*% chol(S0))
sample_fc <- replicate(m, drop(mu + rnorm(d) %*% chol(S)))

# if no additional parameters are provided, the weighted scores are the same as
# the unweighted scores:
es_sample(y = obs, dat = sample_fc) # energy score
twes_sample(y = obs, dat = sample_fc) # threshold-weighted energy score
owes_sample(y = obs, dat = sample_fc) # outcome-weighted energy score

mmds_sample(y = obs, dat = sample_fc) # Gaussian kernel score
twmmds_sample(y = obs, dat = sample_fc) # threshold-weighted Gaussian kernel score
owmmds_sample(y = obs, dat = sample_fc) # outcome-weighted Gaussian kernel score

vs_sample(y = obs, dat = sample_fc) # variogram score
twvs_sample(y = obs, dat = sample_fc) # threshold-weighted variogram score
owvs_sample(y = obs, dat = sample_fc) # outcome-weighted variogram score


# the outcome-weighted scores are undefined if none of dat are between a and b
# this can lead to NaNs in some of the scores calculated below, particularly
# if the thresholds are extreme, or if the dimension is large


# emphasise outcomes greater than 0 in all dimensions
twes_sample(y = obs, dat = sample_fc, a = 0)
owes_sample(y = obs, dat = sample_fc, a = 0)
twmmds_sample(y = obs, dat = sample_fc, a = 0)
owmmds_sample(y = obs, dat = sample_fc, a = 0)
twvs_sample(y = obs, dat = sample_fc, a = 0)
owvs_sample(y = obs, dat = sample_fc, a = 0)

# this can also be done more explicitly by setting a = rep(0, d)
twes_sample(y = obs, dat = sample_fc, a = rep(0, d))
owes_sample(y = obs, dat = sample_fc, a = rep(0, d))

# a should also be specified fully if the threshold changes in each dimension
a <- rnorm(d)
twes_sample(y = obs, dat = sample_fc, a = a)
owes_sample(y = obs, dat = sample_fc, a = a)
twmmds_sample(y = obs, dat = sample_fc, a = a)
owmmds_sample(y = obs, dat = sample_fc, a = a)
twvs_sample(y = obs, dat = sample_fc, a = a)
owvs_sample(y = obs, dat = sample_fc, a = a)

# emphasise outcomes smaller than 0 in all dimensions
twes_sample(y = obs, dat = sample_fc, b = 0)
owes_sample(y = obs, dat = sample_fc, b = 0)
twmmds_sample(y = obs, dat = sample_fc, b = 0)
owmmds_sample(y = obs, dat = sample_fc, b = 0)
twvs_sample(y = obs, dat = sample_fc, b = 0)
owvs_sample(y = obs, dat = sample_fc, b = 0)

# emphasise outcomes between (-1, -1, -1) and (1, 1, 1)
twes_sample(y = obs, dat = sample_fc, a = -1, b = 1)
owes_sample(y = obs, dat = sample_fc, a = -1, b = 1)
twmmds_sample(y = obs, dat = sample_fc, a = -1, b = 1)
owmmds_sample(y = obs, dat = sample_fc, a = -1, b = 1)
twvs_sample(y = obs, dat = sample_fc, a = -1, b = 1)
owvs_sample(y = obs, dat = sample_fc, a = -1, b = 1)

# emphasise outcomes between (-2, 0, -1) and (0, 2, 1)
a <- c(-2, 0, -1)
b <- c(0, 2, 1)
twes_sample(y = obs, dat = sample_fc, a = a, b = b)
owes_sample(y = obs, dat = sample_fc, a = a, b = b)
twmmds_sample(y = obs, dat = sample_fc, a = a, b = b)
owmmds_sample(y = obs, dat = sample_fc, a = a, b = b)
twvs_sample(y = obs, dat = sample_fc, a = a, b = b)
owvs_sample(y = obs, dat = sample_fc, a = a, b = b)


# values of a cannot be larger than the corresponding values of b
twes_sample(y = obs, dat = sample_fc, a = c(0, 0, 0), b = c(0, 0, 1))
twes_sample(y = obs, dat = sample_fc, a = c(0, 0, 0), b = c(0, 0, 0)) # error
twes_sample(y = obs, dat = sample_fc, a = c(0, 0, 0), b = c(1, 1, -1)) # error

# a and b must be of the same length (and of the same length as y)
owmmds_sample(y = obs, dat = sample_fc, a = c(0, 0), b = 1) # error
owmmds_sample(y = obs, dat = sample_fc, a = c(0, 0), b = c(1, 1)) # error


# alternative custom weight and chaining functions can also be used

# Example 1: the default weight function with an alternative chaining function
# the default weight function is 
# w(z) = 1{a[1] < z[1] < b[1], ..., a[d] < z[d] < b[d]}
# the default chaining function is 
# v(z) = (min(max(z[1], a[1]), b[1]), ..., min(max(z[d], a[d]), b[d]))
a <- -2
b <- 2
weight_func <- function(x) as.numeric(all(x > a & x < b))
chain_func <- function(x) pmin(pmax(x, a), b)
owes_sample(y = obs, dat = sample_fc, a = a, b = b)
owes_sample(y = obs, dat = sample_fc, weight_func = weight_func)
twes_sample(y = obs, dat = sample_fc, a = a, b = b)
twes_sample(y = obs, dat = sample_fc, chain_func = chain_func)

# consider an alternative chaining function: v(z) = z if w(z) = 1, else v(z) = 0
chain_func <- function(x) x*weight_func(x)
twes_sample(y = obs, dat = sample_fc, chain_func = chain_func)


# Example 2: a mulivariate Gaussian weight function with mean vector mu and 
# diagonal covariance matrix sigma
mu <- rep(0, d)
sigma <- diag(d)
weight_func <- function(x) prod(pnorm(x, mu, diag(sigma)))
# the corresponding chaining function is
chain_func <- function(x){
 (x - mu)*pnorm(x, mu, diag(sigma)) + (diag(sigma)^2)*dnorm(x, mu, diag(sigma))
}

owvs_sample(y = obs, dat = sample_fc, a = mu)
owvs_sample(y = obs, dat = sample_fc, weight_func = weight_func)
twvs_sample(y = obs, dat = sample_fc, a = mu)
twvs_sample(y = obs, dat = sample_fc, chain_func = chain_func)
}