はじめに

こんにちは@KokiUenoです．

本記事は，予測分布の不確実性と精度を定量的に評価する指標でよく用いられるCRPS　（Continuous Ranked Probability Score）について、特に連続物理量$x$に対する調査結果をまとめます

この記事の対象者

CRPSの概要を学びたい
CRPSの無限大計算を実際にはどのように計算しているか知りたい

本文

ここでは，CRPSの概要と実際の計算についてざっくりと紹介していきます．

1. CRPSの概要

1. CRPSとは

CRPS （Continuous Ranked Probability Score）は，確率予測の統計検証指標の一つであり，予測分布の不確実性と精度を定量的に評価する指標として使用されています．
具体的には，CRPS は『観測値$x$に対する累積分布関数$A(x)$』と『予測された累積分布$F(x)$』を比較評価する指標です．

2. CRPSの一般式

こちらのAppendixを参考にして，連続物理量$x$に対するCRPSは次式で定義されます．
$$
\text{CRPS} = \int_{-\infty}^{\infty} \left[ F_i(x) - A_i(x) \right]^2 dx
$$
ここで，Nは標本数，$F_i(x)$と$A_i(x)$はそれぞれ予測と観測の累積分布関数（CDF: Cumulative Distribution Function）です．

3. CRPS-sum

最終的には，N個の予測や観測がある場合，以下のようにCRPSの平均値（CRPS-sum）をとります．
$$
\text{CRPS-sum} = \frac{1}{N} \sum_{i=1}^{N} \int_{-\infty}^{\infty} \left[ F_i(x) - A_i(x) \right]^2 dx
$$

4. CRPSの可視化

ここでは，CRPSを可視化してイメージを掴んでみます．

例として，ある工場で出荷されたボルトの寿命についてCRPSで評価することを考えましょう．ボルトの経年劣化による寿命が5年以上である確率はほぼ100%である状況を考えます．この場合，観測されたボルトの寿命の累積分布関数$A(x)$は，寿命が5年未満のときには0に近く，5年を超えると急激に上昇し，最終的には1に達します．これは，「ボルトの寿命が5年以上である確率はほぼ100%である」という状況を反映しています．

よって，観測分布$A(x)$は以下のように表せます．
$$
A(x) = \begin{cases}
0, & \text{if } x < 5 \newline
1, & \text{if } x \geq 5
\end{cases}
$$

次に，ボルトの寿命を予測し，その予測分布$F(x)$は正規分布に従うとすると，以下のように表せます．
$$
F(x) \sim \text{Logistic}(\mu = 5, \sigma = \frac{1}{k})
$$

$$
F(x) = \frac{1}{1 + e^{-k(x - 5)}}
$$

ここで、$\mu$ は平均、$\sigma$ は分散を表します．
$k$はスケールを表し，今回は$k=1$とします．

この時，CRPSは以下のように表現することができます．[（参考）]

まとめとして，図の薄ピンクの範囲がCRPSであり，『この面積を小さくする＝積分結果を小さくする』ことを目標にすることで，予測分布$P(x)$を実際の観測値の分布$A(x)$に近づけることができます．そのため，CRPSは予測の精度を評価するための指標の一つとして使用されています．不確実性や精度の評価として使用される理由としては，定量的な評価が可能なためです．

2. 計算

1. CSDI論文を元に近似計算を考える

上記の連続値$x$における積分範囲が${-\infty}$から${\infty}$までであるため、累積分布関数全体にわたって評価を行っていることがわかります．
しかし，${\infty}$計算を実際に行うことはできないため，計算時には何かしらの近似計算を行なっているはずです．

そこで今回は，NeurIPS2021で採択されたCSDIをベースに考えていきます．
※実際には，時系列予測のCRPS-sum計算時には，特徴量数や系列数における正規化処理などもありますが，今回は簡単のため以下のように考えます．CRPS単体の計算方法に変わりはありません．

$$
\text{CRPS-sum}(F^{-1}, x) = \frac{1}{N} \sum_{i=1}^{N} \int_{0}^{1} 2\Lambda_\alpha(F_i^{-1}(\alpha), x) d\alpha
$$

（CSDIのpaperのAppendix【E.3】式(16)ではN=1のCRPSが記載されています．例えば，electricityデータセットのCRPS-sumには370系列の特徴量があるので，N=370の平均値が0.017(0.000)であることを示しています．）

この式は単一の予測だけではなく，予測分布全体の精度を評価するために用いられます．予測分布が観測値に対してどれだけ近いかを分布の各ポイントに対して評価し，それを平均する形でCRPSを算出しています．具体的には，分布の逆累積分布関数$F^{-1}$（quantile function）を使用して，観測値との誤差を累積して評価しています．

2. 逆累積分布関数

逆累積分布関数とは？
逆累積分布関数（ICDF; Inverse Cumulative Distribution Function ） = 分位関数（Quantile Function）とは何か，先ほどのボルトの例を用いて説明します．
累積分布関数CDFは$F(x)$の$x$へ時間を入れることで，ボルトが寿命である（故障する）確率を得ることができました．
一方で，逆累積分布関数ICDFでは$F^{-1}(α)$の$α$へ全体の何%かを入れることで，どのくらいの時間でボルトが故障するかを知ることができます．
なぜ逆累積分布関数を使用するのか？
ここに近似計算の全てが詰まっています．CDFのまま計算するには，$-\infty \sim \infty$の範囲の積分が必要となりますが，ICDFを使用することで，0から1の有限範囲の計算に収めることができます．

3. ヒンジロス計算

次に，ヒンジロス$\Lambda$計算の部分を見ていきます．ヒンジロス部分は以下のように展開でき，分位$α$に基づいて，予測された値$F^{−1}(α)$と，観測された値$x$の間の距離を計算しています．具体的には，予測が観測値よりも大きい場合と小さい場合で異なる重みを与えています．（論文内では0.05）

$$
\Lambda_\alpha(F_i^{-1}(\alpha), x) =
\begin{cases}
\alpha \cdot (F_i^{-1}(\alpha) - x), & \text{if } F_i^{-1}(\alpha) > x \newline
(1-\alpha) \cdot (x - F_i^{-1}(\alpha)), & \text{if } F_i^{-1}(\alpha) \leq x
\end{cases}
$$

4. 連続値から離散値へ

実際の$\text{CRPS}$計算では以下の式で連続値の処理を離散値計算として近似しています．
以下の近似式で使用されている$\sum$は$\text{CRPS-sum}$で平均するための$\sum$ではなく，0 ≦ $\alpha$ ≦ 1の連続範囲で積分する計算を，19個の離散値として計算して平均を取っています．これにより，$\alpha=5 \sim 95$ %の5%ずつ分けた四分位数それぞれの範囲で計算して平均を取ることが可能になります．$\Lambda$の前にある2は，分位点を均等に配置する際に(1/2n)を行うため，計算を簡単にするため正数倍の計算をするために行うようです．
$$
\text{CRPS}(F^{-1}, x) \approx \sum_{j=1}^{19} \frac{2\Lambda_{j \times 0.05}(F^{-1}(j \times 0.05), x)}{19}
$$

おわりに

今回は，不確実性の定量的な評価指標として使用されているCRPS　（Continuous Ranked Probability Score）についての概要と，実際の計算について調査結果を公開しました．
今後は不確実性についての記事を公開していく予定です．

ここまで読んでいただけた方はいいねとストックよろしくお願いします．
@KokiUeno をフォローいただけるととてもうれしく思います．

また明日の記事でお会いしましょう！