気象情報解析特論　第9回　統計検定と推定の考え方

担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import math
from matplotlib.colors import Normalize
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter
import matplotlib.ticker as mticker
import random

下準備：SST・降水のデータセットの読み込みと関数の定義

今回は，海面水温（SST）と降水（Precipitation）のデータを用います。

デトレンドした降水偏差のデータ： http://web.is.ocha.ac.jp/~tsubasa/env_info/detrended_preca_gpcp.npz

（ソース：https://psl.noaa.gov/data/gridded/data.gpcp.html ）

In [2]:
loadfile = 'detrended_ssta_OISST.npz' # デトレンド済みの海面水温偏差の入力ファイル名を定義
ssta_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します 
ssta = ssta_dataset['ssta'] # 海面水温偏差を変数sstaに保存
lon2_sst = ssta_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2_sst = ssta_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = ssta_dataset['y'] # 年(year)を変数yに保存
# 今回は1982年から2019年のデータを用いる
ssta = ssta[:, :, (1982 <= y)*(y <= 2019)]

loadfile = 'detrended_preca_gpcp.npz' # 海面更正気圧偏差の入力ファイル名を定義
preca_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します 
preca = preca_dataset['preca'] # 海面更正気圧偏差を変数precaに保存
lon2_prec = preca_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2_prec = preca_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = preca_dataset['y'] # 年(year)を変数yに保存
m = preca_dataset['m'] # 月(month)を変数mに保存

# 今回は1982年から2019年のデータを用いる
preca = preca[:, :, (1982 <= y)*(y <= 2019)]
m = m[(1982 <= y)*(y <= 2019)]
y = y[(1982 <= y)*(y <= 2019)] # y自身のサイズが変わってしまうので，yは一番最後に書き換えないとダメ

# 領域平均をとる関数
def aave(west, east, south, north, lon2, lat2, var = ssta):
    var = var[(west<=lon2[:, 1])*(lon2[:, 1]<=east),:, :]
    var = var[:, (south<=lat2[1, :])*(lat2[1, :]<=north), :]
    aave_var = np.nanmean(np.nanmean(var, 0), 0)
    return aave_var

# 回帰図を計算する関数    
def reg_map(ref_time, var):
    # z500のサイズをそれぞれ変数imt, jmt, tmtに保存
    [imt, jmt, tmt] = var.shape 

    # 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
    a_var = np.zeros((imt, jmt)) # 回帰係数a
    b_var = np.zeros((imt, jmt)) # 切片b
    
    # np.polyfitがエラーを吐かないようにするために，陸地の場所(nanが入っている)に一度ゼロを入れておきたい
    # （nanmeanのようなnanpolyfitという関数はない）
    # is_land_grids_3Dは，sstaの値がnanのところだけTrueが入っているような3次元配列（360x180x456）
    is_land_grids_3D = (np.isnan(var)==True) 
    var[is_land_grids_3D]=0

    # 回帰図の計算
    # 経度方向にimt回，緯度方向にjmt回forループを回す
    for ii in range(0, imt):
        for jj in range(0, jmt):
            [a_var[ii, jj], b_var[ii, jj]] = np.polyfit(ref_time, var[ii, jj, :], 1)

    # さっきゼロにしておいた陸地の場所にもう一度nanを戻す
    # is_land_grids_2Dは，sstaの1982年1月の値がnanのところだけTrueが入っているような2次元配列（360x180）
    var[is_land_grids_3D]=np.nan
    is_land_grids_2D = np.squeeze(is_land_grids_3D[:, :, 0])
    a_var[is_land_grids_2D]=np.nan
    b_var[is_land_grids_2D]=np.nan
    
    return a_var

