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

2.1 基本的な統計量の算出

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

2.3 コンポジット解析

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

---

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

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

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

　データ解析においては、データの読み込み と 解析で、ライブラリを使い分けることも多い。 例えばcsvファイルであればpandasで読み込むのが楽だが、解析に使えるメソッドは基本的な統計解析に限られる。 そこでファイルの読み込みはpandasで行い、データを数値配列(ndarray)に変換した上で、汎用性の高いNumpyやScipyを 使って詳細な解析行う、といった使い方をする（演習ではそこまで複雑な解析をしないのでその必要はないが）。 ともかく、データの型によって使えるメソッドが限られるという点には十分注意する必要がある。

---

【準備】 複数のファイルを結合して長期時系列データを作る（1.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

# 定数設定

def read_csv(START_YEAR,END_YEAR,STATION):
#START_YEAR = 1960
#END_YEAR = 2023
    DATA_DIR = "../../DB_2020/data/"+STATION+"/"+STATION+"_" #<---【注】"data/AMeDAS/Kyoto/Kyoto_"に変更のこと
    ENCODING = "shift-jis"

    # データを格納するリスト
    dataframes = []

    for year in range(START_YEAR, END_YEAR + 1):
        file_path = DATA_DIR + str(year) + ".csv"

        df = pd.read_csv(
            file_path,
            encoding=ENCODING,
            header=2,
            index_col=0,
            parse_dates=[0]
        )

        # 余計なヘッダー行を除去
        df = df.iloc[2:, :]

        dataframes.append(df)

    # すべての年のデータを結合（時間軸方向に）
    dfm = pd.concat(dataframes, axis=0)
    
    return dfm
    

dfm_kyoto = read_csv(1961,2023,"Kyoto")
dfm_maizuru = read_csv(1961,2023,"Maizuru")

temp_kyoto = dfm_kyoto['平均気温(℃)']
temp_maizuru = dfm_maizuru['平均気温(℃)']

---

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(ndarray)	中央値
np.quantile(ndarray,q)	データのパーセント点を0から1の範囲(q)で求める

　なお統計量の計算にあたっては、pandasとは異なり欠損値は考慮されない。欠損値を考慮するにはnp.nansum()、nanmean() などの関数を別途用いる。

参考） ndarrayのメソッド

# 一括表示
print("describeの結果")
print(temp_kyoto.describe())

describeの結果
count    23004.000000
mean        15.871835
std          8.657889
min         -3.400000
25%          7.900000
50%         16.200000
75%         23.400000
max         32.800000
Name: 平均気温(℃), dtype: float64

# 平均の計算
# (1) メソッドを使う
mean = temp_kyoto.mean()
print(mean)

15.871835332985546

# (2) pandasで読み込んだデータを ndarray（数値配列）に変換し、numpyの関数を使って計算
temp_kyoto_v = np.array(temp_kyoto.values)
np.nanmean(temp_kyoto_v)

15.871835332985569

# (3) ループを回して地道に計算
t_sum = 0
i = 0
for t in temp_kyoto_v:
    if(t < 999): #異常値処理
        t_sum += t
        i +=1
    else:
        print(t)

print(t_sum/i)

nan
nan
nan
nan
nan
nan
15.871835332985546

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

---

◆ 相関係数: [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}

# [1] 相関係数の計算（pandas）
corr_t = temp_kyoto.corr(temp_maizuru)

# [2] 欠損値を考慮した NumPy による相関係数の計算
mask1 = ~temp_kyoto.isnull()  # 値が存在する箇所を True に： "np.logical_not()"と同義
mask2 = ~temp_maizuru.isnull()
mask = mask1 & mask2          # 両方に値がある箇所のみ True

# 欠損値を除いたデータで相関係数を計算
corr_t_np = np.corrcoef(temp_kyoto[mask], temp_maizuru[mask])

# 散布図の描画
fig, ax = plt.subplots()
ax.scatter(temp_kyoto, temp_maizuru)
plt.xlabel("Kyoto Temperature")
plt.ylabel("Maizuru Temperature")
plt.title("Temperature Correlation between Kyoto and Maizuru")
plt.grid(True)
plt.show()

# 結果出力
print('相関係数（pandas）:', corr_t)
print('相関係数（NumPy） :', corr_t_np[0, 1])

相関係数（pandas）: 0.991306351269801
相関係数（NumPy） : 0.991306351269801

相関があるからと言って、因果関係があるとは言えません。いわゆる疑似相関に注意しましょう。

---

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

---

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

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

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

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

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

---

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

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

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

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

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

---

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

  • freq: 時系列の基準頻度を表す（"A", "M", "3D"など）。下表参照。基準頻度の前に整数を書けば、倍数が作られる（e.g.,"3D": 3日）。
  • axis=: 対象とする軸。デフォルトはaxis=0
  • closed=: ダウンサンプリングのとき、どちらの端のデータを含めるか。right(left):右端（左端）を含める。

文字	説明
A	毎年
M	毎月
D	毎日
H	毎時
T (or min)	毎分
S	毎秒

