1.5.11.1. Sprogø 気象局の最大風速予測

練習問題の目的は50年間の最大風速を予測することです, ただし, 予測はそのような期間に渡る測定がされていない場合に行ないます. 入手可能なデータはデンマークの Sprogø 気象局での21年間のデータだけです. まず最初に統計的な手段を与えて, scipy.interpolate モジュールから関数を描きます. そして最後に興味ある読者は全く異なる方法を使い結果を生データから計算することになるでしょう.

1.5.11.1.1. 統計的方法

年間の最大値は正規分布にフィットできると考えられます. しかしこの関数は予測に使うことはできません, なぜならこの関数は風速の最大値から確率を与える関数だからです. 50年間の最大値を知るには逆の方法が必要となり, 定義された確率から結果を求める必要があります. これが分位点関数の役割です, そして目的はこの関数を見つけることになります. 現在のモデルでは, 50年間の年間最大値は上位 2 % の分位点として定義されると考えられます.

定義から分位点関数は累積分布関数の逆関数です. 後者は年間の最大値の分布関数を表わします. 練習問題では, i 年の累積確率 p_i を p_i = i/(N+1), 測定年数 N = 21 とします. これで各年で測定された最大風速の累積確率を計算できます. scipy.interpolate モジュールを使ってこれらの測定点から分位点関数をフィットするのが有効です. 最終的に50年の最大値は 2% 分位点の累積分布関数から評価できます.

1.5.11.1.2. 累積確率の計算

年間風速の最大値は, 既に計算されて max-speeds.npy ファイルに numpy 形式で保存されています, なので numpy を使って読み込みます:

>>> import numpy as np

>>> max_speeds = np.load('intro/summary-exercises/examples/max-speeds.npy')

>>> years_nb = max_speeds.shape[0]

前の節の累積確率 p_i の定義に従って, 対応する値は:

>>> cprob = (np.arange(years_nb, dtype=np.float32) + 1)/(years_nb + 1)

そして, これらは与えられた風速をフィットすると仮定されます:

>>> sorted_max_speeds = np.sort(max_speeds)

1.5.11.1.3. UnivariateSpline による予測

この節では分位点関数を, 点からスプラインを表現する UnivariateSpline クラスを使って予測します. デフォルトの挙動では3次スプラインを使い, 点は信頼度に応じて異なる重みが与えられます. 変種に InterpolatedUnivariateSpline と LSQUnivariateSpline があり, これらは誤差の抑制に違いがあります. 2次元のスプラインが欲しい場合には BivariateSpline クラスファミリーが提供されます. これらの1次元, 2次元スプラインは FITPACK Fortran サブルーチンを使います, その結果としてスプラインを表現, 評価する splrep, splev を経由する低レベルライブラリーにアクセスすることができます. またより単純な目的のために FITPACK パラメータを使わない補間関数も提供されています(interp1d, interp2d, barycentric_interpolate を見てください).

3次スプラインでデータを正確にフィットできているようなので, Sprogø の最大風速のために UnivariateSpline を使います:

>>> from scipy.interpolate import UnivariateSpline

>>> quantile_func = UnivariateSpline(cprob, sorted_max_speeds)

さてこれで分位点関数が確率の全領域に渡って評価されます:

>>> nprob = np.linspace(0, 1, 1e2)

>>> fitted_max_speeds = quantile_func(nprob)

現在のモデルでは50年間の最大風速は上位 2% の分位点で定義されます. 結果として累積確率の値は:

>>> fifty_prob = 1. - 0.02

その結果50年の間に起きる暴風の風速が推測されます:

>>> fifty_wind = quantile_func(fifty_prob)

>>> fifty_wind      

array(32.97989825...)

結果を Matplotlib の figure にまとめます:

intro/summary-exercises/auto_examples/images/plot_cumulative_wind_speed_prediction_1.png

解答: Python source file

1.5.11.1.4. 演習 Gumbell 分布

興味ある読者は21年間の風速測定データを使って練習問題を作りたいと思っているでしょう. 測定周期は90分です(元の周期は10分でしたが練習問題の準備を簡単にするためファイルサイズを小さくしました). データは sprog-windspeeds.npy に numpy 形式で保存されています. 練習問題が終わるまで, 作図のソースコードはみないで下さい.

  • まずは年間での最大値を入手して matplotlib の棒グラフで表示しましょう.

intro/summary-exercises/auto_examples/images/plot_sprog_annual_maxima_1.png

解答 Python source file

  • 次は累積確率 p_i に対して定義される Gumbell 分布 -\log( -\log(p_i) ) を線形分位点関数のフィッティングに使ってみましょう(UnvariateSpline の次数を定義するのを忘れないように). 年間の最大値と Gumbell 分布は以下のような図になるはずです.

intro/summary-exercises/auto_examples/images/plot_gumbell_wind_speed_prediction_1.png

解答 Python source file

  • 最後に50年の間に起きる最大風速として 34.23 m/s を得るはずです.