0
# Continuous Scores
1

2
Metrics for evaluating single-valued continuous forecasts including error measures, bias calculations, efficiency indices, and correlation coefficients. All functions support flexible dimension handling, optional weighting, and angular data processing.
3

4
## Capabilities
5

6
### Mean Squared Error (MSE)
7

8
Calculates the mean squared error between forecast and observations, the most fundamental error metric in forecast verification.
9

10
```python { .api }
11
def mse(
12
    fcst: FlexibleArrayType,
13
    obs: FlexibleArrayType,
14
    *,
15
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
16
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
17
    weights: Optional[xr.DataArray] = None,
18
    is_angular: Optional[bool] = False,
19
) -> XarrayLike:
20
    """
21
    Calculate Mean Squared Error.
22

23
    Args:
24
        fcst: Forecast data
25
        obs: Observation data
26
        reduce_dims: Dimensions to reduce (collapse)
27
        preserve_dims: Dimensions to preserve (all others reduced)
28
        weights: Optional weights for weighted MSE calculation
29
        is_angular: Handle circular data (e.g. wind direction)
30

31
    Returns:
32
        Mean squared error values
33

34
    Formula:
35
        MSE = (1/n) * Σ(forecast_i - observed_i)²
36
    """
37
```
38

39
**Usage Example:**
40

41
```python
42
from scores.continuous import mse
43
import xarray as xr
44

45
# Basic MSE calculation
46
mse_value = mse(forecast, observations)
47

48
# MSE with dimension reduction
49
temporal_mse = mse(forecast, observations, reduce_dims="time")
50

51
# Weighted MSE (e.g. area weighting)
52
weighted_mse = mse(forecast, observations, weights=area_weights)
53

54
# Angular MSE for directional data
55
wind_dir_mse = mse(forecast_dir, observed_dir, is_angular=True)
56
```
57

58
### Root Mean Squared Error (RMSE)
59

60
Square root of MSE, providing error values in the same units as the original data.
61

62
```python { .api }
63
def rmse(
64
    fcst: FlexibleArrayType,
65
    obs: FlexibleArrayType,
66
    *,
67
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
68
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
69
    weights: Optional[xr.DataArray] = None,
70
    is_angular: Optional[bool] = False,
71
) -> FlexibleArrayType:
72
    """
73
    Calculate Root Mean Squared Error.
74

75
    Args:
76
        fcst: Forecast data
77
        obs: Observation data
78
        reduce_dims: Dimensions to reduce
79
        preserve_dims: Dimensions to preserve
80
        weights: Optional weights
81
        is_angular: Handle circular data
82

83
    Returns:
84
        Root mean squared error values
85

86
    Formula:
87
        RMSE = √[(1/n) * Σ(forecast_i - observed_i)²]
88
    """
89
```
90

91
### Mean Absolute Error (MAE)
92

93
Mean absolute deviation between forecast and observations, less sensitive to outliers than MSE.
94

95
```python { .api }
96
def mae(
97
    fcst: FlexibleArrayType,
98
    obs: FlexibleArrayType,
99
    *,
100
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
101
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
102
    weights: Optional[xr.DataArray] = None,
103
    is_angular: Optional[bool] = False,
104
) -> FlexibleArrayType:
105
    """
106
    Calculate Mean Absolute Error.
107

108
    Args:
109
        fcst: Forecast data
110
        obs: Observation data
111
        reduce_dims: Dimensions to reduce
112
        preserve_dims: Dimensions to preserve
113
        weights: Optional weights
114
        is_angular: Handle circular data
115

116
    Returns:
117
        Mean absolute error values
118

119
    Formula:
120
        MAE = (1/n) * Σ|forecast_i - observed_i|
121
    """
122
```
123

124
### Bias Metrics
125

126
#### Additive Bias (Mean Error)
127

128
Average difference between forecast and observations, indicating systematic over- or under-prediction.
129

