気象情報解析特論　第10回　平均値の差のt検定

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

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
from scipy import signal
from scipy.stats import t

下準備：東京の気温のデータセットの読み込み

今回は，東京の気温のデータを用います。

In [2]:
# 東京の気温のデータセットの読み込み
tokyo_temp = np.genfromtxt("Tokyo_temp.csv",  # ファイルのパスを書く
                  delimiter=",",    # 区切り文字
                  usecols=(0, 1, 2) # 読み込みたい列番号
                 )
y = tokyo_temp[:, 0]
m = tokyo_temp[:, 1]
temp = tokyo_temp[:, 2]

# 今回は1880年から2019年の140年分のデータを用いる
temp = temp[(1880 <= y)*(y <= 2019)]
m = m[(1880 <= y)*(y <= 2019)]
y = y[(1880 <= y)*(y <= 2019)]

# 気候値の計算（基準は140年の各月の平均）
temp_clim = np.zeros((12))
for mm in range(1, 13): 
    temp_clim[mm-1] = np.nanmean(temp[m==mm], 0)

# 偏差の計算
tempa = np.zeros((temp.shape))
for yy in range(1880, 2020):
    for mm in range(1, 13):
        tempa[(y==yy)*(m==mm)] = temp[(y==yy)*(m==mm)] - temp_clim[mm-1]
        
# 今回は温暖化に興味があるのでデトレンドはしない

昭和時代と平成時代で比べると，東京は温暖化した？

たとえば皆さんが，昭和生まれのおじさんやおばさんに「昔はもっと涼しくてねぇ」という話をされたとします（ちなみに神山もギリギリ昭和生まれです）。このとき，皆さんはどう答えますか？「気のせいですよ」とか言いましょうか？

%%% あなた方もどうせ令和の子たちと同じ話をするのです，覚悟しましょう。

昭和時代（昭和35年-昭和64年つまり1960-1989年の30年とします）と平成時代（平成2年-平成31年つまり1990-2019年とします）では，東京の気候条件は変化したのでしょうか？

まず真っ先に思いつくのは，二つの時代で気温または気温偏差の平均をとって比べることです（簡単ですが一応コンポジット解析ですね）。

In [3]:
## 1960-1989年と1990-2019年の気温偏差データを抜き出す
tempa_1960 = tempa[(1960 <= y)*(y <= 1989)]
m_1960 = m[(1960 <= y)*(y <= 1989)]
y_1960 = y[(1960 <= y)*(y <= 1989)]
tempa_1990 = tempa[(1990 <= y)*(y <= 2019)]
m_1990 = m[(1990 <= y)*(y <= 2019)]
y_1990 = y[(1990 <= y)*(y <= 2019)]

print("Showa", np.mean(tempa_1960, 0))
print("Heisei", np.mean(tempa_1990, 0))

Showa 0.5838492063492068
Heisei 1.5663492063492066

つまり，1880年から2019年の平均を基準として，昭和が0.58°C，平成が1.57°C暖かかったということです。確かに昭和時代は涼しかったということで，おじさんおばさんの感覚は正しかったようです。

しかしこの数字から，なんらかの外的要因があって気候が変化したのか，気候条件は完全に同一だったにも関わらずたまたまランダムな変動で平均値に差が生まれたのか については何も判断できません。たとえば，「サイコロを30回振る」という試行を全く同じサイコロで2回繰り返したとしても，それぞれの30回で出た目の平均値には若干の差があるでしょう。しかし，平均値にそのような差が存在するからといって，2回の試行の間にサイコロに細工がされたとは普通は考えないですよね。

このように，全く外的要因が同一だったとしても，統計的乱雑さ（stochasticity）によって差が生まれることは日常的にも多いです。同じ学力の高校生が入試を受けても，その日の体調や出題の得意不得意など，「運要素」「ガチャ」と言われるような，コントロール不可能なstochasticityで合否が分かれることもあるでしょう。

話は戻って，気温偏差のグラフを見てみましょう。

In [4]:
## 1960-1989年の気温を青，1990-2019年の気温を赤で描画
plt.plot(tempa_1960, 'b', label = '1960-1989')
plt.plot(tempa_1990, 'r', label = '1990-2019')
plt.plot(0*tempa_1990, 'k--', label = 'Mean of 1880-2019')
plt.ylabel('Temperature Anomalies in Tokyo (°C)')
plt.legend()
plt.show()

