feat: Add Scaled CRPS as an accuracy metric

CHANGELOG.md

@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.13.0] ; 2025-11-07
+
+- Add Scaled CRPS as an accuracy metric.
+
 ## [0.12.0] ; 2025-07-23
 
 ### Added

docs/user-guide/metrics.md

@@ -16,6 +16,7 @@ mmm-eval calculates several key metrics across different validation tests:
 - **MAPE (Mean Absolute Percentage Error)**: Average percentage error between predictions and actual values
 - **SMAPE (Symmetric Mean Absolute Percentage Error)**: Symmetric version of MAPE that treats over and underestimation equally
 - **R-squared**: Proportion of variance explained by the model
+- **CRPS (Continuous Ranked Probability Score)**: Distributional accuracy metric
 
 ### Stability Metrics
 
@@ -63,6 +64,19 @@ R² = 1 - (Σ(y_i - ŷ_i)² / Σ(y_i - ȳ)²)
 - **Scale**: 1 = perfect fit, 0 = no predictive power
 - **Benchmark**: > 0.8 is generally good
 
+### CRPS (Continuous Ranked Probability Score)
+
+```python
+CRPS = ((1/n) * Σ|y_i - ŷ| - (0.5/n) * Σ|y_i - x_i|) / y_i)
+```
+where `x_i` and `y_i` are samples from the posterior predictive distribution
+
+
+**Interpretation**:
+- **Lower is better**: 0% = perfect predictions, no uncertainty
+- **Scale**: Expressed as a percentage, e.g. 15.0 rather than 0.15
+
+
 ## Test-Specific Metrics
 
 ### Holdout Accuracy Test Metrics
@@ -72,6 +86,7 @@ Metrics calculated on out-of-sample predictions using train/test splits.
 - **MAPE**: Overall prediction accuracy
 - **SMAPE**: Symmetric prediction accuracy
 - **R-squared**: Model fit quality
+- **CRPS**: Distribution prediction accuracy
 
 ### In-Sample Accuracy Test Metrics
 
@@ -80,6 +95,8 @@ Metrics calculated on in-sample predictions using the full dataset.
 - **MAPE**: Model fit accuracy
 - **SMAPE**: Symmetric model fit accuracy
 - **R-squared**: Model fit quality
+- **CRPS**: Distribution prediction accuracy
+
 
 ### Cross-Validation Metrics
 
@@ -89,6 +106,8 @@ Metrics calculated on in-sample predictions using the full dataset.
 - **Std SMAPE**: Consistency of symmetric accuracy across folds
 - **Mean R-squared**: Average out-of-sample fit
 - **Std R-squared**: Consistency of fit across folds
+- **Mean CRPS**: Average distribution prediction accuracy
+
 
 ### Refresh Stability Metrics
 
@@ -118,6 +137,7 @@ model picking up spurious correlations.
 - **MAPE < 15%**: Good prediction accuracy
 - **SMAPE < 15%**: Good symmetric prediction accuracy
 - **R-squared > 0.8**: Strong model fit
+- **CRPS < 15%**: Good distributional accuracy
 - **Low perturbation sensitivity**: Robust to input noise
 - **Low placebo ROI (≤ -50%)**: Correctly identifies spurious features and assigns them low effect sizes
 
@@ -134,6 +154,7 @@ model specification to the problem at hand.
 | MAPE | < 5% | 5-10% | 10-15% | > 15% |
 | SMAPE | < 5% | 5-10% | 10-15% | > 15% |
 | R-squared | > 0.9 | 0.8-0.9 | 0.6-0.8 | < 0.6 |
+| CRPS | < 5% | 5-10% | 10-15% | > 15% |
 | Parameter Change | < 5% | 5-10% | 10-20% | > 20% |
 | Perturbation Change | < 5% | 5-10% | 10-15% | > 15% |
 | Placebo ROI | ≤ -50% | -50% to -25% | -25% to 0% | > 0% |
