1. 背景：相関はどの時間スケールを見ているのか
前回までに、Niño3.4 海面水温偏差と SOI の間には強い負の相関があることを確認した。Niño3.4 が高いとき SOI が低くなる傾向は、El Niño / La Niña に伴う大気海洋相互作用を反映している。
しかし、相関係数を 1 つだけ計算しても、その相関が月々の短周期変動によるものなのか、ENSO の 2〜7 年スケールによるものなのか、さらに長い変動によるものなのかは分からない。
r = cov(x,y) / (σx σy)
相関係数 r は 2 つの変数が一緒に増減する程度を表す。ただし、時系列では各月の値が完全に独立ではないため、見かけのデータ数 N をそのまま信用してはいけない。
この回では、FFT で時間スケールを分離し、相関係数だけでなく、有効自由度と p 値まで計算して「その相関をどの程度信じてよいか」を考える。
2. 自由度と p 値：相関を信じる前に
このページで検定したいのは、傾きではなく、Niño3.4 と SOI の相関係数である。傾きの検定は 09 の SST トレンドで詳しく扱う。07 では、相関係数 r が 0 と区別できるか、そして時系列では自由度をどう考えるべきかに集中する。
相関係数 r は何を見ているか
相関係数は、2 つの変数が平均から同じ向きにずれるか、反対向きにずれるかを表す量である。
r = Σ(xi − x̄)(yi − ȳ) / √{ Σ(xi − x̄)² Σ(yi − ȳ)² }
| 量 | 意味 |
|---|---|
| r > 0 | 片方が大きいとき、もう片方も大きくなりやすい。 |
| r < 0 | 片方が大きいとき、もう片方は小さくなりやすい。 |
| r ≈ 0 | 直線的な関係が弱い。 |
Niño3.4 と SOI では、El Niño 時に Niño3.4 が正、SOI が負になりやすいため、負の相関が期待される。
相関係数の検定で考える仮説
相関係数を計算しただけでは、その値が偶然出たものなのか、本当に関係があると考えてよいのかは分からない。そこで、次の帰無仮説を置く。
H0: ρ = 0
H1: ρ ≠ 0
H1: ρ ≠ 0
ここで ρ は母集団における本当の相関であり、r は手元のデータから計算した標本相関係数である。検定では、「本当は相関がない」と仮定したときに、観測された r ほど極端な値がどれくらい起こりうるかを調べる。
相関係数 r から t 値を作る
独立なデータが N 個あると仮定できる場合、相関係数の t 値は次で計算する。
t = r √{ (N − 2) / (1 − r²) }
この式は、「相関係数の大きさ」を「その不確かさ」で割った量と考えるとよい。|r| が大きいほど、また N が大きいほど、|t| は大きくなりやすい。
なぜ自由度は N − 2 なのか
ここは非常に重要である。単に「データが N 個あるから自由度も N」とは考えない。相関係数を計算する時点で、すでにデータから情報を使っているからである。
1. まず平均を引いている
相関係数では、元の値そのものではなく、平均との差を使う。
xi − x̄, yi − ȳ
平均を引いた量には、必ず次の制約がある。
Σ(xi − x̄) = 0
Σ(yi − ȳ) = 0
Σ(yi − ȳ) = 0
つまり、平均との差は N 個あっても、全部が自由に決まるわけではない。たとえば N−1 個の偏差が決まると、最後の 1 個は「合計が 0 になるように」自動的に決まる。
2. 相関係数は、標準化した変数どうしの回帰の傾きでもある
ここで、元の変数 x と y を、それぞれ平均 0、標準偏差 1 になるように標準化する。
Xi = (xi − x̄) / sx
Yi = (yi − ȳ) / sy
Yi = (yi − ȳ) / sy
標準化すると、単位や変動幅の違いが取り除かれる。つまり X も Y も、平均 0、標準偏差 1 の無次元量になる。
この標準化した X と Y に対して、原点を通る直線を考える。
Yi = c Xi + ei
ここで c は、この直線の傾きである。最小二乗法でこの傾き c を求めると、
c = ΣXiYi / ΣXi2
となる。これは、通常の回帰の傾き
傾き = X と Y の共変動 / X のばらつき
という形である。
一方、標準化した X は標準偏差が 1 なので、
ΣXi2 = N − 1
である。また、相関係数は標準化した変数を使うと、
r = ΣXiYi / (N − 1)
と書ける。したがって、
c = ΣXiYi / ΣXi2
= ΣXiYi / (N − 1)
= r
= ΣXiYi / (N − 1)
= r
ここで言いたいこと
相関係数 r は、ただの「関係の強さ」だけではなく、標準化した X と Y の間に直線をあてはめたときの傾きとしても見られる。
したがって、相関係数が 0 かどうかを調べることは、標準化した 2 変数の間に「直線的な傾きがあるかどうか」を調べることと対応している。
相関係数 r は、ただの「関係の強さ」だけではなく、標準化した X と Y の間に直線をあてはめたときの傾きとしても見られる。
したがって、相関係数が 0 かどうかを調べることは、標準化した 2 変数の間に「直線的な傾きがあるかどうか」を調べることと対応している。
3. 直線を決めると、残差には 2 つの制約が入る
相関を検定するときには、標準化した 2 変数の関係を、直線的な関係として評価している。直線をあてはめた後の残差を
ei = Yi − cXi
と書く。この残差は N 個ある。
e1, e2, e3, …, eN
何も条件がなければ、残差は N 個すべてが自由にばらつけるように見える。しかし、最小二乗法で直線を決めた後の残差は、勝手な値をとれるわけではない。必ず次の 2 つの条件を満たす。
Σei = 0
ΣXiei = 0
ΣXiei = 0
| 制約 | 意味 | 自由度への影響 |
|---|---|---|
| Σei = 0 | 残差全体が上側または下側に偏らない。 | 残差の合計が 0 になる必要があるので、最後の 1 個は自動的に決まる。 |
| ΣXiei = 0 | X が大きい側だけ、または小さい側だけに残差が偏らない。 | X との偏りも 0 になる必要があるので、さらに 1 個分の自由さがなくなる。 |
4. なぜ 2 つの制約があると N−2 になるのか
自由度とは、ざっくり言えば「条件を満たしながら、まだ自由に動ける値の数」である。
たとえば、残差が 5 個あるとする。
e1, e2, e3, e4, e5
もし何の制約もなければ、5 個すべてを自由に決められるので、自由度は 5 である。
しかし、1 つ目の制約
Σei = 0
があると、たとえば e1, e2, e3, e4 を決めた時点で、e5 は「合計が 0 になるように」自動的に決まる。したがって、自由度は 1 つ減る。
N → N − 1
さらに、2 つ目の制約
ΣXiei = 0
も同時に満たさなければならない。この条件は、残差が X の大きい側や小さい側に偏らないことを要求している。これにより、もう 1 個分の自由さが失われる。
N − 1 → N − 2
ここが N−2 の意味
残差は N 個あるが、
1つ目の制約 Σei = 0 によって 1 つ、
2つ目の制約 ΣXiei = 0 によってさらに 1 つ、
自由に動ける数が減る。したがって、相関係数の t 検定で使う自由度は N−2 になる。
残差は N 個あるが、
1つ目の制約 Σei = 0 によって 1 つ、
2つ目の制約 ΣXiei = 0 によってさらに 1 つ、
自由に動ける数が減る。したがって、相関係数の t 検定で使う自由度は N−2 になる。
注意
ここでの N−2 は、「データが 2 個消える」という意味ではない。
データは N 個あるが、直線関係を評価するために 2 つの条件を満たす必要があるため、独立にばらつける情報量が 2 個分減る、という意味である。
ここでの N−2 は、「データが 2 個消える」という意味ではない。
データは N 個あるが、直線関係を評価するために 2 つの条件を満たす必要があるため、独立にばらつける情報量が 2 個分減る、という意味である。
5. t 値の式はこう見ればよい
相関係数 r が 0 からどれくらい離れているかを見るには、r をその不確かさで割る。
t = 相関係数 / 相関係数の不確かさ
相関係数の不確かさは、独立なデータ数が増えるほど小さくなる。一方で、|r| が 1 に近づくほど、点は直線に近く並ぶので不確かさも小さくなる。この関係を整理すると、次の形になる。
t = r √{ (N − 2) / (1 − r²) }
学生向けに一言でいうと
相関係数の検定では、2 変数の平均的な位置と直線的な関係をデータから決めている。 そのため、見かけのデータ数 N から 2 つ分を差し引き、自由度を N−2 とする。
相関係数の検定では、2 変数の平均的な位置と直線的な関係をデータから決めている。 そのため、見かけのデータ数 N から 2 つ分を差し引き、自由度を N−2 とする。
t 値から p 値を求める
p 値は、帰無仮説 H0: ρ = 0 が正しいと仮定したとき、観測された |t| 以上に極端な値が出る確率である。両側検定では次のように計算する。
p = 2 × { 1 − Ft( |t| ; df ) }
ここで F_t は t 分布の累積分布関数である。Python では stats.t.cdf を使う。
t = r * np.sqrt((N - 2) / (1 - r**2))
p = 2 * (1 - stats.t.cdf(np.abs(t), df=N - 2))注意
p 値は「相関が存在する確率」ではない。
「本当は相関がない」と仮定したときに、今回のような相関が偶然に出る確率である。
p 値は「相関が存在する確率」ではない。
「本当は相関がない」と仮定したときに、今回のような相関が偶然に出る確率である。
時系列では N をそのまま使えない
ここからがこのページの重要点である。月ごとの Niño3.4 や SOI は、隣り合う月どうしが似ている。特にローパスやバンドパス後の時系列はなめらかになり、見かけのデータ数は多くても、実質的な独立情報量はかなり小さくなる。
たとえば 519 か月のデータがあっても、24–84か月のバンドパスをかけると、独立な情報は「519 個」あるわけではない。数年スケールの波を見ているため、実質的な波の個数はずっと少なくなる。
有効自由度 Neff
このページでは、lag-1 自己相関を使って有効自由度を近似する。
Neff = N × (1 − r1x r1y) / (1 + r1x r1y)
- r1x：Niño3.4 の lag-1 自己相関
- r1y：SOI の lag-1 自己相関
- Neff：実質的な独立データ数
自己相関が小さい場合、r1x r1y は小さくなり、Neff は N に近づく。一方、どちらの時系列も強く自己相関している場合、Neff は大きく減る。
r1x = np.corrcoef(x[:-1], x[1:])[0, 1]
r1y = np.corrcoef(y[:-1], y[1:])[0, 1]
Neff = N * (1 - r1x * r1y) / (1 + r1x * r1y)このページのコードでやっていること
したがって、このページでは普通の N ではなく、Neff を相関検定の式に入れる。
t = r √{ (Neff − 2) / (1 − r²) }
t = r * np.sqrt((Neff - 2) / (1 - r**2))
p = 2 * (1 - stats.t.cdf(np.abs(t), df=Neff - 2))07 の結論
フィルター後に相関係数の絶対値が大きくなっても、同時に自由度は小さくなる。したがって、フィルター後の相関は r だけでなく、Neff と p を一緒に見る。
フィルター後に相関係数の絶対値が大きくなっても、同時に自由度は小さくなる。したがって、フィルター後の相関は r だけでなく、Neff と p を一緒に見る。
3. クロススペクトルへ進むための周波数の考え方
クロススペクトル解析では、2 つの時系列が「どの周期帯で一緒に変動しているか」を見る。今回のページでは、まずその準備として、FFT により時系列を周波数成分へ分解し、24か月以上、24か月未満、24〜84か月という時間スケールに分ける。
| 処理 | 残す周期 | 意味 |
|---|---|---|
| ローパス | 24か月以上 | ゆっくりした変動を残す |
| ハイパス | 24か月未満 | 短周期変動を残す |
| バンドパス | 24〜84か月 | ENSO らしい 2〜7 年周期を残す |
周期と周波数は逆数の関係にある。周期 24か月は周波数 1/24 cycle/month、周期 84か月は 1/84 cycle/month である。
4. 使用データ
使用するデータは、これまでの回で使ったものと同じである。ファイルはこの HTML と同じディレクトリに置く。
| ファイル | 使う列 | 内容 |
|---|---|---|
| sstoi.indices | ANOM34 | Niño3.4 領域の海面水温偏差 |
| SOI_timeseries_updated.txt | SOI | 南方振動指数 |
5. 計算の流れ
- SOI と Niño3.4 を読み込む
- FFT を使って時系列を周波数領域へ変換する
- 24か月を境にローパス・ハイパスへ分ける
- 24–84か月だけを残すバンドパスフィルターを作る
- オリジナルとバンドパス後の相関係数を比較する
- lag-1 自己相関から有効自由度を計算する
- 有効自由度を使って p 値を計算する
6. Pythonでよく出る機能の意味
| 書き方 | 意味 |
|---|---|
| np.fft.fft(x) | 時系列を周波数成分へ分解する。 |
| np.fft.fftfreq(n, d=dt) | FFT の各成分に対応する周波数を作る。 |
| np.fft.ifft(X) | 周波数成分から時系列へ戻す。 |
| np.abs(freq) | 正負両側の周波数を同じように扱う。 |
| mask | 残す周波数を True、消す周波数を False とする配列。 |
| pd.merge(..., on="time") | 2つの時系列を共通の時刻でそろえる。 |
| np.corrcoef(x, y) | 相関係数を計算する。 |
| stats.t.cdf(...) | t 分布を使って p 値を計算する。 |
STEP 1
7. データの読み込み
ここでは、SOI と Niño3.4 を読み込む。06 より少し難しくするため、列名・読み込み条件・取り出す列も一部穴埋めにしている。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
# =========================================================
# 0. 設定
# =========================================================
soi_file = "SOI_timeseries_updated.txt"
nino_file = "sstoi.indices"
# カットオフ周期（月）
cutoff_period_months = ______ # 24か月
# 図の保存名
fig_soi_name = "______"
fig_nino_name = "______"
# =========================================================
# 1. データ読み込み
# =========================================================
# ---- SOI ----
soi_df = pd.read_csv(______, sep=______)
soi_time = pd.to_datetime(dict(year=soi_df[______], month=soi_df[______], day=15))
soi = soi_df[______].values.astype(float)
# ---- Niño3.4 ----
# sstoi.indices は ANOM 列名が重複するので、列名を自分で付ける
colnames = [
"YR", "MON",
"NINO12", "ANOM12",
"NINO3", "ANOM3",
"NINO4", "ANOM4",
"______", "______"
]
nino_df = pd.read_csv(
nino_file,
sep=r"\s+",
skiprows=______,
names=______
)
nino_time = pd.to_datetime(dict(year=nino_df[______], month=nino_df[______], day=15))
nino34 = nino_df[______].values.astype(float)
STEP 2
8. FFTフィルター関数
ローパス・ハイパス・バンドパスの本体である。ここでは、周期を周波数に変換して、残す周波数を mask で指定する。
# =========================================================
# 2. FFTフィルタ関数
# =========================================================
def fft_filter(x, dt=1.0, cutoff_period=24.0):
"""
x : 1次元時系列
dt : サンプリング間隔（月）
cutoff_period : カットオフ周期（月）
return:
x_low : ローパス成分（period >= cutoff_period）
x_high : ハイパス成分（period < cutoff_period）
freq : 周波数 [cycle/month]
X : FFT係数
"""
x = np.asarray(______, dtype=float)
# 欠損があれば線形補間
if np.any(~np.isfinite(______)):
s = pd.Series(x)
x = s.interpolate(limit_direction="both").values
n = len(______)
# 平均を別に保持
x_mean = np.mean(______)
x0 = ______ - ______
# FFT
X = np.fft.fft(______)
freq = np.fft.fftfreq(______, d=______) # cycle/month
# 周期→周波数
# period = 1/f なので、 cutoffより長周期を残すには |f| <= 1/cutoff
fcut = 1.0 / ______
low_mask = np.abs(freq) <= ______
high_mask = np.abs(freq) > ______
# ローパス
X_low = ______ * ______
x_low = np.fft.ifft(______).real + ______
# ハイパス
X_high = ______ * ______
x_high = np.fft.ifft(______).real
return ______, ______, ______, ______
def fft_bandpass_filter(x, dt=1.0, min_period=24.0, max_period=84.0):
"""
FFTベースのバンドパスフィルター
min_period : 残したい最短周期（月）
max_period : 残したい最長周期（月）
例：
min_period = 24
max_period = 84
→ 2〜7年周期だけ残す
"""
x = np.asarray(______, dtype=float)
if np.any(~np.isfinite(______)):
x = pd.Series(x).interpolate(limit_direction="both").values
n = len(______)
x_mean = np.mean(______)
x0 = ______ - ______
X = np.fft.fft(______)
freq = np.fft.fftfreq(______, d=______)
f_low = 1.0 / ______ # 長周期側
f_high = 1.0 / ______ # 短周期側
band_mask = (np.abs(freq) >= ______) & (np.abs(freq) <= ______)
X_band = ______ * ______
x_band = np.fft.ifft(______).real
return ______, ______, ______, ______
考えること
周期が長い成分を残したいとき、周波数の条件は大きくなるのか、小さくなるのか。
周期が長い成分を残したいとき、周波数の条件は大きくなるのか、小さくなるのか。
STEP 3
9. フィルターの実行とスペクトル
SOI と Niño3.4 に同じフィルターを適用し、パワースペクトル用に正の周波数だけを取り出す。
# =========================================================
# 3. 実行
# =========================================================
soi_low, soi_high, freq_soi, X_soi = fft_filter(
soi, dt=1.0, cutoff_period=cutoff_period_months
)
nino_low, nino_high, freq_nino, X_nino = fft_filter(
nino34, dt=1.0, cutoff_period=cutoff_period_months
)
#-- バンドパス フィルター！！！！！
band_min_period = ______
band_max_period = ______
soi_band, freq_soi_band, X_soi_band, mask_soi_band = fft_bandpass_filter(
soi,
dt=1.0,
min_period=______,
max_period=______
)
nino_band, freq_nino_band, X_nino_band, mask_nino_band = fft_bandpass_filter(
nino34,
dt=1.0,
min_period=______,
max_period=______
)
# =========================================================
# 4. パワースペクトル計算用関数
# =========================================================
def calc_power_spectrum(freq, X):
"""
正の周波数側だけ返す
"""
power = np.abs(______) ** 2
pos = ______ > 0
return ______[pos], ______[pos]
freq_soi_pos, power_soi = calc_power_spectrum(freq_soi, X_soi)
freq_nino_pos, power_nino = calc_power_spectrum(freq_nino, X_nino)
# 周期（月）に変換
period_soi = 1.0 / ______
period_nino = 1.0 / ______
# =========================================================
# 5. SOI 図
# =========================================================
fig, axes = plt.subplots(2, 1, figsize=(14, 8), constrained_layout=True)
# 時系列
axes[0].plot(______, ______, color="0.65", lw=1.0, label="Original SOI")
axes[0].plot(______, ______, color="red", lw=2.0,
label=f"Low-pass (period >= {cutoff_period_months:.0f} months)")
axes[0].plot(______, ______, color="blue", lw=1.2,
label=f"High-pass (period < {cutoff_period_months:.0f} months)")
axes[0].axhline(0, color="k", lw=0.8)
axes[0].set_title("SOI: FFT-based low-pass / high-pass filtering")
axes[0].set_ylabel("SOI")
axes[0].legend(loc="upper right")
# スペクトル
axes[1].plot(period_soi, power_soi, color="black", lw=1.5)
axes[1].axvline(cutoff_period_months, color="red", ls="--", lw=1.2,
label=f"Cutoff = {cutoff_period_months:.0f} months")
axes[1].set_xscale("log")
axes[1].set_xlabel("Period (months)")
axes[1].set_ylabel("Power")
axes[1].set_title("SOI power spectrum")
axes[1].legend()
fig.savefig(fig_soi_name, dpi=300, bbox_inches="tight", facecolor="white")
plt.show()
print(f"Saved: {fig_soi_name}")
STEP 4
10. ローパス・ハイパス図
SOI と Niño3.4 それぞれについて、元データ、ローパス、ハイパス、パワースペクトルを描く。
# =========================================================
# 6. Niño3.4 図
# =========================================================
fig, axes = plt.subplots(2, 1, figsize=(14, 8), constrained_layout=True)
# 時系列
axes[0].plot(nino_time, nino34, color="0.65", lw=1.0, label="Original Niño3.4 anomaly")
axes[0].plot(nino_time, nino_low, color="red", lw=2.0,
label=f"Low-pass (period >= {cutoff_period_months:.0f} months)")
axes[0].plot(nino_time, nino_high, color="blue", lw=1.2,
label=f"High-pass (period < {cutoff_period_months:.0f} months)")
axes[0].axhline(0, color="k", lw=0.8)
axes[0].set_title("Niño3.4: FFT-based low-pass / high-pass filtering")
axes[0].set_ylabel("Niño3.4 anomaly (°C)")
axes[0].legend(loc="upper right")
# スペクトル
axes[1].plot(period_nino, power_nino, color="black", lw=1.5)
axes[1].axvline(cutoff_period_months, color="red", ls="--", lw=1.2,
label=f"Cutoff = {cutoff_period_months:.0f} months")
axes[1].set_xscale("log")
axes[1].set_xlabel("Period (months)")
axes[1].set_ylabel("Power")
axes[1].set_title("Niño3.4 power spectrum")
axes[1].legend()
fig.savefig(fig_nino_name, dpi=300, bbox_inches="tight", facecolor="white")
plt.show()
print(f"Saved: {fig_nino_name}")
# 追加：バンドパスフィルターの結果の図
fig, axes = plt.subplots(2, 1, figsize=(14, 8), constrained_layout=True)
axes[0].plot(______, ______, color="0.65", lw=1.0, label="Original SOI")
axes[0].plot(soi_time, soi_band, color="purple", lw=2.0,
label="Band-pass: 24–84 months")
axes[0].axhline(0, color="k", lw=0.8)
axes[0].set_title("SOI: FFT-based ENSO band-pass filtering")
axes[0].set_ylabel("SOI")
axes[0].set_xlabel("Time")
axes[0].legend(loc="upper right")
axes[1].plot(nino_time, nino34, color="0.65", lw=1.0, label="Original Niño3.4 anomaly")
axes[1].plot(nino_time, nino_band, color="purple", lw=2.0,
label="Band-pass: 24–84 months")
axes[1].axhline(0, color="k", lw=0.8)
axes[1].set_title("Niño3.4: FFT-based ENSO band-pass filtering")
axes[1].set_ylabel("Niño3.4 anomaly (°C)")
axes[1].legend(loc="upper right")
plt.savefig("fft_bandpass_nino34_soi.png", dpi=300, bbox_inches="tight", facecolor="white")
plt.show()
SOI のローパス・ハイパスフィルター。
Niño3.4 のローパス・ハイパスフィルター。
STEP 5
11. バンドパス図と相関
24–84か月の ENSO 帯域だけを残し、バンドパス後の Niño3.4 と SOI の相関を計算する。
# =========================
# バンドパス後の共通期間をそろえる
# =========================
df_band = pd.DataFrame({
"time": ______,
"nino": ______,
}).merge(
pd.DataFrame({
"time": ______,
"soi": ______,
}),
on=______,
how=______
)
# NaN除去
df_band = df_band.______()
x = df_band[______].values
y = df_band[______].values
# =========================
# 相関・回帰
# =========================
r = np.corrcoef(______, ______)[0, 1]
R2 = ______**2
coef = np.polyfit(______, ______, 1)
yfit = np.polyval(______, ______)
print("=== Band-pass correlation ===")
print(f"r = {r:.2f}")
print(f"R^2 = {R2:.2f}")
print(f"SOI = {coef[0]:.2f} × Niño3.4 + {coef[1]:.2f}")
# =========================
# 描画
# =========================
plt.figure(figsize=(6,6))
plt.scatter(x, y, s=25, alpha=0.6)
plt.plot(x, yfit, color="red", lw=2)
plt.xlabel("Niño3.4 anomaly (band-pass)")
plt.ylabel("SOI (band-pass)")
plt.title("Band-pass filtered Niño3.4 vs SOI")
plt.text(
0.05, 0.95,
f"r = {r:.2f}\nR² = {R2:.2f}",
transform=plt.gca().transAxes,
verticalalignment="top"
)
plt.grid()
plt.tight_layout()
plt.savefig("scatter_bandpass_nino34_soi.png", dpi=300)
plt.show()
lags = np.arange(-12, 13)
r_lag = []
for lag in lags:
if lag > 0:
rtmp = np.corrcoef(x[lag:], y[:-lag])[0,1]
elif lag < 0:
rtmp = np.corrcoef(x[:lag], y[-lag:])[0,1]
else:
rtmp = np.corrcoef(x, y)[0,1]
r_lag.append(rtmp)
plt.plot(lags, r_lag)
plt.axhline(0, color='k')
plt.xlabel("Lag (months)")
plt.ylabel("Correlation")
plt.title("Lag correlation (band-pass)")
plt.show()
#-----有意性の確認
SOI と Niño3.4 の 24–84か月バンドパス成分。
バンドパス後の散布図。r は非常に大きいが、ここで終わらず自由度を見る。
STEP 6
12. 有効自由度と p 値
フィルター後の時系列はなめらかになるため、隣り合うデータが独立でなくなる。ここでは lag-1 自己相関から有効自由度を見積もる。
# =========================================================
# 8. オリジナルとバンドパス後の相関・有意性の比較
# =========================================================
def effective_dof(x, y):
"""
lag-1自己相関を使った有効自由度の推定
Bretherton et al. 型の簡易補正
"""
x = np.asarray(______, dtype=float)
y = np.asarray(y, dtype=float)
r1 = np.corrcoef(______, ______)[0, 1]
r2 = np.corrcoef(______, ______)[0, 1]
N = len(x)
Neff = ______ * (1 - ______ * ______) / (1 + ______ * ______)
# 念のため下限を設定
Neff = max(3, Neff)
return Neff, r1, r2
def corr_significance(r, N):
"""
相関係数 r の t 検定
N は通常のサンプル数でも、有効自由度でもよい
"""
t = ______ * np.sqrt((______ - 2) / (1 - ______**2))
p = 2 * (1 - stats.t.cdf(np.abs(______), df=______ - 2))
return t, p
def corr_regression_stats(x, y, use_neff=False):
"""
相関・回帰・有意性をまとめて計算
"""
x = np.asarray(______, dtype=float)
y = np.asarray(y, dtype=float)
mask = np.isfinite(x) & np.isfinite(y)
x = x[mask]
y = y[mask]
N = len(x)
r = np.corrcoef(______, ______)[0, 1]
R2 = ______**2
slope, intercept = np.polyfit(______, ______, 1)
if use_neff:
Neff, r1, r2 = effective_dof(x, y)
t, p = corr_significance(r, Neff)
else:
Neff = N
r1 = np.nan
r2 = np.nan
t, p = corr_significance(r, N)
return {
"x": x,
"y": y,
"N": N,
"Neff": Neff,
"r": r,
"R2": R2,
"slope": slope,
"intercept": intercept,
"t": t,
"p": p,
"r1_x": r1,
"r1_y": r2,
}
普通の相関検定で N=519 をそのまま使うと、有意性を過大評価する可能性がある。バンドパス後は Neff が一桁程度まで落ちることがある。
STEP 7
13. 結果表示と散布図比較
最後に、オリジナルとバンドパス後で、N、Neff、r、R²、p を並べて比較する。
# =========================================================
# 8-1. オリジナル：timeで共通期間をそろえる
# =========================================================
df_org = pd.DataFrame({
"time": ______,
"nino": nino34
}).merge(
pd.DataFrame({
"time": ______,
"soi": soi
}),
on=______,
how=______
).dropna()
org_stats = corr_regression_stats(
df_org["nino"].values,
df_org["soi"].values,
use_neff=______
)
# =========================================================
# 8-2. バンドパス後：timeで共通期間をそろえる
# =========================================================
df_bp = pd.DataFrame({
"time": ______,
"nino": nino_band
}).merge(
pd.DataFrame({
"time": ______,
"soi": soi_band
}),
on=______,
how=______
).dropna()
bp_stats = corr_regression_stats(
df_bp["nino"].values,
df_bp["soi"].values,
use_neff=______
)
# =========================================================
# 8-3. 結果表示
# =========================================================
print("===== ORIGINAL =====")
print(f"N = {org_stats['N']}")
print(f"Neff = {org_stats['Neff']:.1f}")
print(f"lag-1 autocorr Niño3.4 = {org_stats['r1_x']:.3f}")
print(f"lag-1 autocorr SOI = {org_stats['r1_y']:.3f}")
print(f"r = {org_stats['r']:.3f}")
print(f"R^2 = {org_stats['R2']:.3f}")
print(f"SOI = {org_stats['slope']:.2f} × Niño3.4 + {org_stats['intercept']:.2f}")
print(f"t = {org_stats['t']:.2f}")
print(f"p = {org_stats['p']:.3e}")
print("\n===== BAND-PASS 24–84 months =====")
print(f"N raw = {bp_stats['N']}")
print(f"Neff = {bp_stats['Neff']:.1f}")
print(f"lag-1 autocorr Niño3.4 = {bp_stats['r1_x']:.3f}")
print(f"lag-1 autocorr SOI = {bp_stats['r1_y']:.3f}")
print(f"r = {bp_stats['r']:.3f}")
print(f"R^2 = {bp_stats['R2']:.3f}")
print(f"SOI = {bp_stats['slope']:.2f} × Niño3.4 + {bp_stats['intercept']:.2f}")
print(f"t = {bp_stats['t']:.2f}")
print(f"p = {bp_stats['p']:.3e}")
# =========================================================
# 8-4. オリジナル vs バンドパス 散布図比較
# =========================================================
fig, axes = plt.subplots(1, 2, figsize=(12, 5), constrained_layout=True)
# ---- Original ----
x = org_stats["x"]
y = org_stats["y"]
xline = np.linspace(np.min(______), np.max(______), 200)
yline = org_stats[______] * xline + org_stats[______]
axes[0].scatter(______, ______, s=25, alpha=0.55)
axes[0].plot(xline, yline, color="red", lw=2)
axes[0].axhline(0, color="k", lw=0.8)
axes[0].axvline(0, color="k", lw=0.8)
axes[0].grid(True, alpha=0.4)
axes[0].set_title("Original")
axes[0].set_xlabel("Niño3.4 SST anomaly")
axes[0].set_ylabel("SOI")
axes[0].text(
0.05, 0.95,
f"r = {org_stats['r']:.2f}\n"
f"R² = {org_stats['R2']:.2f}\n"
f"p = {org_stats['p']:.1e}",
transform=axes[0].transAxes,
va="top",
bbox=dict(facecolor="white", edgecolor="none", alpha=0.8)
)
# ---- Band-pass ----
x = bp_stats["x"]
y = bp_stats["y"]
xline = np.linspace(np.min(______), np.max(______), 200)
yline = bp_stats["slope"] * xline + bp_stats["intercept"]
axes[1].scatter(x, y, s=25, alpha=0.55)
axes[1].plot(xline, yline, color="red", lw=2)
axes[1].axhline(0, color="k", lw=0.8)
axes[1].axvline(0, color="k", lw=0.8)
axes[1].grid(True, alpha=0.4)
axes[1].set_title("Band-pass filtered 24–84 months")
axes[1].set_xlabel("Niño3.4 anomaly (band-pass)")
axes[1].set_ylabel("SOI (band-pass)")
axes[1].text(
0.05, 0.95,
f"r = {bp_stats['r']:.2f}\n"
f"R² = {bp_stats['R2']:.2f}\n"
f"Neff = {bp_stats['Neff']:.1f}\n"
f"p = {bp_stats['p']:.1e}",
transform=axes[1].transAxes,
va="top",
bbox=dict(facecolor="white", edgecolor="none", alpha=0.8)
)
plt.savefig("scatter_original_vs_bandpass_nino34_soi.png",
dpi=300, bbox_inches="tight", facecolor="white")
plt.show()
print("Saved: scatter_original_vs_bandpass_nino34_soi.png")
オリジナルとバンドパス後の比較。バンドパス後は r が大きい一方で、Neff が大きく減る。
14. まとめ
- 相関係数は、2 つの時系列が一緒に変動する程度を表す。
- FFT フィルターにより、月々の変動、低周波変動、ENSO 帯域の変動を分けて考えられる。
- フィルター後は相関が強く見えることがあるが、自由度は小さくなる。
- したがって、r だけでなく、Neff と p を一緒に見る必要がある。
15. 解答の表示
パスワードを入力すると、穴埋めなしの完全スクリプトを表示する。
16. 自由度の近似式について
このページでは、相関係数の検定で、見かけのデータ数 N をそのまま使わず、自己相関を考慮した有効自由度を使った。
Neff = N × (1 − r1x r1y) / (1 + r1x r1y)
この式は突然出てきたように見えるが、より一般的な有効サンプルサイズの式に、いくつかの近似を入れたものである。
1. まず、なぜ有効自由度が必要なのか
通常の相関係数の検定では、データが互いに独立であることを仮定する。しかし、SOI や Niño3.4 のような気候時系列では、隣り合う月の値が似ている。さらにローパスやバンドパスをかけると、時系列はよりなめらかになり、隣り合う値はますます独立ではなくなる。
したがって、見かけ上は N 個の月データがあっても、統計的に独立な情報は N 個より少ない。
この「独立なデータ数に換算した値」を Neff と考える。
この「独立なデータ数に換算した値」を Neff と考える。
2. 一般的な有効サンプルサイズの式
2つの時系列 x と y の相関を検定するとき、自己相関を考慮した有効サンプルサイズは、一般には次のように書ける。
Neff ≈ N / [ 1 + 2 Σk=1N−1 (1 − k/N) ρx(k) ρy(k) ]
| 記号 | 意味 |
|---|---|
| ρx(k) | x の lag k における自己相関 |
| ρy(k) | y の lag k における自己相関 |
| (1 − k/N) | 有限長の時系列で、lag が大きくなるほど使えるペア数が減ることを表す補正 |
この式は、「1か月後、2か月後、3か月後……」というすべての lag の自己相関の影響を足し合わせて、実質的な独立データ数を小さくする式である。
3. Bretherton et al. (1999) の書き方
Bretherton et al. (1999) では、上と同じ考え方を、正の lag と負の lag をまとめて次のような形で書いている。
T*XY = T / Στ=−(T−1)T−1 (1 − |τ|/T) ρX(τ) ρY(τ)
この形では、τ = 0 の項が 1 になるため、分母の中に 1 + 2Σ が含まれていると考えればよい。見た目は違うが、正の lag だけで書いた式と同じ考え方である。
4. AR(1) 近似を入れる
このページのコードでは、すべての lag の自己相関を個別に推定するのではなく、lag-1 自己相関だけで代表させる。これは、時系列が AR(1)、つまり赤色雑音的にふるまうと仮定していることに対応する。
ρx(k) ≈ r1xk
ρy(k) ≈ r1yk
ρy(k) ≈ r1yk
ここで、
q = r1x r1y
とおくと、
ρx(k) ρy(k) ≈ qk
となる。
5. 十分長い時系列として近似する
さらに、時系列が十分長く、自己相関がデータ長よりかなり短い lag で減衰すると考える。このとき、有限長補正
1 − k/N
をおおよそ 1 とみなし、和の上限も無限大まで伸ばす。
Neff ≈ N / [ 1 + 2 Σk=1∞ qk ]
等比級数の公式より、
Σk=1∞ qk = q / (1 − q)
なので、
Neff ≈ N / [ 1 + 2q/(1−q) ]
分母を整理すると、
1 + 2q/(1−q)
= (1−q)/(1−q) + 2q/(1−q)
= (1+q)/(1−q)
= (1−q)/(1−q) + 2q/(1−q)
= (1+q)/(1−q)
したがって、
Neff ≈ N × (1−q)/(1+q)
最後に q = r1x r1y を戻すと、
Neff ≈ N × (1 − r1x r1y) / (1 + r1x r1y)
これが、このページの Python コードで使っている式である。
6. なぜ「近似式」なのか
この式は便利だが、厳密に常に正しい式ではない。主な近似は次の3つである。
| 近似 | 内容 | 注意点 |
|---|---|---|
| AR(1) 近似 | lag-1 自己相関だけで全 lag の自己相関を代表させる。 | ENSO のように周期性が強い場合、AR(1) だけでは表しきれないことがある。 |
| 十分長い時系列の近似 | 1 − k/N を 1 に近いとみなす。 | 短い時系列では、この近似が粗くなる。 |
| 定常性の仮定 | 平均や分散、自己相関構造が時間とともに大きく変わらないと仮定する。 | 強いトレンドや非定常性がある場合は注意が必要。 |
7. このページでの使い方
このページでは、教育用・実用的な近似として、次のように計算している。
r1x = np.corrcoef(x[:-1], x[1:])[0, 1]
r1y = np.corrcoef(y[:-1], y[1:])[0, 1]
Neff = N * (1 - r1x * r1y) / (1 + r1x * r1y)その後、相関係数の t 検定では、N の代わりに Neff を使う。
t = r √{ (Neff − 2) / (1 − r²) }
df = Neff − 2
df = Neff − 2
重要
Neff は「有効サンプルサイズ」に近い量であり、t 分布に入れる自由度そのものは Neff − 2 である。
たとえば Neff = 7 の場合、相関検定で使う自由度はおよそ 5 である。
Neff は「有効サンプルサイズ」に近い量であり、t 分布に入れる自由度そのものは Neff − 2 である。
たとえば Neff = 7 の場合、相関検定で使う自由度はおよそ 5 である。
8. 07での解釈
バンドパス後の Niño3.4 と SOI では、相関係数の絶対値が大きくなることがある。しかし、バンドパスによって時系列はなめらかになり、自己相関も強くなるため、Neff は大きく減る。
結論
フィルター後の相関は、r だけを見ると強く見える。
しかし、独立な情報量は少なくなる。
したがって、r、Neff、p をセットで見る必要がある。
フィルター後の相関は、r だけを見ると強く見える。
しかし、独立な情報量は少なくなる。
したがって、r、Neff、p をセットで見る必要がある。
参考文献
この節の考え方は、自己相関時系列に対する effective sample size の議論に基づく。Bretherton et al. (1999) では、時系列の自己相関を考慮した有効サンプルサイズの一般式と、赤色雑音・Markov 過程の場合の簡略式が説明されている。