# 日本周辺の偏差を描画する関数
def draw_field_map(draw_var, lon2, lat2, fig_title, \
                   lon_min=0, lon_max=360.1, lat_min=-90, lat_max=90, color_map = 'BrBG'):

    # 描画する枠を作る
    fig = plt.figure(figsize=(8.0, 4.0))
    
    c_lon = 180
    
    # 枠の中に絵を入れる（図を1枚しか書かないときは1, 1, 1でOK）
    ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(central_longitude=c_lon))
    
    # 降水っぽいカラーバーを指定
    cm = plt.get_cmap(color_map) 
    
    # カラーバーの範囲
    vmin = 0.05
    vmax = 0.5
    vint = 0.05

    # 色で塗られた等高線を描く
    cs = ax.contourf(lon2, lat2, draw_var, \
                    np.hstack([np.arange(-vmax, -vmin+0.001, vint), np.arange(vmin, vmax+0.001, vint)]), \
                    cmap=cm, norm=Normalize(vmin=-vmax, vmax=vmax),extend='both', \
                    transform=ccrs.PlateCarree())
                    #colors=['black'], transform=ccrs.PlateCarree())

    # 海岸線を書く（これが気象場を描画するときにめちゃくちゃ重要になる）
    ax.coastlines(lw=0.5,color='gray',resolution='50m')

    # 軸のラベルの間隔（写経でOK）
    dlon,dlat=60,30
    xticks=np.arange(60,360.1,dlon)
    yticks=np.arange(-60,60.1,dlat)
    ax.set_xticks(xticks,crs=ccrs.PlateCarree())
    ax.set_yticks(yticks,crs=ccrs.PlateCarree())

    # 軸のラベルのフォーマット（写経でOK）
    latfmt=LatitudeFormatter()
    lonfmt=LongitudeFormatter(zero_direction_label=True)
    ax.xaxis.set_major_formatter(lonfmt)
    ax.yaxis.set_major_formatter(latfmt)
    ax.axes.tick_params(labelsize=12)
    
    # カラーバーをつける
    plt.colorbar(cs, shrink=0.6)
    
    # タイトルをつける
    plt.title(fig_title)

    # 描画範囲の指定
    # 正距円筒図法であることを明示するために2つ目の引数が必要
    # 経度の指定では，c_lon±180の値以外を書くと全領域が表示される
    # 0, 360と書くと360が0と同じ数字だと思われて地図が棒になる
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
    
    plt.show()

回帰図で出たシグナルは全部信用していいの？

今回からは，何かを発見するための解析手法というよりは，何かを発見したあとにそれがどれだけ確かな結果なのかを確かめる方法を学びます。

いま，たとえば降雨偏差のNiño3.4に対する回帰図から，世界各地へのエルニーニョ南方振動(ENSO)の影響を調べようとしたとします。

In [3]:
nino34 = aave(190, 240, -5, 5, lon2_sst, lat2_sst)
a_prec_on_nino34 = reg_map(nino34, preca)
draw_field_map(a_prec_on_nino34, lon2_prec, lat2_prec, 'Regression map 1982-2019')

このとき，太平洋上に広がる湿潤偏差や，海洋大陸（インドネシア〜フィリピン〜パプアニューギニアくらいの場所をこう言います）上に広がる乾燥偏差は，とても強いシグナルなので，なんだか信じて良さそうな気がしてしまいます。

一方，この図で日本付近を拡大して，我が国はENSOでどのような影響を受けるかを調べてみると，どうなるでしょうか？

In [4]:
draw_field_map(a_prec_on_nino34, lon2_prec, lat2_prec, 'Regression map 1982-2019', \
               lon_min=125, lon_max=150, lat_min=30, lat_max=50)

この図だけから判断すると，「エルニーニョのときには関東南部は多雨傾向」と判断したくなります。

しかし，熱帯太平洋からこんなに遠く離れた地域の気候を判断するには，根拠が少し弱い気がします。本当に信じて良いのでしょうか？

今日は，そのような「微妙な偏差を信じていいか」を客観的に評価してみようというお話をします。

二分割テストと統計検定

もっともシンプルなテストは，今持っているデータを二つに分けて，二つとも同じ結果を得られればOK，とする方法です。

未来の30年も同じ結果になる？

神山が大学院に入ったのは，衛星観測が始まってから約30年後でしたので，30年分のデータしか持っていませんでした。そして当時，神山の指導教員であるDennis Hartmannからよく「君の解析結果は面白い。でも次の30年のデータでも同じ結果になるかな？」と聞かれていました。「そんなもん，30年後にならなければわからないじゃないか！」と一瞬は思います（笑）。