AE-JAN, AE-FEB,... | 一年に１度、指定した月の月末ごと

（余談）[Numpy]

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

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

---

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

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

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

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

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

# 年平均気温データの作成
t_annual = temp_kyoto.resample('A').mean()
yr = t_annual.index.year

# 線形回帰関数の定義
def f_trend(x, a, b):
    return a * x + b

# 回帰直線のフィッティング
popt, pcov = curve_fit(f_trend, yr, t_annual)
fitting = f_trend(yr, *popt)  # popt[0]: 傾き, popt[1]: 切片

# トレンドを10年あたりで表現
trend10 = popt[0] * 10
print(f'Trend = {trend10:.4f} K/decade')

# プロット描画
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(yr, t_annual, label='Annual-mean Temperature')
ax.plot(yr, fitting, label='Linear Trend', linestyle='--', color='red')
ax.set_xlabel('Year')
ax.set_ylabel('Temperature (°C)')
ax.set_title('Annual Mean Temperature and Linear Trend')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

Trend = 0.2976 K/decade

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

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

---

【課題 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属性を取得して、各インデックスの月を求める。

# 月ごとの平均気温を算出
key = temp_kyoto.index.month
group = temp_kyoto.groupby(key)
temp_monthly = group.mean()

# プロット描画
fig, ax = plt.subplots(figsize=(10, 4))
temp_monthly.plot(ax=ax, label="Monthly Mean", marker='o')

# 軸ラベル・凡例など設定
ax.set_xlabel('Month')
ax.set_ylabel('Temperature (°C)')
ax.set_title('Monthly Mean Temperature')
ax.set_xticks(range(1, 13))
ax.legend()
ax.grid(True)


plt.show()

---

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

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

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

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

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

#真偽値配列の作成
mask_mon = (temp_kyoto.index.month == 1)

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

print(temp_kyoto_jan)

4.448796722990274

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

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

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

## ここから描画

# プロット描画
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(x,mon_ave,label='Monthly Mean', marker='o')

# 軸ラベル・凡例など設定
ax.set_xlabel('Month')
ax.set_ylabel('Temperature (°C)')
ax.set_title('Monthly Mean Temperature')
ax.set_xticks(range(1, 13))
ax.legend()
ax.grid(True)

plt.show()

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

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

---

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

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

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

【課題 2-8】6-8月の期間において、気圧が高い日/低い日で気温がどれほど違うか計算してみよ。ここで、気圧が「平均+1σ（σ：標準偏差）以上の日」を「気圧が高い日」、「平均-1σ以下の日」を「気圧が低い日」とする（＋2点）

【課題 2-9】１時間毎のデータを使って、気温と気圧の平均的な日変化を求め、結果について簡単に考察せよ。京都と富士山のデータを比較すると面白い（＋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軸（時間）も最初と最後をトリミングしたものを用いる必要があるので注意。

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日移動平均を計算し、元のデータ時系列と併せて図示せよ。

# 移動平均ウィンドウ（31日）
wd = 31

# (1) [pandas] rolling を使った方法
t_mv1 = temp_kyoto.rolling(wd, center=True).mean()
x1 = t_mv1.index
print(x1)

# (2) [NumPy] convolve を使った方法
temp_kyoto_v = temp_kyoto.values  # 値のみ抽出
w = np.ones(wd) / wd        # 重み関数（平均）
t_mv2 = np.convolve(temp_kyoto_v, w, mode='valid')  # 畳み込みで移動平均を計算

# x軸（時間）をトリミング
#wd2 = wd // 2
wd2 = int((wd-1)/2)
x2 = temp_kyoto.index[wd2:-wd2]  # mode='valid' に合わせる
#print(x2)  # 必要なら確認

# プロット描画（最初の90日間のみ表示）
fig, ax = plt.subplots(figsize=(10, 4))

ax.plot(x1[:90], temp_kyoto[:90], label='(0) Original')
ax.plot(x1[:90], t_mv1[:90], label='(1) 31-day MA (rolling)')
ax.plot(x2[:90], t_mv2[:90], label='(2) 31-day MA (convolve)', linestyle='--')

# 軸・凡例・装飾
ax.set_xlabel('Date')
ax.set_ylabel('Temperature (°C)')
ax.set_title('31-Day Moving Average: pandas vs NumPy')
ax.legend()
ax.grid(True)

plt.tight_layout()
plt.show()

DatetimeIndex(['1961-01-01', '1961-01-02', '1961-01-03', '1961-01-04',
               '1961-01-05', '1961-01-06', '1961-01-07', '1961-01-08',
               '1961-01-09', '1961-01-10',
               ...
               '2023-12-22', '2023-12-23', '2023-12-24', '2023-12-25',
               '2023-12-26', '2023-12-27', '2023-12-28', '2023-12-29',
               '2023-12-30', '2023-12-31'],
              dtype='datetime64[ns]', name='年月日', length=23010, freq=None)

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

### High-pass filter

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

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