【気象予報】予報精度の評価指標を計算してみる

# 1. 概要
気象予報がどれくらい当たったかを検証するための指標がいくつかあります．

降水の例で言えば，

- 適中率
- 空振り率
- 見逃し率

などが挙げられます．

本記事では，超簡易的な機械学習モデルで予測した「最高気温」「最低気温」「降水確率」「降水の有無」の予報精度をさまざまな指標で検証してみます．

# 2. データと機械学習モデル
超簡易的なデータを使って，超簡易的な機械学習モデルで学習させました．

今回使用したデータとモデルについて概要を説明します．

## データ
気象庁の「[過去の気象データダウンロード](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php)」から取得したデータを使用しました．

> - 地点：京都
> - 期間：2022年11月26日～2023年11月26日
> - データの種類：時別値

観測項目は，

- 気温
- 降水量
- 日照時間
- 風向・風速
- 現地気圧
- 相対湿度
- 蒸気圧
- 露点温度
- 天気

を用いました．

風向・風速の情報から，風速の東西成分と南北成分を算出しました．
天気はone-hot Encodingで変換しました．

また，1時間～9時間前の各観測項目を特徴量に追加しました．
いわゆる「ラグ特徴量」ってやつです．

全部で特徴量は151個です．

## 機械学習モデル

### 気温予測
以下を入出力とする重回帰モデルです．

- 入力：時別データ（特徴量151個）
- 出力：24時間後の気温

最高気温や最低気温は単純に，各日の気温の最大値・最小値と定義します．

### 降水予測
以下を入出力とするロジスティック回帰モデルです．

- 入力：時別データ（特徴量151個）
- 出力：翌日の日降水量が0.5mm以上かどうか（確率）

各日の降水確率は，24個の出力値の平均としました．
分類のしきい値は0.5としています．

# 3. 予報精度の評価

全データの80%を訓練用データ，残り20%を評価用データとしました．
（もちろん，データリークしないよう，評価データは訓練データよりも後の時間のデータとしています．）

下図は評価用データに対する予測結果です．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/bf24411c-17d9-89cc-3cae-58a4d25920a3.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/82099d58-ea08-5d1b-efc2-23960750fd87.png)

これらの予測に対して，精度を評価してみましょう．

## 気温予測の評価 
### 平均誤差（ME; Mean Error）

平均誤差は，予測値を$\hat{y_i}$，実況値を$y_i$とすると

```math
\text{ME}(\boldsymbol{\hat{y}},\boldsymbol{y}) = \frac{1}{N}\sum_{i=1}^{N}(\hat{y_i}-y_i)
```

と定義されます．

予報値が上にブレやすいか，下にブレやすいかがわかりやすい指標です．

ただし正の値と負の値が打ち消しあって$\text{ME}=0$となり，実際には誤差が生じているのにあたかも誤差がないように見えてしまう，という問題点があります．

今回の予測モデルでの平均誤差は，以下のようになりました．

```
Highest Temperature
- ME   : -0.2472207344355263
Lowest Temperature
- ME   : 0.6470339790754622
```

最高気温の予測は，やや下にブレやすく，一方で最低気温はやや上にブレやすいことがわかります．


### 2乗平均平方根誤差（RMSE; Root Mean Square Error）

2乗平均平方根誤差は，予測値を$\hat{y_i}$，実況値を$y_i$とすると

```math
\text{RMSE}(\boldsymbol{\hat{y}},\boldsymbol{y}) = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(\hat{y_i}-y_i)^2}
```

と定義されます．機械学習をやったことがある人ならお馴染みですね．

平均誤差（ME）とは違い，二乗の計算が入るので，どれくらい予報と実況にズレが生じたのかがわかります．

$|\hat{y_i}-y_i|<1$の場合，2乗すると値が小さくなるので微小ノイズに対して鈍感です．
逆に$|\hat{y_i}-y_i|>1$の場合，2乗すると値が大きくなるので大きな誤差に対して敏感です．

今回の予測モデルでの2乗平均平方根誤差は，以下のようになりました．

```
Highest Temperature
- RMSE : 2.1101543396002227
Lowest Temperature
- RMSE : 2.0555428536482445
```

今回の線形モデルだと，最高気温も最低気温もだいたい2℃くらいズレるといったところでしょうか．


## 降水予測の評価

### 適中率
適中率といってもいくつか種類があります．