ですが，もし15年前に同じ解析をしていたならば，データを15年しか持っていないので，きっと15年分で解析したでしょう。そして，Dennisは「次の15年のデータでも同じ結果になるかな？」と聞いたかもしれません。ところがいま我々は，その次の15年のデータを持っています。当然，前半の15年と後半の15年という独立なデータで同じ解析をしてみて，もし全く同じ結果を得ることが出来たとすれば，「次の30年でもおそらく同じ結果を得るだろう」と言っても説得力が増します（二分割テスト）。 逆に，前半の15年でしか同じ結果を得られなければ，30年全体の結果だと思っていたものは前半の15年に強く影響されて現れただけの結果に過ぎないことが分かります。

以上の考察から分かるように，上記の質問でDennisが言いたかったのは，もちろん本当に次の30年を待てと言うことではなく，一つの時系列データで結果を得ただけで満足してはいけない ということなのです。

（上の例では，時系列を前半と後半に分けましたが，データを2つに分ける方法は無数にあります。偶数年と奇数年で分けるとか。また，三分割テストや四分割テストも同じように考えることもできます）

ENSOの降水影響でやってみる

前置きが長くなりましたが，二分割テストの重要性がわかっていただけたところで，実際にやってみましょう。

いまは1982-2019年の38年分のデータがありますから，1982-2000年の前半19年と，2001-2019年の後半19年に分けて，先ほどと同様の回帰図を描いてみましょう。

In [5]:
a_prec_on_nino34_1sthalf = reg_map(nino34[(1982 <= y)*(y <= 2000)], preca[:, :, (1982 <= y)*(y <= 2000)])
a_prec_on_nino34_2ndhalf = reg_map(nino34[(2001 <= y)*(y <= 2019)], preca[:, :, (2001 <= y)*(y <= 2019)])

draw_field_map(a_prec_on_nino34_1sthalf, lon2_prec, lat2_prec, 'Regression map 1982-2000')
draw_field_map(a_prec_on_nino34_2ndhalf, lon2_prec, lat2_prec, 'Regression map 2001-2019')

二つの独立したデータで，同じような降水の回帰図が出てきました。 これなら，この回帰図は概ね信用して良さそうに思えます。

ただ，よくみると二つの図で違う場所もありますね？日本付近を拡大してみましょう。

In [6]:
draw_field_map(a_prec_on_nino34_1sthalf, lon2_prec, lat2_prec, 'Regression map 1982-2000', \
               lon_min=125, lon_max=150, lat_min=30, lat_max=50)
draw_field_map(a_prec_on_nino34_2ndhalf, lon2_prec, lat2_prec, 'Regression map 2001-2019', \
               lon_min=125, lon_max=150, lat_min=30, lat_max=50)

冒頭の図では「エルニーニョのときには関東南部は多雨傾向」と判断したくなりましたが，こちらの図を見た後では印象が違うと思います。

ご覧の通り，「関東南部が多雨傾向だったと思ったのは，実は2000年以降のENSOイベントでたまたまそうだった」ということに過ぎず，前半19年ではそのような傾向がほとんど見られないことがわかります。

これを少し難しい言い方をすると「1982年-2019年の38年分のデータからそのまま回帰図を描いて得られた仮説が，二分割テストを行うことで棄却された」ということになります。このように，なんらかの統計的考察によって仮説の妥当性を評価する方法を統計的仮説検定あるいは単に統計検定(statistical test)と言います。

モンテカルロ法と統計的推定

上記の二分割テストの結果を見てみると，たとえば「エルニーニョのときは九州地方に雨が降りやすい」という仮説は棄却されません。

仮にこのシグナルが妥当だとして，実際にはどの程度の雨が降る傾向があると考えて良いのでしょうか？これをもう少し定量的に計算してみましょう。

まず真っ先に思いつく方法は，先ほどの回帰図について，九州の領域平均を計算することです。

