Calculate the Brier score (BS)

    The Brier score (BS) scores binary forecasts $k \in \{0, 1\}$,

    ..math:
        BS(p, k) = (p_1 - k)^2,

    where $p_1$ is the forecast probability of $k=1$.

    Parameters
    ----------
    observations, forecasts : array_like
        Broadcast compatible arrays of forecasts (probabilities between 0 and
        1) and observations (0, 1 or NaN).

    Returns
    -------
    out : np.ndarray
        Brier score for each forecast/observation.

    References
    ----------
    Jochen Broecker. Chapter 7 in Forecast Verification: A Practitioner's Guide
        in Atmospheric Science. John Wiley & Sons, Ltd, Chichester, UK, 2nd
        edition, 2012.
        https://drive.google.com/a/climate.com/file/d/0B8AfRcot4nsIYmc3alpTeTZpLWc
    Tilmann Gneiting and Adrian E. Raftery. Strictly proper scoring rules,
        prediction, and estimation, 2005. University of Washington Department of
        Statistics Technical Report no. 463R.
        https://www.stat.washington.edu/research/reports/2004/tr463R.pdf

def brier_score(observations, forecasts):
    """
    Calculate the Brier score (BS)

    The Brier score (BS) scores binary forecasts $k \in \{0, 1\}$,

    ..math:
        BS(p, k) = (p_1 - k)^2,

    where $p_1$ is the forecast probability of $k=1$.

    Parameters
    ----------
    observations, forecasts : array_like
        Broadcast compatible arrays of forecasts (probabilities between 0 and
        1) and observations (0, 1 or NaN).

    Returns
    -------
    out : np.ndarray
        Brier score for each forecast/observation.

    References
    ----------
    Jochen Broecker. Chapter 7 in Forecast Verification: A Practitioner's Guide
        in Atmospheric Science. John Wiley & Sons, Ltd, Chichester, UK, 2nd
        edition, 2012.
        https://drive.google.com/a/climate.com/file/d/0B8AfRcot4nsIYmc3alpTeTZpLWc
    Tilmann Gneiting and Adrian E. Raftery. Strictly proper scoring rules,
        prediction, and estimation, 2005. University of Washington Department of
        Statistics Technical Report no. 463R.
        https://www.stat.washington.edu/research/reports/2004/tr463R.pdf
    """
    machine_eps = np.finfo(float).eps
    forecasts = np.asarray(forecasts)
    if (forecasts < 0.0).any() or (forecasts > (1.0 + machine_eps)).any():
        raise ValueError('forecasts must not be outside of the unit interval '
                         '[0, 1]')
    observations = np.asarray(observations)
    if observations.ndim > 0:
        valid_obs = observations[~np.isnan(observations)]
    else:
        valid_obs = observations if not np.isnan(observations) else []
    if not set(np.unique(valid_obs)) <= {0, 1}:
        raise ValueError('observations can only contain 0, 1, or NaN')
    return (forecasts - observations) ** 2