130
```python { .api }
131
def additive_bias(
132
    fcst: XarrayLike,
133
    obs: XarrayLike,
134
    *,
135
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
136
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
137
    weights: Optional[XarrayLike] = None,
138
    is_angular: Optional[bool] = False,
139
) -> XarrayLike:
140
    """
141
    Calculate additive bias (mean error).
142

143
    Args:
144
        fcst: Forecast data
145
        obs: Observation data
146
        reduce_dims: Dimensions to reduce
147
        preserve_dims: Dimensions to preserve
148
        weights: Optional weights
149
        is_angular: Handle circular data
150

151
    Returns:
152
        Additive bias values
153

154
    Formula:
155
        Bias = (1/n) * Σ(forecast_i - observed_i)
156
    """
157

158
# Alias for additive_bias
159
def mean_error(
160
    fcst: XarrayLike,
161
    obs: XarrayLike,
162
    *,
163
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
164
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
165
    weights: Optional[XarrayLike] = None,
166
    is_angular: Optional[bool] = False,
167
) -> XarrayLike:
168
    """Alias for additive_bias."""
169
```
170

171
#### Multiplicative Bias
172

173
Ratio-based bias measure useful for positive-valued variables.
174

175
```python { .api }
176
def multiplicative_bias(
177
    fcst: XarrayLike,
178
    obs: XarrayLike,
179
    *,
180
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
181
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
182
    weights: Optional[XarrayLike] = None,
183
) -> XarrayLike:
184
    """
185
    Calculate multiplicative bias.
186

187
    Args:
188
        fcst: Forecast data (positive values)
189
        obs: Observation data (positive values)
190
        reduce_dims: Dimensions to reduce
191
        preserve_dims: Dimensions to preserve
192
        weights: Optional weights
193

194
    Returns:
195
        Multiplicative bias values
196

197
    Formula:
198
        Multiplicative Bias = mean(forecast) / mean(observation)
199

200
    Notes:
201
        - Values > 1 indicate forecast over-prediction
202
        - Values < 1 indicate forecast under-prediction
203
        - Perfect forecast has multiplicative bias = 1
204
    """
205
```
206

207
#### Percent Bias (PBIAS)
208

209
Percentage-based bias calculation commonly used in hydrology.
210

211
```python { .api }
212
def pbias(
213
    fcst: XarrayLike,
214
    obs: XarrayLike,
215
    *,
216
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
217
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
218
    weights: Optional[XarrayLike] = None,
219
) -> XarrayLike:
220
    """
221
    Calculate percent bias.
222

223
    Args:
224
        fcst: Forecast data
225
        obs: Observation data
226
        reduce_dims: Dimensions to reduce
227
        preserve_dims: Dimensions to preserve
228
        weights: Optional weights
229

230
    Returns:
231
        Percent bias values
232

233
    Formula:
234
        PBIAS = 100 * Σ(forecast_i - observed_i) / Σ(observed_i)
235

236
    Notes:
237
        - Positive values indicate model overestimation bias
238
        - Negative values indicate model underestimation bias
239
        - Perfect forecast has PBIAS = 0
240
    """
241
```
242

243
### Efficiency Indices
244

245
#### Kling-Gupta Efficiency (KGE)
246

247
Decomposed performance metric combining correlation, variability ratio, and bias ratio.
248

249
```python { .api }
250
def kge(
251
    fcst: xr.DataArray,
252
    obs: xr.DataArray,
253
    *,
254
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
255
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
256
    scaling_factors: Optional[Union[list[float], np.ndarray]] = None,
257
    include_components: Optional[bool] = False,
258
) -> XarrayLike:
259
    """
260
    Calculate Kling-Gupta Efficiency.
261

262
    Args:
263
        fcst: Forecast data
264
        obs: Observation data
265
        reduce_dims: Dimensions to reduce
266
        preserve_dims: Dimensions to preserve
267
        scaling_factors: Component scaling factors [ρ, α, β]
268
        include_components: Return individual components
269

270
    Returns:
271
        KGE efficiency values (and components if requested)
272

273
    Formula:
274
        KGE = 1 - √[(s_ρ*(ρ-1))² + (s_α*(α-1))² + (s_β*(β-1))²]
275

276
    Components:
277
        - ρ: Correlation coefficient
278
        - α: Variability ratio (σ_forecast/σ_observation)
279
        - β: Bias ratio (μ_forecast/μ_observation)
280
        - s_ρ, s_α, s_β: Scaling factors (default: [1,1,1])
281

282
    Notes:
283
        - Perfect forecast has KGE = 1
284
        - Default benchmark (observation mean) has KGE ≈ -0.41
285
    """
286
```
287