In [7]:
# 九州地方を128E-132E, 31N-34Nと定義
a_prec_on_nino34_1sthalf_kyushu = aave(128, 132, 31, 34, lon2_prec, lat2_prec,\
                                       var = a_prec_on_nino34_1sthalf[:, :, np.newaxis])[0]
a_prec_on_nino34_2ndhalf_kyushu = aave(128, 132, 31, 34, lon2_prec, lat2_prec,\
                                       var = a_prec_on_nino34_2ndhalf[:, :, np.newaxis])[0]

# 領域平均。単位：(mm/day)/°C，つまりNiño3.4海域が1°C暖かくなると何mm/dayだけ雨が降りやすくなるか
print(a_prec_on_nino34_1sthalf_kyushu)
print(a_prec_on_nino34_2ndhalf_kyushu)

0.25135053821213016
0.24676175683968127

このように，Niño3.4海域が1°C暖かくなると，

・1982-2000年のデータだと0.251 mm/day

・2001-2019年のデータだと0.247 mm/day

だけ雨が振りやすくなるという結果になっています。

このように，サンプルとして取り出す月によって結果が変わってしまうことをサンプリング変動(sampling variability) といいます。そして，サンプリング変動の範囲内のどこかに「真の値」があると考えて，統計的考察でそれを見積もる方法が統計的推定です。

もっともシンプルな統計的推定は，モンテカルロ法 を用いる方法です。モンテカルロ法とは，サンプリングや数値計算などを乱数で行う手法の総称です。

ここでは，1982年-2019年の38年間の中から，無作為に19年分のデータを取り出して，先ほどと同様に回帰図の領域平均を計算するのを繰り返します（ここでは半分の大きさの部分データを取り出すが，必ずしも半分ではなくてもよい）。38年間から19年分のデータを取り出す場合の数は，

$$\ _{38} C_{19} = \frac{38!}{19!19!} = 35,345,263,800$$

で，約350億通りもあります。これらのうち無作為に1000通りを選んで，Niño3.4に対する回帰係数を計算して領域平均を計算することで，九州の降雨へのENSOの影響を見積もろうというのです（手計算だと絶対やりたくないですが，コンピュータは得意です！）。

では，早速やってみましょう。

In [8]:
# 回帰図の計算量を抑えるために，先にデータの領域を九州だけに絞っておく
preca_kyushu = preca[(128<=lon2_prec[:, 1])*(lon2_prec[:, 1]<=132),:, :]
preca_kyushu = preca_kyushu[:, (31<=lat2_prec[1, :])*(lat2_prec[1, :]<=34), :]
lon2_prec_kyushu = lon2_prec[(128<=lon2_prec[:, 1])*(lon2_prec[:, 1]<=132),:]
lon2_prec_kyushu = lon2_prec_kyushu[:, (31<=lat2_prec[1, :])*(lat2_prec[1, :]<=34)]
lat2_prec_kyushu = lat2_prec[(128<=lon2_prec[:, 1])*(lon2_prec[:, 1]<=132),:]
lat2_prec_kyushu = lat2_prec_kyushu[:, (31<=lat2_prec[1, :])*(lat2_prec[1, :]<=34)]

# 1000個の要素を持つ配列を用意しておく
a_prec_on_nino34_MC_kyushu = np.zeros(1000)

# 以下の操作を1000回繰り返す
for cc in range(0, 1000):
    # 1982から2019の自然数から重複のないように無作為に19個取り出す
    years_sample = random.sample(range(1982, 2020), 19)
    # years_sampleに含まれる年のみTrueとなる配列
    year_boolean = np.isin(y, years_sample)
    # years_booleanでTrueの月のデータのみを使って，回帰係数を計算
    a_prec_on_nino34_MC = reg_map(nino34[year_boolean], preca_kyushu[:, :, year_boolean])
    # 九州における回帰係数の領域平均を保存
    a_prec_on_nino34_MC_kyushu[cc] = aave(128, 132, 31, 34, lon2_prec_kyushu, lat2_prec_kyushu,\
                                          var = a_prec_on_nino34_MC[:, :, np.newaxis])[0]

それでは，散布図を描いてみましょう。

