### Is your feature request related to a problem? Please describe.
For probabil…istic predictors that predict a (univariate) **distribution** $D$ rather than give a point estimate, the continuous ranked probability score between $D$ and $y$ is
```math
CRPS(D, y) = \int_{{supp}(D)} (F_D(x) - 1[x \leq y])^2 dx.
```
See here:
https://en.wikipedia.org/wiki/Scoring_rule#Continuous_ranked_probability_score
It is a proper scoring rule for probabilistic inference, and useful in time series inference. When $D$ is Gaussian, it has a closed-form formula devised in [1]. I would like to see it implemented, since Gaussian predictors are ubiquitous in statistics (i.e. the entire confidence interval theory is based on Gaussian sampling distributions).
[1]: Gneiting, T., Raftery, A.E., Westveld III, A.H. and Goldman, T., 2005. Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation. Monthly weather review, 133(5), pp.1098-1118.
### Describe the solution you'd like.
Here is my implementation. I hope the adaptation to the SciPy coding standards is easy:
```python
from scipy.stats import norm
import numpy as np
def crps_gaussian(mu, sigma, y):
"""
Compute the Continuous Ranked Probability Score (CRPS) for a Gaussian forecast.
Gneiting, T., Raftery, A.E., Westveld III, A.H. and Goldman, T., 2005. Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation. Monthly weather review, 133(5), pp.1098-1118.
Parameters
----------
mu : float or array-like
Mean of the Gaussian forecast.
sigma : float or array-like
Standard deviation of the Gaussian forecast (must be positive).
y : float or array-like
Observed value(s).
Returns
-------
crps : float or numpy.ndarray
The CRPS value(s).
"""
mu = np.asarray(mu)
sigma = np.asarray(sigma)
y = np.asarray(y)
z = (y - mu) / sigma
phi = norm.pdf(z)
Phi = norm.cdf(z)
# Closed-form expression for CRPS
crps = sigma * (z * (2 * Phi - 1) + 2 * phi - 1/np.sqrt(np.pi))
return crps
```
### Describe alternatives you've considered.
Computing the integral for generic distributions using Monte-Carlo simulation, using the property that
```math
CRPS(D, y) = E _{x\sim D} [\lvert x - y \rvert] - \frac{1}{2} E_{x,x'\sim D} [|x - x'|]
```
In this sense, it generalizes the absolute distance for probabilistic forecasts.
However, the closed-form formula for Gaussians makes it extremely efficient for this special case. We don't need Monte-Carlo simulation.
### Additional context (e.g. screenshots, GIFs)
CRPS between the Gaussian distributions $\mathcal{N}(0.8, 0.01)$, $\mathcal{N}(0.8, 0.05)$, and $\mathcal{N}(0.8, 0.1$ between various values of $y$ between 0 and 1:
```python
mu = 0.8
ys = np.linspace(0, 1, 1000)
plt.plot(ys, np.abs(mu - ys), label='Absolute error')
plt.plot(ys, crps_gaussian(mu, 0.01, ys), label='CRPS ($\sigma=0.01$)')
plt.plot(ys, crps_gaussian(mu, 0.05, ys), label='CRPS ($\sigma=0.05$)')
plt.plot(ys, crps_gaussian(mu, 0.1, ys), label='CRPS ($\sigma=0.1$)')
plt.axhline(y=0, color='k', linestyle='--')
plt.xlim([0.6, 1])
plt.ylim([-0.1, 0.3])
plt.legend()
plt.show()
```
