2. 時系列データの解析（統計解析・時系列解析）

2.1 基本的な統計量の算出

2.2 ダウンサンプリングとトレンド推定

2.3 コンポジット解析

2.4 【参考】移動平均（時間フィルターの基礎）

　前章では、ファイルの読み込みと基本的な描画を紹介した。 続いて本章では、長期時系列データセットを用いた基礎的な時系列解析・統計解析を実習する。

　本演習では、主にpandasのメソッドを解析に用いることにする。

　ただし、インデックスがタグ付けされていないデータ配列(ndarray)を用いた解析例(Numpyの関数を用いた手法）についても 参考までに併記している。なお、Numpy関数を用いる場合にでもpandasのオブジェクトを引数としてそのまま渡せる場合が多い。

【準備】 複数のファイルを結合して長期時系列データを作る（1.1と同じ）

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 = 1960; yen = 2020
dlist = []
for yy in range(yst,yen+1,1): # range(start,stop,step)は "start"から"stop"（ただし,stopは含まない）まで増分"step"の等差数列を作る
    cyr = str(yy) # 文字列化
    file = "../../DB_2020/data/Kyoto/Kyoto_" + cyr  +".csv" #<---【注】"data/AMeDAS/Kyoto/Kyoto_"にパスを変更のこと
    
    df = pd.read_csv(file,encoding = "shift-jis",header=2,index_col=0,parse_dates=[0])
    # "encoding=":文字コードを指定、"header=2":2行目をヘッダー行（変数情報）に指定、"index_col=0": 0列目をインデックス行（時刻）に指定
    # ""parse_dates =[0]": 「日時」の型として読み込む（これにより、年・月・曜日といった情報を簡単に取り出せる） 
    
    df = df.iloc[2:,:] # ヘッダーの余計な部分（空白行など）を除いて上書き
    # これで一年分のデータのかたまりができた

    dlist.append(df) # listに要素(df)を追加していくことで

# pandas.concatで、データのかたまり（41年分）を結合する
dfm = pd.concat(dlist,axis=0) # axis=0 (時間軸方向）に結合
    
temp = dfm['平均気温(℃)']

2.1 基本的な統計量の算出・相関係数

◆ 統計量を一括して表示する: describe()　メソッド

Series、DataFrame型のデータについてdescribe()メソッドを使うと、基本的な統計量を一覧にして示してくれる。

◆ 個別の統計量を計算する

pandasでは個々の統計量を計算するメソッドも多種用意されている。下表に代表例を示す：

pandasメソッド 説明
要素の数 count
sum 合計
mean 平均
median 中央値
quantile 分位値（0から1の範囲）
min, max 最小値，最大値
idxmin, idxmax 最小値、最大値が得られた要素のラベル（注１）
var(ddof=) 分散（ddof = 0: 標本分散, 1: 不偏分散,...; デフォルトは1）
std(ddof=) 標準偏差（ddof = 0, 1,..)
cov(other) otherとの共分散

（注） 複数のインデックス（日時）で同じ最小（大）値を示す場合には、"最も小さい"ラベル（日付で言えば"若い"ラベル）しか 結果を返さないので注意が必要である。これを確かめるには、data[data == data.min()]などとして（「1.2.1 条件を指定したデータの抜き出し」を参照）、最小値を取るデータを抜き出して調べる必要がある。

（参考） [Numpy]

Numpy においても、一部については同様のメソッドが用意されている(sum,mean,min,max,var,std など)。 その他に関しては、Numpyの関数を用いる：

Numpy関数 説明
np.median(array) 中央値
np.quantile(array,q) データのパーセント点を0から1の範囲(q)で求める

　なお統計量の計算にあたっては、pandasとは異なり欠損値は考慮されない。欠損値を考慮するにはnansum、nanmean などのメソッドを別途用いる。

参考） ndarrayのメソッド

In [3]:
# 一括表示
print("describeの結果")
print(temp.describe())