In [9]:
plt.scatter(range(0, 1000), a_prec_on_nino34_MC_kyushu)
plt.xlabel('# of test') # 横軸は何番目の試行か
plt.ylabel('Regression Coefficient [mm/day/°C]')# 縦軸は回帰係数
plt.show()

これを見ると，どの19年を取り出すかによって九州地方の降雨に対するENSOの影響がだいぶ違っていることがわかります。

この分布を度数分布表（ヒストグラム）に直してみましょう。

In [10]:
# 0.2から0.6まで，ビン幅0.05でヒストグラムを描く
plt.hist(a_prec_on_nino34_MC_kyushu, np.arange(-0.2, 0.6, 0.05))
plt.xlabel('Regression Coefficient [mm/day/°C]') # 横軸は回帰係数
plt.ylabel('Counts') # 縦軸は何回の試行が該当したか

# 下位2.5%（今の例だと下から25番目）の値を計算し，赤線で表示
percentile_25 = np.percentile(a_prec_on_nino34_MC_kyushu, 2.5)
print(percentile_25)
plt.vlines(percentile_25, 0, 200, 'r')
# 上位2.5%（今の例だと上から25番目）の値を計算，青線で表示
percentile_975 = np.percentile(a_prec_on_nino34_MC_kyushu, 97.5)
print(percentile_975)
plt.vlines(percentile_975, 0, 200, 'b')

plt.show()

0.015614772944002752
0.4273576631585275

このヒストグラムによると，下位2.5%と上位2.5%を外れ値とみなせば，九州地方における降雨のNiño3.4に対する回帰係数$a$は95%の確率で

$$0.016 \ \mathrm{mm/day/°C} < a < 0.427 \ \mathrm{mm/day/°C} $$

の範囲に入ることが統計的に推定できます！この範囲のことを， 回帰係数の95%信頼区間 といいます。

また今のケースでは，推定範囲が正の値に収まっているかどうかを確かめることにもなっているので，統計検定を行ったと考えることもできます。今回の状況の場合は，モンテカルロテストの結果，九州地方の降雨には95%の信頼度で有意なENSOの影響がある などと言います。

（逆に，もし推定範囲が正の値に収まっていない場合は，「九州地方の降雨には有意なENSOの影響は認められない」という結論になります。）

ただし，統計的に有意な結果が出ればそれで良いというわけではありません。そもそも，95%という数字は（研究業界では非常によく使われますが）それほど特別な意味があるわけではありません。実際，最近では95%の信頼度で有意でありさえすれば論文が出てしまう風潮を憂えて，「信頼区間を示しましょう」という流れもあるようです（詳しく知りたい方は「ASA声明」で検索してください）。

実際の研究では，統計検定の結果に頼るだけでなく，起こっている物理のメカニズムを理論やデータ解析で検証したり，数値シミュレーションで再現したりするなど，様々な方法で仮説を検証していくことが大切です。

課題9（締切：次回の授業開始前まで）

この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。

A問題. アメリカ南東部のメキシコ湾沿いの地域（いわゆるGulf Statesと呼ばれる場所）の降雨は，ENSOの影響を受けていると言えますか？二分割テストの結果から考察してください。

B問題. 中国華南地方（福建省，海南省，広東省を含む地域）の降雨は，どの程度ENSOの影響を受けているでしょうか？モンテカルロ法によって，回帰係数の領域平均を95%信頼区間で推定することによって考察してください。

C問題. 下図のSSTのNiño3.4への回帰図を見ると，「北大西洋の海面水温はエルニーニョ時に冷える傾向にある」という仮説が立ちます。この仮説の妥当性を確かめてください。

D問題. 神山がシアトルに住んでいたとき，シアトルで雨が降っている日は東京でも雨が降っていることが多いような気がしました。気のせいですか？

提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「気象情報解析特論9，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。

In [11]:
# C問題の回帰図
a_sst_on_nino34 = reg_map(nino34, ssta)
draw_field_map(a_sst_on_nino34, lon2_sst, lat2_sst, 'Regression map 1982-2019', color_map = 'seismic')