　近年、豪雨の被害に関するニュースをよく耳にするように思います。そこで、極値統計と離散の状態空間モデルを組み合わせて、1日の降水量の年間最大値が増加しているかを調べてみます。（例えば，https://www.tandfonline.com/doi/full/10.1080/01621459.2020.1750414 のように，確率過程の極値の分布も多くの研究がありますが，今回は実装が簡単な状態空間モデルを用います．）

　一般化極値分布のStanでの実装は、このサイトを参考にさせていただきました。降水量のデータは気象庁のサイトからダウンロードしました。上位$r$個の確率分布については、西郷達彦・有本彰雄『Rによる極値統計学』を参考にしました。

モデリング

　次のグラフは、1876年から2020年までの東京の最大日降水量の推移とヒストグラムです。

　左のグラフを見ると、わずかに上昇傾向にあるようにも見えます。今回はこの最大日降水量が増加傾向にあるか調べます。

　まず、独立同一分布に従う確率変数の最大値は、サンプルサイズが十分に大きいとき、極値分布と呼ばれる分布族に属する分布に従うことが知られています。極値分布の一種である一般化極値分布の累積分布関数は次です。

F(x ; \mu, \sigma, \xi)=\exp \left\{-\left[1+\xi\left(\frac{x-\mu}{\sigma}\right)\right]^{-1 / \xi}\right\}

　下にパラメータを変化させた一般化極値分布の確率密度関数を示しました。

おおよそ、$\mu$は分布の位置を変化させ、$\sigma$は分布の広がりを変化させ、$\xi$は分布の形状を変化させることが分かります。

　以上を踏まえて、年間の最大日降水量を次のようにモデル化します。

\mu_t \sim \mathcal{N}(\mu_{t-1}, 1) \quad t=2, \ldots, T \\
y_t \sim GEV(\mu_t, \sigma, \xi) \quad t=1, \ldots, T

$y_t$が$t$年の最大日降水量であり、$GEV(\mu, \sigma, \xi)$は一般化極値分布です。このモデルは、一般化極値分布のパラメータ$\mu$を状態変数として、プロセスノイズを正規分布、観測ノイズを一般化極値分布とした状態空間モデルと理解することができます。状態変数が急激に変化することはないと考え、プロセスノイズの標準偏差を$1$としています。

解析

　パラメータと状態変数をMCMCにより推定するコードは以下です。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pystan

code = """
    functions {
        real gev_lpdf (real y, real mu, real sigma, real xi){
            real z;
            z =  1 + (y - mu) * xi/sigma;
            return -log(sigma) - (1 + 1/xi)*log(z) - pow(z,-1/xi);
        }
    }
    data {
        int<lower=0> N;
        vector[N] y;
    }
    parameters {
        real<lower = 0> mu[N];
        real sigma;
        real<lower = 0, upper = 1>  xi;
    }
    model {
        for (i in 1:N)
            y[i] ~ gev(mu[i], sigma, xi);
        for (i in 2:N)
            mu[i] ~ normal(mu[i-1], 1);
    }
    generated quantities {
        real p[N];
        real predict[N];
        for (i in 1:N) {
            p[i] = uniform_rng(0, 1);
            predict[i] = (xi != 0) ? (mu[i] - (sigma/xi) *(1-(-log (p[i]))^ -xi)) : (mu[i] - sigma*-log(p[i]));
        }
    }
"""

# df['降水量']は最大日降水量のデータ
data = {
    'N': len(df['降水量']),
    'y': df['降水量'].values
}

sm = pystan.StanModel(model_code=code)
fit = sm.sampling(data=data, iter=1000, chains=4)

　すべてのパラメータでRhatが$1.0$となり、よく収束しました。

　結果を可視化してみます。

la = fit.extract(permuted=True)
quar1 = np.quantile(la['predict'], 0.05, axis=0)
quar2 = np.quantile(la['predict'], 0.5, axis=0)
quar3 = np.quantile(la['predict'], 0.95, axis=0)

fig, axes = plt.subplots(figsize=(12, 4), ncols=2)
fig.suptitle("Observed Data and Prediction")
axes[0].set_title("Trainsition")
axes[0].plot(df['年'], df['降水量'], label='Obeserved')
axes[0].plot(df['年'], quar2, c='tab:orange', label='Median')
axes[0].fill_between(df['年'], quar1, quar3, facecolor='tab:orange', alpha=0.3, label='90% Interval')
axes[0].set_ylabel('Rainfall[mm]')
axes[0].set_xlabel('Year')
axes[0].legend()
axes[1].set_title("Median of Prediction")
axes[1].plot(df['年'], quar2, c='tab:orange')
axes[1].set_ylabel('Rainfall[mm]')
axes[1].set_xlabel('Year')
plt.show()

　左のグラフは、観測値と予測値の中央値・90%区間を示したものです。それらしい結果が出ていることがわかります。また、左の図では、予測値の変化が分かりにくいので、右側に予測の中央値の推移のみを示しました。この図から、近年では最大日降水量が増加傾向にあり、1960年代と比べると10mm程度増加していることがわかります。

