function [CRPS,mCRPS] = crps(frcst,y)
%CRPS  Continuous ranked probability score for ensemble forecasts
%   [crps,mcrps] = CRPS(frcst,y), where frcst is the ensemble forecast and 
%     y are the observations that the ensemble is verified against. 
%     'frcst' has dimensions [ens_members,time] and 'y' is [time] and 
%     note that the time dimension can be 1.
%
%   Calculates the continuous ranked probability score, assuming that the
%   ensemble is normally distributed. Note that the score is negatively
%   oriented, so that smaller values are better. The score is also 'proper'
%   meaning that it does not reward hedging (or under-forecasting) to get a 
%   better skill score.
%
%   Formula taken from Equation (5):
%
%   Calibrated Probablistic Forecasting Using Ensemble Model Output
%   Statistics and Minimum CRPS Estimation. Gneiting, Raftery,
%   Westveld, Goldman. Monthly Weather Review 2004
%
%   http://journals.ametsoc.org/doi/pdf/10.1175/MWR2904.1 
%
%   If you want to convert CRPS into a skill score (defined such that 0
%   indicates no improvement over a baseline and 1 indicates a perfect
%   forecast), then compute both a baseline and forecast and then: 
%   1 - forecast_score/baseline_score
%
%   Nathan Steiger, December 2015

% Parameters of the forecast
mu = mean(frcst); mu=mu(:);
sig = std(frcst); sig=sig(:);
% Standardized y
sy = (y(:) - mu)./sig;
% Precomputations
npdf = normpdf(sy);
ncdf = normcdf(sy);
pi_inv = 1./sqrt(pi);
% Compute CRPS for each time
CRPS = sig.*(sy.*(2.*ncdf - 1) + 2.*npdf - pi_inv);
% Mean CRPS over time
mCRPS=mean(CRPS);

end

