C. 信頼区間の推定について（コンポジット・トレンド）

　統計解析を行った結果にについては、統計的有意性の評価（エラーバーをつける）必要がある。

C.1 トレンドの信頼区間推定（2.2節の補足）

C.2 コンポジット解析の信頼区間推定（2.3節の補足）

【準備】 これまで用いてきたものと同じデータを用いる

In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import datetime
import scipy.stats as stats
from scipy.optimize import curve_fit
In [2]:
yst = 1980; yen = 2018
list = []
for yy in range(yst,yen+1,1):
    cyr = str(yy) # 文字列化
    file = "../../data/Kyoto/Kyoto_" + cyr  +".csv"
    df = pd.read_csv(file,encoding = "shift-jis",header=2,index_col=0,parse_dates=[0])
    df = df.iloc[2:,:] # ヘッダーの余計な部分を除く
    
    list.append(df) # listに要素(df)を追加していく

# pandas.concatで結合する
dfm = pd.concat(list,axis=0) # axis=0 (時間軸方向）に結合
    
temp = dfm['平均気温(℃)']

C.1 トレンドの信頼区間推定（2.2節の補足）

　t-検定を用いて2.2節で求めたトレンド（回帰係数）の信頼区間（エラーバー）を推定してみよう。 求まった回帰係数の標準誤差: s.e.($\beta$)は

\begin{align} s.e.(\beta) = \sqrt{ \frac{s_e^2}{\sum(x_i - \bar x)^2} } \end{align}

と表せる（詳細略）。ここで$s_e^2$は、求めた回帰直線$(\hat Y_i=\beta X_i + \alpha)$で説明できない残差分散で、

\begin{align} s_e^2 = \frac{1}{N-2}\sum \bigl ( y_i - (\beta x_i + \alpha) \bigr ) ^2 \end{align}

である。

　母回帰係数を$\hat \beta$とすると、 \begin{align} t = \frac{\hat \beta - \beta}{s.e.(\beta)} \end{align} が自由度$N-2$の$t-$分布に従うことを利用して信頼区間を求める。

すなわち、例えば$\hat \beta$の95%信頼区間は \begin{align} \beta + t(N-2, 0.975)\cdot s.e.(\beta) < \hat \beta < \beta + t(N-2,0.025)\cdot s.e.(\beta) \end{align} として求まる。

◆ $t-$値の計算: [関数] scipy.stats.t()

  • [使い方]: rv = scipy.stats.t(df=): 自由度dfの$t-$分布をとるオブジェクト(rv)を作成.
  • そのオブジェクトについて、以下のメソッドを使って値を取り出す。
    • .isf(q,..): 上側確率がqとなる$t-$値を返す。
    • ...

（参考）「scipy.stats.t」について: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.t.html

【注意】自由度が十分に大きいときには、$t(0.025) \sim 2$、$t(0.975) \sim -2$であるの。この値は覚えておくと良い。

◆ 自由度について

　自由度をいくつに取るかというのはなかなか悩ましい問題である。気象データの変動は、時間的にも空間的にも相関を 持っているのでこれらを考慮して自由度を決める必要があるからである。

　ここでの例で言えば、解析に用いたデータの数（〜30日ｘ年数）を自由度とすると、明らかに過大評価である。 なぜなら、昨日の気温と今日の気温は相関を持っているからである。そこでここでは、まずは各年１個のデータを 作り、各々のデータが独立であると仮定してデータ年数を自由度としている （もちろん年々変動もあるのでこれが正しいとは言い切れないが）。

　他に自由度の見積もり方としては例えば、自己相関係数を計算して、それらが小さくなる（0.2-0.3)ラグ期間を 特徴的な時間スケールとして、データの長さをこれらのスケールで割ったものとする手法がある（松山・谷本、2005)。

【実践】以上の準備をしたうえで2.2節で算出したトレンドの信頼区間を推定してみよう

In [3]:
t_annual = temp.resample('A-DEC',closed='right').mean()

# トータル時間数を示す値の配列をつくる--------
x = t_annual.index # 時刻インデックスを取得
x0 = x[0] # 基準時刻
dt_index = (x - x0).days # 日数に変換
dt = dt_index.values #日数の値を取得->x_iの配列
#--------------------------------------

# 関数を定義する
def f_trend(x, a, b):
    return a*x + b

popt,pcov =  curve_fit(f_trend,dt,t_annual)
# popt[0]:a(回帰直線の傾き)、popt[1]:b(y切片) 

fitting = popt[0]*dt + popt[1]

# トレンド(10年あたりで表す)
trend10 = popt[0]*10*365.25
print('Trend = ',trend10,'K/decade')

# 以下、信頼区間の計算
N = yen-yst+1 #年数

# 標準誤差を計算
se2 = ((t_annual - fitting)**2).sum()/(N-2) #残差分散
x2 = dt.var()*N
error = np.sqrt(se2/x2)*365.25*10

# トレンド（回帰直線の傾き）の95%信頼区間
rv = stats.t(df=N-2) # N:年数
print('t_value(0.025) =', rv.isf(0.025)) #t-値 (Nが十分に大きければ～２) 