　同様の解析を、福岡と札幌のデータに適用した結果が以下です。

　以上の結果から、福岡と札幌でも最大日降水量は増加傾向にあると考えられます。

上位3個のデータを用いたパラメータ推定

　上では、各年の最大値のみを用いてパラメータを推定していましたが、最大値はばらつきやすく、推定量の分散が大きくなります。この問題を解決するために、上位$r$個のデータを用いるパラメータ推定手法があります。

　上位$r$個を上から順に$y_{1}, y_{2}, \ldots, y_{r}$とすると、上位$r$個の確率密度関数は次です。

f\left(y_{1}, y_{2}, \ldots, y_{r}\right) = \begin{cases} \frac{1}{\sigma^{r}}\exp \left[ \left\{1+\xi\left(\frac{y_{r}-\mu}{\sigma}\right)\right\}^{-\frac{1}{\xi}} \right] \prod_{i=1}^{r}\left\{1+\xi\left(\frac{y_{i}-\mu}{\sigma}\right)\right\}^{-\frac{1}{\xi}-1}, & \xi \neq 0 \\\ \frac{1}{\sigma^{r}} \exp \left[-\exp \left\{-\left(\frac{y_{r}-\mu}{\sigma}\right)\right\}\right] \prod_{i=1}^{r} \exp \left\{-\left(\frac{y_{i}-\mu}{\sigma}\right)\right\}, & \xi=0
\end{cases}

　そこで、上のモデルを拡張して、

y_{t1}, y_{t2}, \ldots, y_{tr} \sim rGEV(\mu_t, \sigma, \xi) \quad t=1, \ldots, T

とします。$y_{t1}, y_{t2}, \ldots, y_{tr}$は$t$年の日降水量の上位$r$個であり、$rGEV$は上位$r$個の確率分布です。
　また、$\mu_t$の推移の標準偏差も推定することにし、切断正規分布を事前分布として設定します。すなわち、

\mu_t \sim \mathcal{N}(\mu_{t-1}, \sigma_\mu) \quad t=2, \ldots, T \\
\sigma_\mu \sim \mathcal{N}_+(0, 2)

とします。ここで、$\mathcal{N}_+$は、$0$で切断された正規分布です。

　ちなみに、今回のデータには時系列の相関があるため、例えば、連日豪雨が続くと上位の値が偏る可能性があります。しかし、東京の観測史上最大の1958年のデータを見ても、最大を記録した前日や翌日の降水量はせいぜい$30mm$程度であることから、上位数個のデータには大きな影響がないと考えられます。

　これを実装すると次のようになります。今回は$r=3$としました。

code = """
    functions {
      real rgev_lpdf (vector y, int r, real mu, real sigma, real xi){
          real logp;
          logp =  -r*log(sigma) - pow(1+(y[r]-mu)*xi/sigma, -1/xi);
          for (i in 1:r) {
              real z;
              z =  1 + (y[i] - mu) * xi/sigma;
              logp += - (1 + 1/xi)*log(z);
          }
          return logp;
      }
    }
    data {
        int<lower=0> N;
        int<lower=0> r;
        vector[r] y[N];
    }
    parameters {
        real<lower = 0> mu[N];
        real<lower=0> sigma;
        real<lower = 0, upper = 1>  xi;
        real<lower=0> sigma_mu;
    }
    model {
        for (i in 1:N)
          y[i] ~ rgev(r, mu[i], sigma, xi);
        for (i in 2:N)
          mu[i] ~ normal(mu[i-1], sigma_mu);
        sigma_mu ~ normal(0, 2);
    }
    generated quantities {
        real p[N];
        real predict[N];
        for (i in 1:N) {
            p[i] = uniform_rng(0, 1);
            predict[i] = (xi != 0) ? (mu[i] - (sigma/xi) *(1-(-log (p[i]))^ -xi)) : (mu[i] - sigma*-log(p[i]));
      }
    }
"""
topr = 3
str_topr = [f'Top{i+1}' for i in range(topr)]

# df_topr[str_topr]は、各年の上位r個のデータ
data = {
    'N': len(df_topr),
    'r': topr,
    'y': df_topr[str_topr].values
}
sm = pystan.StanModel(model_code=code)
fit = sm.sampling(data=data, iter=1000, chains=4, seed = 1234)

　可視化のコードは上と同様のため省略します。

　上とほぼ同様に、近年では年最大日降水量が増加傾向にあると言えそうです。ただし、上位$r$個のデータを用いた解析は、最大値のみを用いた解析と比べ、Rhatがわずかに悪くなったため、より安定的に解析するためには、実装やモデリングをさらに工夫する必要がありそうです。

まとめ

　パラメータが緩やかに時間変化する一般化極値分布で、東京、福岡、札幌の最大日降水量をモデリングしたところ、最大日降水量は増加傾向にあることが分かりました。もっと多くの地点で調べなければ結論は出せませんが、日本全体で最大日降水量が増加傾向にあると言えそうです。

　上位$3$個のデータを用いた解析でも、類似の結果が得られましたが、より安定的に解析するためには実装やモデリングを工夫する必要がありそうです。