実際，このグラフを見ると，なんとなく昭和時代と比べて平成時代の温度は上昇しているように見えます。しかし前回の授業で，たとえ気候条件がほぼ同じと考えられる30年間でも，そのうちどの15年をサンプルとして取り出すかによって回帰係数がかなり変動しうる（=サンプリング変動）ことを確かめました。そう考えると今回のケースでも，サンプリング変動のせいで取り出した30年にたまたま差があるように見えるだけなのか，本当に気候条件になんらかの外的な要因（=外部強制）に変化があったのかについては，これを見ただけではわかりません。

そこで，平均値の差について統計検定を行うことで，もし昭和時代と平成時代の気候条件が同じだったと仮定すると，サンプリング変動でこの程度の差が生まれてしまうことはありえない（と考えていいくらい確率の低いことである）のか をテストします。

平均値の差のt検定

前回のようにモンテカルロテストをやっても良いのですが，モンテカルロテストをするにはある程度大量のデータがある必要があります。そこで，今回はより汎用的な手法である t検定 を勉強します。

t検定は，母集団（=気象の文脈では，外部強制が同一の気候の時代が無限に続くと考えたときの気象データ全体） の分布が正規分布に従うと考えて良い場合のうち，母集団から取り出せるサンプル（標本）の数（=実際に利用できる独立なデータの数）が30程度よりも少ない場合に用いることができる強力な方法です。母集団が正規分布に従うなら，その正規分布の両側2.5%にデータがはみ出すような現象は起こりにくい，はみ出さない現象は偶然でも起こりうる，のように「偶然だと仮定した際の起こりにくさ」を判断することができます。ただし，サンプルが少ない場合，そのサンプルから見積もることのできる統計量は，正規分布ではなくt分布に従うので，そちらと比較することになります（詳しくは学部の統計学を復習してください）。

t検定の手順

t検定の具体的なやり方は次の通りです。

I. 要求する有意水準を決める

有意水準とは，どの程度の確率で起こることを「十分確率が低い」と判定するかの基準です。大気海洋研究コミュニティでは，多くの場合0.05（信頼度95%）を用います。しかし，降水のようにノイズが大きいデータや，得られるサンプル数自体が限定的な場合は，0.1（信頼度90%）で妥協することもあります。

発展的な注：有意水準について詳しく知りたい場合は，「第一種の過誤」についてググってください。

II. 帰無仮説をつくる

帰無仮説とは，棄却するためにわざと作る仮説です。今回の場合だと，「昭和時代と平成時代の気温の母平均の差はゼロである」です（母平均=母集団の平均のこと）。

III. 対立仮説をつくる

帰無仮説を棄却することによって，示したい仮説のことです。帰無仮説の否定となるようにつくります。今回の場合だと，「昭和時代と平成時代の気温の母平均の差はゼロでない」です。

IV. 検定統計量を計算する

帰無仮説が正しいと仮定するとt分布に従う量（=検定統計量）を計算します。この検定統計量は，検定したい状況に応じて計算すべき数式が違うので，統計学の教科書などを適宜参照してください。

今回の状況では，

・昭和時代の30年間の各月と平成時代の30年の各月には対応がない

・昭和時代の気温の分散と平成時代の気温の分散は同じとは限らない

というケースに該当すると考えられるので，

$$t = \frac{\overline{x_1}-\overline{x_2}}{\sqrt{\frac{\sigma_1^2}{N_1^*}+\frac{\sigma_2^2}{N_2^*}}}$$

という量を計算します。ただし，$\overline{x_1}$と$\overline{x_2}$はサンプル平均，$\sigma_1$と$\sigma_2$はサンプル標準偏差，$N_1^*$と$N_2^*$は実効的サンプル数（後述）で，添字の1と2はそれぞれ平成時代と昭和時代の気温データから計算された量を表すことにします。

V. 検定統計量と棄却域を比較し，帰無仮説を棄却するかどうかを判断する