print('信頼区間: ',trend10-error*rv.isf(0.025),'< trend <', trend10+error*rv.isf(0.025),'(degrees/decade)')

Trend =  0.35585801503769093 K/decade
t_value(0.025) = 2.0261924630291097
信頼区間:  0.21828239626705082 < trend < 0.49343363380833105 (degrees/decade)

95%信頼区間を考慮してもトレンドは正であり、統計的に有意な上昇トレンドを持っていることがわかる。

C.2. コンポジット解析の信頼区間推定（2.3節の補足）

　コンポジット解析の結果には、$t-$検定を用いた信頼区間推定を行うことが多い（母平均の区間推定）。 推定した平均値の標準誤差(s.e.)は、標本標準偏差($s$)を用いて

\begin{align} s.e. = \frac{s}{\sqrt{n-1}} \end{align}

と表せる。ここで$n$はサンプルの自由度を表す。

　 母平均を$\mu$、標本平均を$\bar x$とすると、$t-$値は \begin{align} t_{\bar x} = \frac{\bar x - \mu}{s.e} \end{align} と表せる。自由度が$n$のとき、これが自由度$n-1$の$t-$分布に従うことが知られている。したがって、例えば母平均$\mu$の95%信頼区間は

\begin{align} t(n-1,0.975) < & t_{\bar x}<t(n-1,0.025) \\ \Longleftrightarrow \ \bar x - t(0.025)\frac{s}{\sqrt{n-1}} < & \mu < \bar x - t(0.975)\frac{s}{\sqrt{n-1}} \end{align}

と求まる。

◆「標準偏差」か「標準誤差（信頼区間）」か？

標準偏差(Standard Deviation)と標準誤差(Standard Error)はしばしば混同しがちである。

  • 標準偏差: 平均のまわりのデータのバラツキ度合いを表す。データ数に依らない。

  • 標準誤差（信頼区間）: 求めた統計量の推定精度を表す指標であり検定にも用いる。データ数($n$)が多いほど小さくなる（つまり、推定精度は高くなる）。 　

　エラーバーで標準誤差を描くか標準誤差を描くかは、実際にそのグラフで何を表現したいかによる。 推定した平均値の精度を表したい場合（例：ある値と有意な差があるかどうかを調べたい場合）には標準誤差を示せば良いし、 平均の周りでデータがどれほどバラついているかを表したい場合には標準偏差を示せば良い。

　以下の例で言えば、標準誤差を描けば「１月の平均値と２月の平均値で有意な差があるか？」を議論できるし、 標準偏差を描けば「各月において気温の変動幅がどの程度か？」を表せる。

【実践】2.3節で算出した月平均気温にエラーバーをつけてみよう

In [4]:
# 月平均データを作成：
# 日平均データをそのまま使うのではなく、各年１個の月平均データを予め作る
#　--> 上記の「自由度について」を参照。
temp_m = temp.resample('M',closed='right').mean()

key = temp_m.index.month
group = temp_m.groupby(key)

# コンポジット値（これは2.3で求めた）
data_MM = group.mean()

### 以下エラーバーの計算----
# コンポジット値まわりの標準偏差を計算（->エラーバー計算のため）
std = group.std()

# 検定用
df = yen-yst+1 #自由度: degrees of freedom（データの標本数=重ねる年数）
rv = stats.t(df=df-1)

## 95%信頼区間の計算
error = std/np.sqrt(df-1)*rv.isf(0.025)
#---


#以下、描画
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)

x = np.arange(1,13)
ax.errorbar(x, data_MM, error, fmt='ro', capsize=4, ecolor='black',label='Monthly Average + 95% level')
ax.errorbar(x+0.15, data_MM, std, fmt='bo', capsize=4, ecolor='blue',label='Monthly Average + std')

ax.set_ylabel('Temperature ($^\circ$C)')
ax.set_xlabel('Month')
ax.legend()

plt.show()

  • 上のグラフから、以下のことが言えそうである：

    • 「2月と3月の気温は有意に異なる」（エラーバーが重なっていない）
    • 「1月と2月の気温差は有意とは言えない」（エラーバーが重なっている）

  • データの期間を変えてみて、データ数が増えるに従ってエラーバーは小さくなることを確かめると良い。

〇 エラーバー付きのマーカープロット: ax.errorbar (x, y, xerr =, yerr=, capsize=, ecolor=)

[オプションの説明]

  • x : x軸の値（x軸のデータ点；サイズ$N$)
  • y : y軸の値（y軸のデータ点; サイズ$N$）
  • yerr : y軸方向のエラーバーの範囲を与える。一次元配列(N,)の時、各データ点からの距離（$+/-$に等距離）を与える。二次元配列(N,2)の時、エラーバーの下限値と上限値を別々に与える。
  • xerr: 同上（x軸方向のエラーバー）
  • capsize: エラーバーの"帽子"の長さ
  • ecolor: エラーバーの色
  • ...