describeの結果
count    22275.000000
mean        15.801688
std          8.648235
min         -3.400000
25%          7.800000
50%         16.100000
75%         23.300000
max         32.800000
Name: 平均気温(℃), dtype: float64
In [4]:
# 各統計量を個別に計算
# 平均
mean = temp.mean()
print(mean)

15.801687991021284

【練習課題】2020年に猛暑日となった日数は？

◆ 相関係数: [pandas] .corr メソッド, [Numpy] np.corrcoef(x,y)

  • [pandas] .corr(other): 各行同士の相関係数を計算し、DataFrame形式で戻してくれる。Series形式の場合は引数に相手(other)が必要
  • [Numpy] np.corrcoef(x,y): x, yの相関係数を計算。結果は、相関係数を要素に持つ行列として返される。

相関係数の定義

\begin{align} R(x:y) = \frac{\sum (x_i - \bar x)(y_i - \bar y)}{ \sqrt{\sum (x_i-\bar x)^2} \sqrt{\sum (y_i - \bar y)^2} } \end{align}
In [5]:
precip = dfm['降水量の合計(mm)']

corr_tp = temp.corr(precip)
corr_tp_np = np.corrcoef(temp,precip)

print('R(pandas) ',corr_tp)
print('R(np)', corr_tp_np)

R(pandas)  0.12218649142960808
R(np) [[nan nan]
 [nan  1.]]

【練習課題】相関係数の定義に沿った方法で計算し、上記の結果と一致することを確かめよ。

【課題 2-1】データの期間内において、最高（最低）気温の最大（小）値を記録したのは何年何月何日か？（＋２点）

【課題 2-2】日降水量の頻度分布を描いてみよ。併せて、「10年に一度の発生頻度の雨」とはどの程度か求めてみよ（＋２点）

（ヒント）分位値メソッド

【課題 2-3】 6-8月の日最高気温について、異なる２つの期間（例：1979-1998と1999-2018）のヒストグラムを重ね合わせることで、極端気象の発現度に違いがあるか調べよ（＋２点）

（Tips1）異なるヒストグラムを重ねる際には、ヒストグラム描画のオプションでalpha=0.5（グラフの透過度）などとすると見やすい。 （Tips2）特定の月のデータの取り出しは2.3節も参照のこと

2.2 ダウンサンプリング と トレンド推定

　気象場の変動は、様々な時間スケールの変動が重なって生じている（2.5節「スペクトル解析」も参照）。 データ解析では、予め着目する時間スケールの変動を抽出して解析を行うことも多い。

　例えば年々〜十年規模変動を見たい場合には年平均のデータがあれば十分だし、 逆に低気圧・高気圧の通過に伴うような数日周期の変動を見たい場合には長い時間スケールの変動を除いた方が見やすい。

　以下では、アメダスの気温データを用いて年平均データを作成した上でトレンドを推定してみよう。

◆ データの間隔を変える: resample メソッド

  • [pandas メソッド] .resample(freq,axis=,closed=):

    • freq: 時系列の基準頻度を表す（"A-DEC", "M", "3D"など）。下表参照。基準頻度の前に整数を書けば、倍数が作られる（e.g.,"3D": 3日）。
    • axis=: 対象とする軸。デフォルトはaxis=0
    • closed=: ダウンサンプリングのとき、どちらの端のデータを含めるか。right(left):右端（左端）を含める。
文字 説明
A-JAN, A-FEB,... 一年に１度、指定した月の月末ごと
M 月末ごと
D 日ごと
H 毎時
T (or min) 毎分
S 毎秒

（余談）[Numpy]

　pandasを用いずに同様の処理を行うのは少々面倒である。例えば年平均のデータを作るにしてみても、 各年の日数分（閏年を考慮）を平均するといった操作が必要になる。