@@ -165,6 +186,7 @@ class CustomMetric(BaseMetric):
 - **Start with MAPE**: Most intuitive for business users
 - **Include SMAPE**: More robust alternative to MAPE for symmetric evaluation
 - **Include R-squared**: Technical measure of fit quality
+- **Add CRPS**: Measure the uncertainty of the predictions
 - **Monitor stability**: Critical for production models
 - **Track performance**: Important for scalability
 

docs/user-guide/tests.md

@@ -39,12 +39,14 @@ The holdout accuracy test evaluates model performance by splitting data into tra
 - **MAPE (Mean Absolute Percentage Error)**: Average percentage error
 - **SMAPE (Symmetric Mean Absolute Percentage Error)**: Symmetric version of MAPE
 - **R-squared**: Proportion of variance explained by the model
+- **CRPS (Continuous Ranked Probability Score)**: Distribution accuracy metric, scaled to sales levels.
 
 #### Interpretation
 
 - **Lower MAPE**: Better model performance
 - **Lower SMAPE**: Better symmetric model performance
 - **Higher R-squared**: Better model fit (0-1 scale)
+- **Lower CRPS**: Less uncertainty
 
 ### In-Sample Accuracy Test
 
@@ -61,12 +63,14 @@ The in-sample accuracy test evaluates model performance by fitting the model on
 - **MAPE (Mean Absolute Percentage Error)**: Average percentage error
 - **SMAPE (Symmetric Mean Absolute Percentage Error)**: Symmetric version of MAPE
 - **R-squared**: Proportion of variance explained by the model
+- **CRPS (Continuous Ranked Probability Score)**: Distribution accuracy metric, scaled to sales levels.
 
 #### Interpretation
 
 - **Lower MAPE**: Better model fit to training data
 - **Lower SMAPE**: Better symmetric model fit
 - **Higher R-squared**: Better explanatory power
+- **Lower CRPS**: Less uncertainty
 - **Comparison with holdout**: Helps identify overfitting (much better in-sample than holdout performance)
 
 ## Cross-Validated Holdout Accuracy Test
@@ -94,6 +98,8 @@ integers)
 - **MAPE**: Out-of-sample prediction accuracy
 - **SMAPE**: Out-of-sample symmetric prediction accuracy
 - **R-squared**: Out-of-sample explanatory power
+- **CRPS**: Out-of-sample predictive distribution accuracy
+
 
 ### Interpretation
 
@@ -215,8 +221,8 @@ modify the thresholds in `mmm_eval/metrics/threshold_constants.py`.
 
 ### Good Model Indicators
 
-- **Holdout Accuracy**: MAPE < 15%, SMAPE < 15%, R-squared > 0.8
-- **In-Sample Accuracy**: MAPE < 10%, SMAPE < 10%, R-squared > 0.9
+- **Holdout Accuracy**: MAPE < 15%, SMAPE < 15%, R-squared > 0.8, CRPS < 15%
+- **In-Sample Accuracy**: MAPE < 10%, SMAPE < 10%, R-squared > 0.9, CRPS < 10%
 - **Cross-Validation**: Out-of-sample MAPE/SMAPE similar to in-sample
 - **Refresh Stability**: Parameter changes < 10%
 - **Perturbation**: ROI changes < 5%
@@ -225,6 +231,7 @@ modify the thresholds in `mmm_eval/metrics/threshold_constants.py`.
 ### Warning Signs
 
 - **Poor Performance**: High MAPE/SMAPE or low R-squared
+- **High Uncertainty**: CRPS significantly higher than MAPE
 - **Overfitting**: Much better in-sample than holdout performance
 - **Unstable Model**: Large parameter changes
 - **Data Issues**: Missing values or extreme outliers

mmm_eval/adapters/base.py

@@ -86,7 +86,7 @@ def fit(self, data: pd.DataFrame) -> None:
         pass
 
     @abstractmethod