288
#### Nash-Sutcliffe Efficiency (NSE)
289

290
Normalized measure of model performance relative to observation mean.
291

292
```python { .api }
293
def nse(
294
    fcst: XarrayLike,
295
    obs: XarrayLike,
296
    *,
297
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
298
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
299
    weights: Optional[XarrayLike] = None,
300
    is_angular: Optional[bool] = False,
301
) -> XarrayLike:
302
    """
303
    Calculate Nash-Sutcliffe Efficiency.
304

305
    Args:
306
        fcst: Forecast data
307
        obs: Observation data
308
        reduce_dims: Dimensions to reduce
309
        preserve_dims: Dimensions to preserve
310
        weights: Optional weights
311
        is_angular: Handle circular data
312

313
    Returns:
314
        Nash-Sutcliffe efficiency values
315

316
    Formula:
317
        NSE = 1 - Σ(observed_i - forecast_i)² / Σ(observed_i - observed_mean)²
318

319
    Notes:
320
        - Perfect forecast has NSE = 1
321
        - Forecast as good as observation mean has NSE = 0
322
        - NSE can be negative (worse than observation mean)
323
        - Range: (-∞, 1]
324
    """
325
```
326

327
### Correlation Coefficients
328

329
#### Pearson Correlation
330

331
Linear correlation coefficient measuring strength of linear relationship.
332

333
```python { .api }
334
def pearsonr(
335
    fcst: XarrayLike,
336
    obs: XarrayLike,
337
    *,
338
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
339
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
340
) -> XarrayLike:
341
    """
342
    Calculate Pearson correlation coefficient.
343

344
    Args:
345
        fcst: Forecast data
346
        obs: Observation data
347
        reduce_dims: Dimensions to reduce
348
        preserve_dims: Dimensions to preserve
349

350
    Returns:
351
        Pearson correlation coefficients
352

353
    Notes:
354
        - Range: [-1, 1]
355
        - Measures linear association strength
356
        - Perfect positive correlation: r = 1
357
        - Perfect negative correlation: r = -1
358
        - No linear relationship: r = 0
359
        - This function doesn't support weights
360
    """
361
```
362

363
#### Spearman Rank Correlation
364

365
Non-parametric correlation based on rank ordering.
366

367
```python { .api }
368
def spearmanr(
369
    fcst: XarrayLike,
370
    obs: XarrayLike,
371
    *,
372
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
373
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
374
) -> XarrayLike:
375
    """
376
    Calculate Spearman rank correlation coefficient.
377

378
    Args:
379
        fcst: Forecast data
380
        obs: Observation data
381
        reduce_dims: Dimensions to reduce
382
        preserve_dims: Dimensions to preserve
383

384
    Returns:
385
        Spearman rank correlation coefficients
386

387
    Notes:
388
        - Range: [-1, 1]
389
        - Measures monotonic association strength
390
        - Less sensitive to outliers than Pearson
391
        - This function doesn't support weights
392
        - Based on rank ordering rather than actual values
393
    """
394
```
395

396
### Advanced Scoring Functions
397

398
#### Quantile Score
399

400
Asymmetric loss function for evaluating quantile forecasts.
401