いっそのこと、用いるデータが何であれ、pandasの型(1次元:Series, 2次元: DateFrame）にしてしまうのが 一番楽なのかも知れない。

◆ 曲線/直線でフィッティング：[関数] scipy.optimize.curvefit

  • 使い方: popt, pcov = scipy.optimize.curvefit(f, xdata, ydata, ....):

    • 関数 ydata = f(xdata)をフィッティングする。関数fは予め自分で定義する（下記の例を参照）

    • 出力値： popt: フィッティングパラメター と pcov: 共分散

【練習課題】年平均のデータを作り時系列で図示せよ。またこのデータに一次関数をフィッティングすることで、トレンドを推定せよ。

In [6]:
## (1) [pandas] resampleメソッド
t_annual = temp.resample('A-DEC').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]
fitting = f_trend(dt,popt[0],popt[1])

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

# 描画
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)
ax.plot(x,t_annual,label='Annual-mean T')
ax.plot(x,fitting,label='Regression')
ax.legend()
plt.show()

Trend =  0.2756897475483593 K/decade

◆ 時刻型のインデックスから、積算時間への変換

これまではx軸の座標（インデックス）を日時(DateTime型；例,"2017-01-01"）として扱ってきたが、 ここでは回帰係数を計算するために日時を開始時刻からの積算時間とする必要がある。 pandasでは便利なことに時刻の差分を計算できる(Timedelta型という）ので、これを用いる。 具体的には、以下の手順を踏んでいる：

  1. 時刻(Datetime型）の差分を取る。これにより、Timedelta型のオブジェクト(dt_index)を作成。

  2. dt_indexの属性である.day（日数）、.seconds（秒数）などを取り出す。

  3. このままでは(pandasの)indexとして見なされているので、最後に.valuesメソッドを使って数値配列(dt)にする（これが$x_i$の配列)。

（参考）TimedeltaIndex型の属性・メソッド一覧はこちら

◆ 信頼区間の推定（詳細はこちらを参照）

本演習では省略するが、実際には求めたトレンドの信頼区間（エラーバー）を推定する必要がある。その幅（誤差範囲）を考慮した上でもトレンドが正で あれば、「統計的に有意な上昇トレンドがある」と（自信を持って）言えることになる。

【課題 2-4】 回帰係数を以下の定義に従って直接計算し、上記の結果に一致することを確かめよ。和の計算にはメソッドを用いて良い（＋２点）

　 $$ \beta = \frac{\sum_{i=1}^N (x_i-\bar x)(y_i - \bar y)}{\sum_{i=1}^N (x_i - \bar x)^2} $$

【課題 2-5】 真夏日/猛暑日の年間日数のトレンドが分かるグラフを作成せよ（＋３点）

　

2.3 コンポジット解析　〜月平均値を作ってみる〜

　ある条件を満たすデータのみを抽出し、その条件下での平均的な特徴を求めることをコンポジット（合成図）解析という（「重ね合わせる」というニュアンス）。気象データ解析においては、ある事象の平均的な構造を調べたり、事象の振る舞いが 外部条件によってどのように変わるかを調べる際によく使われる手法である。

　簡単な例としては：

1）１月の平均気温を求める（条件：「月が１月であること」）

2）気圧が高い日と低い日で気温は有意に異なるか調べる（条件：「気圧がXX以上」「気圧がXX以下」でそれぞれ平均を取って比べる）

などが挙げられる。

【練習課題】 月平均気温の年サイクルを求めよ。

A. Groupby メソッドを使った方法　

　ある条件・基準（「キー」と呼ぶ）にしたがってデータをグループ化し、それぞれのグループごとに関数(sum, meanなど）を適用した 結果を返すことができる。

　キーとしては

  • 値のリストや配列（グループされる軸と同じ長さのもの）
  • データフレームの列名を示す値

  • 軸のインデックス、あるいは、各インデックスのラベルに対して呼び出される関数

    などを用いることができる。

　ここでは、インデックス配列がdatetime型（日付・時刻）であるので、その属性を利用して 「月」のインデックス配列を作成し、それをキーとして与えることでグループ化する（上の例の一番目）。

DatetimeIndexの属性としては、年・月・日などがある（以下、一例）：

属性 意味
year 年
month 月
day 日
hour 時間
dayofweek 曜日（0:月曜日...6:日曜日）
dayofyear 一年の何日目か

（参考）DatetimeIndexの属性・メソッドの一覧はこちら

ここではmonth属性を取得して、各インデックスの月を求める。

In [7]:
# 時間情報を元に条件を決めるのであるから、まずは時間の情報を抜き出す。
date = temp.index #ndex列が時間情報なので
#print(date)

# "月"の配列を作る
key = date.month
#print(key)

# これをキーとしてグループ化
group = temp.groupby(key)

# 各グループで平均を計算
temp_monthly = group.mean()

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

temp_monthly.plot(ax=ax,label="Monthly mean")

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

plt.show()

B. 真偽値関数を使う方法

　より汎用性の高い方法として、1.2.2節で紹介した真偽値型(論理型; bool)の配列を使う方法もあるので紹介しておく。

(Step 1)１月の平均気温を計算してみる

　まずは1月の平均気温を計算してみよう。

　「１月ならTrue、それ以外はFalse」の真偽値配列を作って１月のデータを取り出す。

In [8]:
#真偽値配列の作成
mask_mon = (temp.index.month == 1)

#１月データを抜き出して平均
temp_jan = temp[mask_mon].mean()

print(temp_jan)

4.4226864093072455

（Step 2) Step1 を応用し、１−１２月の各月について月平均を算出して図示する

各月において上記と同様の処理を繰り返せば良い。

In [9]:
mon_ave = np.zeros(12) #大きさ12（12か月分）の配列を用意しておく
for mm in range(1,13):
    mon_ave[mm-1] = temp[temp.index.month == mm].mean()
    # mm-1 としたのは配列の要素が0から始まるため
    
x = np.arange(1,13) #x軸（1-12の月を表す整数配列）

## ここから描画
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)

ax.plot(x,mon_ave,label='Monthly Average')
ax.set_ylabel('Temperature ($^\circ$C)')
ax.set_xlabel('Month')

ax.legend()
plt.show()

◆ 信頼区間の推定（詳細はこちらを参照）

トレンドの解析と同様、コンポジットした値にも信頼区間（エラーバー）をつけるのがマナーである。詳細は上記リンクを参照のこと。

【課題 2-6】上記の結果に年周期の三角関数をフィッティングすることで、季節変化の振幅を求めよ（＋2点）

$\ \ \ $ヒント：三角関数はnp.sin()、np.cos()で与えられる。また引数の単位はラジアンであることに注意。

【課題 2-7】曜日による気温の変動はあるだろうか？曜日によるコンポジット図を作成してみよ [ 参考：藤部 (2014) ]　（＋2点）

【課題 2-8】１時間毎のデータを使って、気温と気圧の平均的な日変化を求めてみよ。京都と富士山のデータを比較すると面白い（＋2点）

$\ \ $参考：１時間毎のAMeDASデータの読み込みはこちら（の最初のセル）を参考にしてください

2.5【参考】移動平均・時間フィルター: [pandas] rolling メソッド, [np] convolve 関数

　移動平均 (Running-mean)とは、一言で言えばデータを均（なら）す処理である。式で書けば、 \begin{equation} y_j = \frac{1}{2K+1}\sum_{k=-K}^{K} x_{j+k} \end{equation} ここで、$x_j$は元の時系列、$y_j$は移動平均した時系列。すなわち、各点$j$において、その周り$\pm K$点の平均を計算することを意味する。

　長周期変動を抜き出すという意味で、移動平均は簡単なローパスフィルターとしても使われる（さらに、これを元の時系列から差し引けば短周期成分を抜き出すハイパスフィルターとなる）。但し、そのフィルター特性は必ずしも良くないので、実際の研究で用いる際には注意を要する。大まかに言って『$n$日移動平均は、ほぼ$2n$日以下の周期を落とすフィルターである』（伊藤・見延, 2010)。（より正確には、移動平均した時系列と元の時系列のパワー（振幅の二乗）を比べると、周期$2n$でパワーが半分になる（周期$2n$以下ではそれ以下））。

