5.9 予測の分布の正確性評価

ここまでの指標は全て点予測の正確性を測るものでした。予測の分布を評価するには、他の指標を使う必要があります。

分位点スコア

前節のGoogle株価事例を考えてみましょう。図5.23は、ナイーブ法による予測の80%区間予測を示しています。

google_fc %>%
  filter(.model == "ナイーブ") %>%
  autoplot(bind_rows(google_2015, google_jan_2016), level=80)+
  labs(y = "米ドル", level = "区間予測",
       title = "Google株価終値")
Google株価の2016年1月のナイーブ予測と80%区間予測

図 5.23: Google株価の2016年1月のナイーブ予測と80%区間予測

この区間予測の下限は、予測の分布の10番目のパーセンタイル(言い換えれば、0.1の分位点)に当たります。ですから、実績値が下限を下回る確率は約10%で、上回る確率は90%と期待しているわけです。実績値をこのパーセンタイルと比較するとしたら、下回るより上回ることの方が多いとみておく必要があるわけです。

より一般的に、将来の時点\(t\)での確率\(p\)の分位点予測に興味があるとして、それを\(f_{p,t}\)と表記しましょう。つまり、観測値\(y_t\)\(f_{p,t}\)を下回る確率は\(p\)だけと期待しているわけです。例えば、10番目のパーセンタイルは\(f_{0.10,t}\)の表記になります。もし、\(t\)時点での観測値を\(y_{t}\)とすると、分位点スコアは、以下のように表せます。 \[ Q_{p,t} = \begin{cases} 2(1 - p) \big(f_{p,t} - y_{t}\big), & \text{if $y_{t} < f_{p,t}$}\\ 2p \big(y_{t} - f_{p,t}\big), & \text{if $y_{t} \ge f_{p,t}$} \end{cases} \] これは「ピンボール損失関数」と呼ばれることもあります。そのグラフがピンボール・テーブル上のボールの軌道に似ているからです。乗数2は省略されることがしばしばありますが、残しておくと解釈が若干容易になります。\(Q_{p,t}\)の値が小さいほど、分位点推計は上手くなされていることを表します。

分位点スコアは絶対誤差のように解釈することができます。実際、\(p=0.5\)の場合、分位点スコア\(Q_{0.5,t}\)は絶対誤差と同値です。\(p\)がそれ以外の値の場合、「誤差」\((y_t - f_{p,t})\)を条件が成立しなかった方の確率をウェートとして加重したものになります。もし\(p>0.5\)なら、観測値が分位点推計を上回る場合、下回る場合よりも大きなペナルティが課されます。もし\(p<0.5\)なら、逆になります。

5.23では、1期先(2016年1月4日)の10%分位点予測は\(f_{0.1,t} = 744.54\)で、観測値は\(y_t = 741.84\)でしたから、分位点スコアは以下のようになります。 \[ Q_{0.1,t} = 2(1-0.1) (744.54 - 741.84) = 4.86 \] 分位点スコアは、accuracy()の中でquantile_score()関数を以下のように指定すれば、簡単に計算できます。

google_fc %>%
  filter(.model == "ナイーブ", Date == "2016-01-04") %>%
  accuracy(google_stock, list(qs=quantile_score), probs=0.10)
#> # A tibble: 1 × 4
#>   .model   Symbol .type    qs
#>   <chr>    <chr>  <chr> <dbl>
#> 1 ナイーブ GOOG   Test   4.86

Winklerスコア

分位点を2、3評価するよりも、区間予測を評価することに興味があることがしばしばあります。 Winkler (1972) が提案したWinklerスコアは、この目的のために設計されています。もし、\(t\)時点の\(100(1-\alpha)\)%区間予測が\([\ell_{\alpha,t}, u_{\alpha,t}]\)なら、Winklerスコアは、区間予測の長さに、観測値が区間予測外に落ちたときのペナルティを加味して、次のように定義されます。 \[ W_{\alpha,t} = \begin{cases} (u_{\alpha,t} - \ell_{\alpha,t}) + \frac{2}{\alpha} (\ell_{\alpha,t} - y_t) & \text{if } y_t < \ell_{\alpha,t} \\ (u_{\alpha,t} - \ell_{\alpha,t}) & \text{if } \ell_{\alpha,t} \le y_t \le u_{\alpha,t} \\ (u_{\alpha,t} - \ell_{\alpha,t}) + \frac{2}{\alpha} (y_t - u_{\alpha,t}) & \text{if } y_t > u_{\alpha,t} \end{cases} \] 観測値が区間予測内に入っていれば、Winklerスコアは単に区間予測の長さです。ですから、区間予測が短ければ、スコアは小さくなります。しかし、観測値が区間予測外に落ちてしまうと、どれだけ外側にはずれたかに比例したペナルティが課されます。