402
```python { .api }
403
def quantile_score(
404
    fcst: FlexibleArrayType,
405
    obs: FlexibleArrayType,
406
    alpha: float,
407
    *,
408
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
409
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
410
    weights: Optional[xr.DataArray] = None,
411
) -> FlexibleArrayType:
412
    """
413
    Calculate quantile score for quantile forecasts.
414

415
    Args:
416
        fcst: Quantile forecast values
417
        obs: Observation values
418
        alpha: Quantile level (e.g., 0.1 for 10th percentile)
419
        reduce_dims: Dimensions to reduce
420
        preserve_dims: Dimensions to preserve
421
        weights: Optional weights
422

423
    Returns:
424
        Quantile scores
425

426
    Formula:
427
        QS = Σ[(α - I(obs_i < fcst_i)) * (fcst_i - obs_i)]
428

429
    Notes:
430
        - Proper scoring rule for quantile forecasts
431
        - Alpha should match the quantile level being forecast
432
        - Lower scores indicate better performance
433
    """
434
```
435

436
#### Flip-Flop Index
437

438
Measures forecast consistency by detecting oscillatory behavior.
439

440
```python { .api }
441
def flip_flop_index(
442
    fcst: xr.DataArray,
443
    obs: xr.DataArray,
444
    *,
445
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
446
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
447
) -> xr.DataArray:
448
    """
449
    Calculate Flip-Flop Index for forecast consistency.
450

451
    Args:
452
        fcst: Forecast data (must have time dimension)
453
        obs: Observation data
454
        reduce_dims: Dimensions to reduce
455
        preserve_dims: Dimensions to preserve
456

457
    Returns:
458
        Flip-flop index values
459

460
    Notes:
461
        - Measures forecast oscillatory behavior
462
        - Higher values indicate more inconsistent forecasts
463
        - Requires time series data
464
        - Used to detect forecast instability
465
    """
466

467
def flip_flop_index_proportion_exceeding(
468
    fcst: xr.DataArray,
469
    obs: xr.DataArray,
470
    threshold: float,
471
    *,
472
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
473
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
474
) -> xr.DataArray:
475
    """
476
    Calculate proportion of flip-flop indices exceeding threshold.
477

478
    Args:
479
        fcst: Forecast data
480
        obs: Observation data
481
        threshold: Threshold value for comparison
482
        reduce_dims: Dimensions to reduce
483
        preserve_dims: Dimensions to preserve
484

485
    Returns:
486
        Proportion values [0, 1]
487
    """
488
```
489

490
## Usage Patterns
491

492
### Standard Workflow
493

494
```python
495
from scores.continuous import mse, rmse, mae, kge, pearsonr
496
from scores.sample_data import continuous_forecast, continuous_observations
497

498
# Load data
499
forecast = continuous_forecast()
500
observations = continuous_observations()
501

502
# Calculate multiple metrics
503
scores_dict = {
504
    'MSE': mse(forecast, observations),
505
    'RMSE': rmse(forecast, observations),
506
    'MAE': mae(forecast, observations),
507
    'KGE': kge(forecast, observations),
508
    'Correlation': pearsonr(forecast, observations)
509
}
510

511
# Display results
512
for metric, value in scores_dict.items():
513
    print(f"{metric}: {value.values:.4f}")
514
```
515

516
### Dimension-aware Scoring
517

518
```python
519
# Multi-dimensional data (time, lat, lon)
520
forecast_3d = forecast.expand_dims({"lat": 10, "lon": 15})
521
observations_3d = observations.expand_dims({"lat": 10, "lon": 15})
522

523
# Compute temporal scores at each grid point
524
spatial_temporal_mse = mse(forecast_3d, observations_3d, reduce_dims="time")
525

526
# Compute spatial scores for each time step
527
temporal_spatial_mse = mse(forecast_3d, observations_3d, reduce_dims=["lat", "lon"])
528

529
# Overall score (all dimensions)
530
overall_mse = mse(forecast_3d, observations_3d)
531
```
532

533
### Weighted Scoring
534

535
```python
536
from scores.functions import create_latitude_weights
537
import numpy as np
538

539
# Create latitude weights for area weighting
540
lats = np.linspace(-90, 90, forecast_3d.sizes["lat"])
541
area_weights = create_latitude_weights(lats)
542

543
# Calculate area-weighted scores
544
weighted_mse = mse(forecast_3d, observations_3d,
545
                   reduce_dims=["time", "lat", "lon"],
546
                   weights=area_weights)
547
```