棄却域は「もし帰無仮説が正しいなら，これ以上検定統計量tがt分布の中心から外れる確率は低いので，帰無仮説は正しくないのだろう（=帰無仮説を「棄却」しよう）」と判断する領域のことです。棄却域は，あらかじめ定めておいた有意水準と統計的自由度（用いたサンプルの合計-2）についてのt分布関数から計算しますが，pythonの組み込み関数で求めることができます。

その比較によって「帰無仮説を棄却する」と判断した場合は，対立仮説を採択します。逆に，「帰無仮説を棄却しない」と判断したとした場合は，対立仮説を採択することはできません（ちなみにこの場合でも，帰無仮説を「棄却しない」だけであって，ただちに帰無仮説の方を「採択できる」わけではありません）。

今の場合は，検定統計量tが棄却域を超えた場合のみ，対立仮説である「昭和時代と平成時代の気温の母平均の差はゼロでない」を採択します。すなわち，サンプリング変動のせいでこのような気温の差が生まれる確率は十分低いと考えられるので，東京の気候は変化しているのだろう，と結論づける ことになります。

時系列データで検定する際に，サンプル数に関する注意

今回の文脈では，母集団から360ヶ月分のサンプルを取り出しているので，独立なサンプルを360個持っていると言いたくなります。

しかし，時系列データがレッドノイズ的である場合は，前の年の情報を引きずるので，個々の年は独立なサンプルであるとは言えません。

たとえば，ローパスフィルタリングをする前とした後の時系列では，持っている情報の量が全然違うのがわかるでしょう。

In [5]:
def butterworth_lowpass(time_series, delta_t, f_cut, order):
    b, a = signal.butter(order, f_cut, btype='low', fs = 1/delta_t)
    time_series_lp = signal.filtfilt(b, a, time_series)
    return time_series_lp

# 例として年々の8月の気温を見てみる
tempa_1990_Aug = tempa_1990[m_1990==8]
# 年々の8月の気温のうち5年周期より長いものをローパスしたもの
tempa_1990_Aug_low = butterworth_lowpass(tempa_1990_Aug, 1, 0.2, 1)

plt.plot(tempa_1990_Aug) # 年ごとにある程度の寒暖の独立な情報がある
plt.plot(tempa_1990_Aug_low) # 数年に1回程度くらいしか寒暖の独立な情報がない
plt.show()

この図を見るとわかるように，8月の気温は年ごとにある程度の寒暖の独立な情報がありますが，それをローパスすると数年に1度くらいしか寒暖の独立な情報がありません。すなわち，強いレッドノイズであればあるほど，時系列が持っている独立な情報が少ないのです。

そこで，実効的サンプル数 を求める以下の公式を用います（Bretherton et al. 1999）。

$$N^* = N \frac{1-r}{1+r}$$

ここで，$N^*$はレッドノイズであることを考慮した後の「実効的サンプル数」，$N$はレッドノイズであることを考慮する前の「元のサンプル数」（いまの例だと360個），$r$は時系列のタイムステップを1つずらした時系列と元の時系列の相関係数（「ラグ1自己相関」）です。

（ちなみに上記の公式は，平均値の検定を行うときのみに用いることができるもので，分散や相関係数を検定する際に用いるべき公式は異なります。）

たとえば仮に，年々の気温が理想的なホワイトノイズ（=前年の気温とは全く独立に翌年の気温が決まる）だった場合を考えてみると，気温の時系列を1年ずらしたものとの相関，すなわちラグ1自己相関$r$は0になりますね。この場合，実効的サンプル数は元のサンプル数に一致します。実際は，気温はレッドノイズ的ですので，$r$は正の値となり，実効的サンプル数は元のサンプル数よりも小さい値となります。

実際にやってみる

それでは以上を踏まえて，1960-1989年と1990-2019年の平均気温に統計的に有意な差があるかどうかをt検定で確かめてみましょう。

I. 有意水準：0.05（信頼度95%）

II. 帰無仮説「昭和時代と平成時代の1月気温の母平均の差はゼロである」

III. 対立仮説「昭和時代と平成時代の1月気温の母平均の差はゼロでない」

IV. 検定統計量tは，以下のように計算：