（補足）より特性の良いフィルターとしては、「Lanczos（ランチョス）フィルター」などが良く用いられる。

  • [pandas] .rolling(wd, center=, win_type =,...): 移動する窓関数。窓の中にあるデータをまとめてしまう。.mean()などと一緒に用いる。デフォルトではデータの端($<$wd/2, n-wd/2$<$)は欠損(NAN)となる。
    • wd: 窓の長さ（要素数；上の式での$2K+1$）
    • center: True なら各indexの前後合わせて wd 個の要素を使う（デフォルトはFalse)
    • win_type:窓の種類（デフォルトは等重み；他にhammingなど）
  • (参考）[Numpy] np.convolve(a,v, mode=): 畳み込み積分（定義は以下）。a: 時系列データ、b:重み関数。実際には、長い方の配列の要素１つ１つに対して短い方の配列の要素全てを掛け合わせたものを足し合わせるという操作を行っている。mode:valid(端っこは計算しない）、same（長い方の配列と同じ）
\begin{align} (a*v)[n] = \sum_{m=-\infty}^{\infty} a[m] v[n-m] \end{align}

　　例) $v=[0.33,0.33,0.33]$の時（三点移動平均）のときは（$0 \le n-m \le 2$のみ計算することに注意すると)、

\begin{align} (a*v)[0] &= a[0]*v[0] \\ (a*v)[1] &= a[1]*v[0] + a[0]*v[1] \\ (a*v)[2] &= a[2]*v[0] + a[1]*v[1] + a[0]*v[2] \\ (a*v)[3] &= a[3]*v[0] + a[2]*v[1] + a[1]*v[2] \\ \cdots \end{align}