以下は「[気象庁｜検証方法の説明](https://www.data.jma.go.jp/fcd/yoho/kensho/explanation.html)」からの引用です．

- **降水の有無**の適中率
    - 総適中回数の全予報数に対する割合
- **「降水あり」** の適中率
    - 予報が｢降水あり｣だった場合だけを取り出して、そのうち実況が｢降水あり｣となった割合
- **「降水なし」** の適中率
    - 予報が｢降水なし｣だった場合だけを取り出して、そのうち実況が｢降水なし｣となった割合

予報と実況の適中/不適中をまとめた表（混同行列といいます）を使って計算してみます．

|  | 予報「降水なし」| 予報「降水あり」 |
|:-:|:-:|:-:|
| **実況「降水なし」**  |  $A$ | $B$  |
| **実況「降水あり」**  |  $C$ | $D$  |

**降水の有無の適中率**は，総適中回数の全予報数に対する割合なので

```math
\frac{A+D}{A+B+C+D}
```
です．機械学習の文脈では「正解率(accuracy)」と呼ばれます．

**「降水あり」の適中率**は，予報が｢降水あり｣だった場合だけを取り出して、そのうち実況が｢降水あり｣となった割合なので

```math
\frac{D}{B+D}
```
です．機械学習の文脈では「適合率(precision)」と呼ばれます．

**「降水なし」の適中率**は，予報が｢降水なし｣だった場合だけを取り出して、そのうち実況が｢降水なし｣となった割合なので

```math
\frac{A}{A+C}
```

です．機械学習の文脈では「特異度(specificity)」と呼ばれます．

今回の予測モデルでは，次のような混同行列になりました．

```
[[58  0]
 [15  1]]
```
これより，適中率は次のように求まります．
```
適中率（降水の有無）: 79.73%
「降水あり」予報の適中率: 100.00%
「降水なし」予報の適中率: 79.45%
```

### 見逃し率，空振り率，捕捉率
以下の定義も「[気象庁｜検証方法の説明](https://www.data.jma.go.jp/fcd/yoho/kensho/explanation.html)」からの引用して説明します．

**見逃し率**は，「見逃し」回数の全予報数に対する割合です．
「見逃し」とは、予報が｢降水なし｣だったのに実況が｢降水あり｣となった場合のことです．
混同行列と照らし合わせると，

```math
\frac{C}{A+B+C+D}
```
です．

**空振り率**は，「空振り」回数の全予報数に対する割合です．
「空振り」とは、予報が｢降水あり｣だったのに実況が｢降水なし｣となった場合のことです．
混同行列と照らし合わせると，

```math
\frac{B}{A+B+C+D}
```
です．

**捕捉率**は，実況が｢降水あり｣となった場合だけを取り出して、そのうち予報も｢降水あり｣
だった割合です．
混同行列と照らし合わせると，

```math
\frac{D}{C+D}
```
です．機械学習の文脈では「再現率(recall)」と呼ばれます．

今回の予測モデルでは，次のような結果となりました．
```
見逃し率: 20.27%
空振り率: 0.00%
捕捉率: 6.25%
```

### スレットスコア，バイアススコア
**スレットスコア**（TS; Threat Score）は，発生数が少ない現象（＝降水）の予測性能を評価するときによく使われます．

```math
\frac{D}{B+C+D}
```

「降水あり」の適中数を分子に持ってくるのはいいとして，分母には「降水なし」の適中数を**除いた**数を持ってきていることがポイントです．

**バイアススコア**（BS; Bias Score）は，「降水あり」と予報した回数と，実況が「降水あり」となった回数の比です．

```math
\frac{C+D}{B+D}
```

バイアススコアが1なら予報と実況で「降水あり」の回数が一致しています．1より大きいと「降水あり」という予報回数が多すぎるということになり，1より小さいと「降水あり」という予報回数が少なすぎるということになります．

今回の予測モデルでは，次のような結果となりました．
```
スレットスコア: 0.06
バイアススコア: 0.06
```

### ブライアスコア
**ブライアスコア**（Brier Score）は，降水確率の評価指標です．
定義は単純で，降水確率$f_t$と実況$o_t$（降水あり＝1，降水なし＝0）の差の二乗平均です．

```math
\frac{1}{N}\sum_{t=1}^{N}(f_t-o_t)^2
```

要は二乗平均誤差です．この指標は天気予報以外ではあまり聞いたことがありません．

なお，今回の予測モデルでのブライアスコアは
```
ブライアスコア: 0.16
```
となりました．

# ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています．

本記事のコードは`06-forecast_eval`の中に入っています．

https://github.com/matsuda-tkm/qiita-calender-2023