In [6]:
## 検定統計量を計算するのに必要な量を計算する関数。引数は時系列。
def calc_statistic(x):
    # サンプル平均（=母集団から360ヶ月分のサンプルを取り出して平均をとったもの）
    x_bar = np.mean(x, 0)
    # サンプル標準偏差（=母集団から360ヶ月分のサンプルを取り出して標準偏差を計算）
    sigma = np.std(x) 
    # 実効的サンプルサイズ（Bretherton et al.の公式によって計算）
    N = x.shape[0]
    r = np.corrcoef(x[0:-1], x[1:])[0, 1]
    N_star = N * (1-r)/(1+r)
    # 中を見たければコメントアウトを外す
    # print('x_bar, sigma, r, N, N_star are')
    # print([x_bar, sigma, r, N, N_star])
    return x_bar, sigma, N_star

## 検定統計量を求める関数
def calc_t(x1, x2):
    # 平成時代の気候の統計量
    x_bar1, sigma1, N_star1 = calc_statistic(x1)
    # 昭和時代の気候の統計量
    x_bar2, sigma2, N_star2 = calc_statistic(x2)
    # 検定統計量
    t_sample = (x_bar1 - x_bar2) / (np.sqrt(sigma1**2/N_star1 + sigma2**2/N_star2))
    return t_sample

t_sample = calc_t(tempa_1990, tempa_1960)
print('\n t =', t_sample)

 t = 8.523360770349363

V. 棄却域は，以下のように計算：

95%信頼度で検定する場合は，両側に2.5%ずつ外れ値を仮定するので（両側検定），t.ppfの引数は0.975となることに注意します。

In [7]:
## 比較するt分布の統計的自由度の計算
x_bar1, sigma1, N_star1 = calc_statistic(tempa_1990)
x_bar2, sigma2, N_star2 = calc_statistic(tempa_1960)
# サンプルサイズの合計から2を引いたものが統計的自由度
dof = N_star1 + N_star2 - 2

# 有意水準と統計的自由度から，棄却域の計算
t_975 = t.ppf(0.975, dof)

print('\n t_975 =', t_975)

 t_975 = 1.9667480845443612

これにより，棄却域は

$$|t| > t_{975} = 1.966...$$

となります。この棄却域に，求めた検定統計量$t$が入っている場合は，帰無仮説を棄却します（仮に検定統計量$t$がt分布に従っていると仮定すると，分布の中心から外れすぎており，そのようなことが起こる確率は十分低いと考えられるため）。

今回の場合は，検定統計量$t = 8.523$が棄却域に入っているため，帰無仮説を棄却し，対立仮説「昭和時代と平成時代の1月気温の母平均の差はゼロでない」を採択します。

すなわち，少なくとも統計的には，サンプリング変動のせいでこのような気温の差が生まれる確率は十分低いと考えられるので，東京の気候は変化していると結論づけることができるのです。

コンポジット図の統計検定

最後に上記の考え方を持いて，コンポジット図の有意性についても同様な統計検定を行って考えてみましょう。

SSTデータセットの読み込み，関数の定義，コンポジット図の計算

今回は，環境情報論第7回で扱ったようなエルニーニョのコンポジット図を考えます。

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

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