区間予測は通常、分位点を使って下限\(\ell_{\alpha,t} = f_{\alpha/2,t}\)と上限\(u_{\alpha,t} = f_{1-\alpha/2,t}\)を構築します。下限と上限の分位点スコアを合計して、\(\alpha\)で割ったものがWinklerスコアです。 \[ W_{\alpha,t} = (Q_{\alpha/2,t} + Q_{1-\alpha/2,t})/\alpha \]

5.23に示した2016年1月4日についての1期先予測の80%区間予測は[744.54, 773.22]で、実績値は741.84でしたので、Winklerスコアは、以下のようになります。 \[ W_{\alpha,t} = (773.22 - 744.54) + \frac{2}{0.2} (744.54 - 741.84) = 55.68 \] Winklerスコアは、accuracy()の中でwinkler_score()関数を以下のように指定すれば、簡単に計算できます。

google_fc %>%
  filter(.model == "ナイーブ", Date == "2016-01-04") %>%
  accuracy(google_stock,
    list(winkler = winkler_score), level = 80)
#> # A tibble: 1 × 4
#>   .model   Symbol .type winkler
#>   <chr>    <chr>  <chr>   <dbl>
#> 1 ナイーブ GOOG   Test     55.7

連続階数化確率スコア (CRPS)

特定の分位点や区間予測ではなく、予測の分布全体に興味あることがしばしばあります。その場合、全ての\(p\)値にわたって分位点スコアを計算して、その平均を取ることで、連続階数化確率スコア Continuous Ranked Probability Score、もしくは、CRPS (Gneiting & Katzfuss, 2014)が得られます。

Google株価の例では、テストセット中の全日についての平均CRPS値を計算しています。CRPS値は予測の分布全体から計算された加重絶対誤差のようなもので、加重には確率を考慮しています。

google_fc %>%
  accuracy(google_stock, list(crps = CRPS)) %>% 
  arrange(crps)
#> # A tibble: 3 × 4
#>   .model   Symbol .type  crps
#>   <chr>    <chr>  <chr> <dbl>
#> 1 ナイーブ GOOG   Test   26.5
#> 2 ドリフト GOOG   Test   33.5
#> 3 平均     GOOG   Test   76.7

この例では、ナイーブ法がドリフト法や平均法よりも、予測の分布については良い予測をしているようです。

スキルスコアを用いた尺度から独立した比較

点予測のときと同様、異なる尺度の系列について、複数の手法から得た予測の分布の正確性を比較できると有益です。その目的のため、点予測では尺度調整済み誤差を用いました。別のアプローチは、スキルスコアを用いることです。スキルスコアは点予測の正確性にも、予測の分布の正確性にも用いることができます。

スキルスコアでは、あるベンチマークとなる手法に対する予測正確性指標を計算します。例えば、ナイーブ法をベンチマークとして、ドリフト法でも予測するならば、ナイーブ法と比較したドリフト法のCRPSスキルスコアは、以下のように計算できます。 \[ \frac{\text{CRPS}_{\text{ナイーブ}} - \text{CRPS}_{\text{ドリフト}}}{\text{CRPS}_{\text{ナイーブ}}} \] これは、ナイーブ法に対してドリフト法がCRPSをどれだけの比率で改善できたかを表しています。accuracy()関数を使えば、簡単に計算できます。

google_fc %>%
  accuracy(google_stock, list(skill = skill_score(CRPS))) %>% 
  arrange(desc(skill))
#> # A tibble: 3 × 4
#>   .model   Symbol .type  skill
#>   <chr>    <chr>  <chr>  <dbl>
#> 1 ナイーブ GOOG   Test   0    
#> 2 ドリフト GOOG   Test  -0.266
#> 3 平均     GOOG   Test  -1.90

ナイーブ法のスキルスコアがゼロなのは、自分自身に対して改善などできないので当然です。他の2つの手法のCRPS値はナイーブ法と比べて大きくなってしまったので、スキルスコアは負値になっています。ドリフト法はナイーブ法と比べ26.6%悪化しています。

skill_score()関数は、ベンチマーク予測がfableオブジェクトに含まれていなくても、CRPSを計算して常に適切なベンチマーク予測と比較します。データに季節性があるなら、使用されるベンチマークはナイーブ法ではなく季節ナイーブ法になります。ベンチマーク予測にも確実に同じ訓練データが使われるよう、accuracy()関数に渡すデータの始点は訓練データの始点と同じにすることが大事です。

skill_score()関数はどんな正確性指標にも使えます。例えば、skill_score(MSE)とすれば、さまざまな系列間でMSE値を比較できます。しかし、テストセットが十分大きくないと、正確性指標、特に分母の正確性指標、が信頼できません。そのため、尺度から独立した点予測の正確性指標としては、MASEやRMSSEの方が多くの場合望ましいでしょう。

参考文献

Gneiting, T., & Katzfuss, M. (2014). Probabilistic forecasting. Annual Review of Statistics and Its Application, 1(1), 125–151. [DOI]
Winkler, R. L. (1972). A decision-theoretic approach to interval estimation. Journal of the American Statistical Association, 67(337), 187–191. [DOI]