このうち、mode=validでは、$(a*v)[2]$からの値のみを返す（すなわち元の配列とは大きさが異なる）。 したがって図示する際には、x軸（時間）も最初と最後をトリミングしたものを用いる必要があるので注意。

In [10]:
a = [1,2,3,4,5]
v = [0.33,0.33,0.33]
print(np.convolve(a,v,mode="valid"))

[1.98 2.97 3.96]

【練習課題】31日移動平均を計算し、元のデータ時系列と併せて図示せよ。

In [11]:
# (1) [pandas] rolling を使った方法 -----
wd = 31
t_mv1 = temp.rolling(wd,center=True).mean()
x1 = t_mv1.index # t_mv1のインデックスをそのまま使う
#print(x1)
# ------------------------------

# (2) [Nump] convolveを使った方法 -----
temp_v = temp.values # 値を取り出す
#print(type(temp_v))

# 重み関数の定義 
w = np.ones(wd)/wd 
t_mv2 = np.convolve(temp_v,w,mode='valid')

# x軸（時間軸）をトリミング
wd2 = int((wd-1)/2)
x2 = temp.index[wd2:-wd2]
print(x2)
#----------------------------

# 図示
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)

#
ax.plot(x1[:90], temp[:90],label='(0) Original')
ax.plot(x1[:90], t_mv1[:90],label='(1) 31-day moving average (rolling)')
#ax.plot(x2[:90], t_mv2[:90],label='(2) Convolution')

ax.legend()
plt.show()

DatetimeIndex(['1960-01-16', '1960-01-17', '1960-01-18', '1960-01-19',
               '1960-01-20', '1960-01-21', '1960-01-22', '1960-01-23',
               '1960-01-24', '1960-01-25',
               ...
               '2020-12-07', '2020-12-08', '2020-12-09', '2020-12-10',
               '2020-12-11', '2020-12-12', '2020-12-13', '2020-12-14',
               '2020-12-15', '2020-12-16'],
              dtype='datetime64[ns]', name='年月日', length=22251, freq=None)

これを元の時系列から差し引けば、ハイパスフィルターとなる。

In [12]:
### High-pass filter

# オリジナル時系列から低周波成分を差し引く
t_high = temp - t_mv1

#図示
fig = plt.figure(figsize=(10,4))
ax = fig.add_subplot(111)
ax.plot(t_high[:90],label='High-pass filtered')
ax.legend()
plt.show()