def aave(west, east, south, north, 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 draw_field(field, fig_title, vmin = -6, vmax = 6, vint = 1):
        
    plt.figure(figsize=(8.0, 4.0))
    cm = plt.get_cmap('seismic')
    cs = plt.contourf(lon2, lat2, field,\
                  cmap=cm, norm=Normalize(vmin=vmin, vmax=vmax),\
                  levels=np.arange(vmin,vmax+vint,vint), extend='both')    
    plt.colorbar(cs)
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    title = fig_title
    plt.title(title)
    plt.xlim(150, 300)
    plt.ylim(-60, 60)

nino34 = aave(190, 240, -5, 5)
draw_field(np.mean(ssta[:, :, nino34>1], 2), \
           'Composited SSTA for months when Niño3.4 > 1 ℃', -3, 3, 0.5)
plt.show()

これは，Niño3.4が1°Cを超えているときのみの海面水温偏差を平均して求めた，エルニーニョ現象のコンポジット図です。

各点についてt検定を繰り返し実行する

上記のコンポジット図のうち，サンプリング変動ではなく，本当にエルニーニョの影響を受けている地点にのみ印をつけてみましょう。すなわち以下のt検定を，全ての地点について繰り返し実行します。

I. 95%信頼度

II. 帰無仮説「コンポジットをとった偏差は0と差がない」

III. 対立仮説「コンポジットをとった偏差は0と差がある」

IV. 検定統計量は，先ほどの例で$x_2$の方が常にゼロであるような時系列（ただしゼロで固定なので，サンプル数$N_2$もゼロ）を考えれば，

$$t = \frac{\overline{x_1}}{\sqrt{\frac{\sigma_1^2}{N_1^*}}}$$

となります。

V. 棄却域は，自由度が$N_1^*-2$となるように各点でt_975を計算し，検定統計量tと比較します。

In [9]:
## 時系列x1に対して，帰無仮説の棄却(True)か棄却しない(False)を返す関数
def t_test(x1):
    # 海面水温の統計量
    x_bar1, sigma1, N_star1 = calc_statistic(x1)
    # 検定統計量
    t_sample = x_bar1 / (np.sqrt(sigma1**2/N_star1))
    # サンプルサイズの合計から2を引いたものが統計的自由度
    dof = N_star1 - 2
    # 有意水準と統計的自由度から，棄却域の計算
    t_975 = t.ppf(0.975, dof)
    # 棄却域に入る時だけTrueを返す
    test_result = (np.abs(t_sample) > t_975)
    return test_result

## 準備
imt, jmt, tmt = ssta.shape
test_result = np.array([-9999, -9999]) # テキトーな2次元ベクトルを用意

## t検定
# 経度方向にimt（=360）回，緯度方向にjmt（=180）回forループを回して，t検定を実行
for ii in range(0, imt):
    for jj in range(0, jmt):
        if t_test(np.squeeze(ssta[ii, jj, nino34>1])):
            # もしt_testで棄却域に入るならば，
            # 経度と緯度の入った配列をtest_resultにくっつける
            test_result = np.vstack((test_result, np.array([lon2[ii, 0], lat2[0, jj]])))
            
    if (ii % 60 == 0):
        print(ii) # ちゃんと計算が進んでるかチェックするために，30回に1回iiを出力する

# 最初につけたテキトーな2次元ベクトルを外す
test_result = test_result[1:, :]

0
60
120
180
240
300

test_resultの中には，t検定で「95%信頼度でエルニーニョの有意な影響を受けている」と判定された経度と緯度の組がベクトルとして入っています。

これらを散布図として描画することで，地図の上に重ねてみましょう。

In [10]:
# コンポジット図を描く
draw_field(np.mean(ssta[:, :, nino34>1], 2), \
           'Composited SSTA for months when Niño3.4 > 1 ℃', -3, 3, 0.5)
# t検定をパスした点の散布図を重ねる
plt.scatter(test_result[:, 0], test_result[:, 1], marker='.', s=0.5, c='black')
plt.show()

いかがでしたか？描画してみると，意外と通過している点は多くないことがわかると思います。少なくとも統計的には，これらの点のみを「エルニーニョの影響がある」点と思っておくのが良さそうです。

一方で，統計にとらわれ過ぎても物理の本質を見失うことがありますので，あくまでも「統計はエビデンスの一つ」くらいに思っておく のがいいでしょう。

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

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

A問題. 2000-2009年と，2010-2019年の各10年間では，東京の気温偏差に有意差はみとめられますか？

B問題. ラニーニャ現象(nino34<-1とでもすれば良いでしょう）について，上記と同様なt検定を行ってください。

C問題. あなたの趣味などに関する好きなデータを用いて，平均値の差のt検定を行ってください。

D問題. 「分散の差のF検定」について調べて，上記と同様に東京の気温偏差の分散に適用してみてください。ただし，レッドノイズである（=自己相関を持つ）時系列の実効的サンプル数は，分散の検定を行う場合は以下の公式を用います（Bretherton et al. 1999）。

$$N^* = N \frac{1-r^2}{1+r^2}$$

ただし，$N^*$は「実効的サンプル数」，$N$は「元のサンプル数」，$r$は「ラグ1自己相関」です。

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