-    def predict(self, data: pd.DataFrame | None = None) -> np.ndarray:
+    def predict(self, data: pd.DataFrame | None = None) -> tuple[np.ndarray, np.ndarray]:
         """Make predictions on new data.
 
         Args:
@@ -97,34 +97,34 @@ def predict(self, data: pd.DataFrame | None = None) -> np.ndarray:
                   the fitted model state instead
 
         Returns:
-            Predicted values
+            Predicted values, tuple of (posterior_mean, posterior_distribution)
 
         """
         pass
 
     @abstractmethod
-    def fit_and_predict(self, train: pd.DataFrame, test: pd.DataFrame) -> np.ndarray:
+    def fit_and_predict(self, train: pd.DataFrame, test: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
         """Fit on training data and make predictions on test data.
 
         Args:
             train: dataset to train model on
             test: dataset to make predictions using
 
         Returns:
-            Predicted values.
+            Predicted values, tuple of (posterior_mean, posterior_distribution)
 
         """
         pass
 
     @abstractmethod
-    def fit_and_predict_in_sample(self, data: pd.DataFrame) -> np.ndarray:
+    def fit_and_predict_in_sample(self, data: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
         """Fit the model on data and return predictions for the same data.
 
         Args:
             data: dataset to train model on and make predictions for
 
         Returns:
-            Predicted values for the training data.
+            Predicted values for the training data, tuple of (posterior_mean, posterior_distribution)
 
         """
         pass

mmm_eval/adapters/meridian.py

@@ -519,11 +519,11 @@ def fit(self, data: pd.DataFrame, max_train_date: pd.Timestamp | None = None) ->
         self.analyzer = Analyzer(self.model)
         self.is_fitted = True
 
-    def _predict_on_all_data(self) -> np.ndarray:
+    def _predict_on_all_data(self) -> tuple[np.ndarray, np.ndarray]:
         """Make predictions on all data provided to fit().
 
         Returns
-            predicted values on data provided to fit().
+            predicted values on data provided to fit(), tuple of (posterior_mean, posterior_distribution)
 
         """
         if not self.is_fitted or self.analyzer is None:
@@ -532,9 +532,11 @@ def _predict_on_all_data(self) -> np.ndarray:
         # shape (n_chains, n_draws, n_times)
         preds_tensor = self.analyzer.expected_outcome(aggregate_geos=True, aggregate_times=False, use_kpi=True)
         posterior_mean = np.mean(preds_tensor, axis=(0, 1))
-        return posterior_mean
+        # A (n sample x n_chain) by n_time matrix
+        posterior_distribution = np.array(preds_tensor).reshape(-1, preds_tensor.shape[-1])
+        return posterior_mean, posterior_distribution
 
-    def predict(self, data: pd.DataFrame | None = None) -> np.ndarray:
+    def predict(self, data: pd.DataFrame | None = None) -> tuple[np.ndarray, np.ndarray]:
         """Make predictions using the fitted model.
 
         This returns predictions for the entirety of the dataset passed to fit() unless
@@ -548,22 +550,22 @@ def predict(self, data: pd.DataFrame | None = None) -> np.ndarray:
             data: Ignored - Meridian uses the fitted model state for predictions.
 
         Returns:
-            Predicted values
+            Predicted values, tuple of (posterior_mean, posterior_distribution)
 
         Raises:
             RuntimeError: If model is not fitted
 
         """
-        posterior_mean = self._predict_on_all_data()
+        posterior_mean, posterior_distribution = self._predict_on_all_data()
 
         # if holdout mask is provided, use it to mask the predictions to restrict only to the
         # holdout period
         if self.holdout_mask is not None:
             posterior_mean = posterior_mean[self.holdout_mask]
+            posterior_distribution = posterior_distribution[:, self.holdout_mask]
+        return posterior_mean, posterior_distribution
 
-        return posterior_mean
-
-    def fit_and_predict(self, train: pd.DataFrame, test: pd.DataFrame) -> np.ndarray:
+    def fit_and_predict(self, train: pd.DataFrame, test: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
         """Fit the Meridian model and make predictions given new input data.
 
         The full dataset must be passed to `fit()`, since making out-of-sample predictions
@@ -575,27 +577,28 @@ def fit_and_predict(self, train: pd.DataFrame, test: pd.DataFrame) -> np.ndarray
             test: Test data
 
         Returns:
-            Predicted values for the test period
+            Predicted values for the test period, tuple of (posterior_mean, posterior_distribution)
 
         """
         train_and_test = pd.concat([train, test])
         max_train_date = train[self.date_column].squeeze().max()
         self.fit(train_and_test, max_train_date=max_train_date)
         return self.predict()
 
-    def fit_and_predict_in_sample(self, data: pd.DataFrame) -> np.ndarray:
+    def fit_and_predict_in_sample(self, data: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
         """Fit the model on data and return predictions for the same data.
 
         Args:
             data: dataset to train model on and make predictions for
 
         Returns:
-            Predicted values for the training data.
+            Predicted values for the training data, tuple of (posterior_mean, posterior_distribution)
 
         """
         # no max train date specified, so predictions are all in-sample
         self.fit(data)
-        return self._predict_on_all_data()
+        posterior_mean, posterior_distribution = self._predict_on_all_data()
+        return posterior_mean, posterior_distribution
 
     def get_channel_roi(
         self,

mmm_eval/adapters/pymc.py

@@ -198,15 +198,15 @@ def fit(self, data: pd.DataFrame) -> None:
         self._channel_roi_df = self._compute_channel_contributions(data)
         self.is_fitted = True
 
-    def predict(self, data: pd.DataFrame | None = None) -> np.ndarray:
+    def predict(self, data: pd.DataFrame | None = None) -> tuple[np.ndarray, np.ndarray]:
         """Predict the response variable for new data.
 
         Args:
             data: Input data for prediction. This parameter is required for PyMC
                 predictions and cannot be None.
 
         Returns:
-            Predicted values
+            Predicted values, tuple of (posterior_mean, posterior_distribution)
 
         Raises:
             RuntimeError: If model is not fitted
@@ -221,39 +221,42 @@ def predict(self, data: pd.DataFrame | None = None) -> np.ndarray:
 
         if InputDataframeConstants.RESPONSE_COL in data.columns:
             data = data.drop(columns=[InputDataframeConstants.RESPONSE_COL])
-        predictions = predictions = self.model.predict(
-            data, extend_idata=False, include_last_observations=True, **self.predict_kwargs
+
+        posterior_prediction = self.model.predict_posterior(
+            data, extend_idata=False, include_last_observations=True, combined=True, **self.predict_kwargs
         )
-        return predictions
+        posterior_prediction = np.array(posterior_prediction)
+        posterior_mean = posterior_prediction.mean(axis=1)
+        return posterior_mean, posterior_prediction.transpose()
 
-    def fit_and_predict(self, train: pd.DataFrame, test: pd.DataFrame) -> np.ndarray:
+    def fit_and_predict(self, train: pd.DataFrame, test: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
         """Fit on training data and make predictions on test data.
 
         Arguments:
             train: training dataset
             test: test dataset
 
         Returns:
-            model predictions.
+            model predictions, tuple of (posterior_mean, posterior_distribution)
 
         """
         self.fit(train)
         return self.predict(test)
 
-    def fit_and_predict_in_sample(self, data: pd.DataFrame) -> np.ndarray:
+    def fit_and_predict_in_sample(self, data: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]:
         """Fit the model on data and return predictions for the same data.
 
         Args:
             data: dataset to train model on and make predictions for
 
         Returns:
-            Predicted values for the training data.
+            Predicted values for the training data, tuple of (posterior_mean, posterior_distribution)
 
         """
         self.fit(data)
         if self.model is None:
             raise RuntimeError("Model must be fit before prediction.")
-        return self.model.predict(data, extend_idata=False, **self.predict_kwargs)
+        return self.predict(data)
 
     def get_channel_roi(
         self,