07 フィルタリング・自由度・P値

穴埋め解答

24.0
24.0
84.0
fcut
fcut
f_low
f_high

完全スクリプト

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.0
band_min_period = 24.0
band_max_period = 84.0

# =========================================================
# 1. データ読み込み
# =========================================================
soi_df = pd.read_csv(soi_file, sep=r"\s+")
soi_time = pd.to_datetime(dict(year=soi_df["YEAR"], month=soi_df["MONTH"], day=15))
soi = soi_df["SOI"].values.astype(float)

colnames = [
    "YR", "MON",
    "NINO12", "ANOM12",
    "NINO3",  "ANOM3",
    "NINO4",  "ANOM4",
    "NINO34", "ANOM34"
]

nino_df = pd.read_csv(
    nino_file,
    sep=r"\s+",
    skiprows=1,
    names=colnames
)

nino_time = pd.to_datetime(dict(year=nino_df["YR"], month=nino_df["MON"], day=15))
nino34 = nino_df["ANOM34"].values.astype(float)

# =========================================================
# 2. フィルター関数
# =========================================================
def fft_filter(x, dt=1.0, cutoff_period=24.0):
    x = np.asarray(x, dtype=float)

    if np.any(~np.isfinite(x)):
        x = pd.Series(x).interpolate(limit_direction="both").values

    n = len(x)
    x_mean = np.mean(x)
    x0 = x - x_mean

    X = np.fft.fft(x0)
    freq = np.fft.fftfreq(n, d=dt)

    fcut = 1.0 / cutoff_period

    low_mask = np.abs(freq) <= fcut
    high_mask = np.abs(freq) > fcut

    x_low = np.fft.ifft(X * low_mask).real + x_mean
    x_high = np.fft.ifft(X * high_mask).real

    return x_low, x_high, freq, X


def fft_bandpass_filter(x, dt=1.0, min_period=24.0, max_period=84.0):
    x = np.asarray(x, dtype=float)

    if np.any(~np.isfinite(x)):
        x = pd.Series(x).interpolate(limit_direction="both").values

    n = len(x)
    x_mean = np.mean(x)
    x0 = x - x_mean

    X = np.fft.fft(x0)
    freq = np.fft.fftfreq(n, d=dt)

    f_low = 1.0 / max_period
    f_high = 1.0 / min_period

    band_mask = (np.abs(freq) >= f_low) & (np.abs(freq) <= f_high)
    x_band = np.fft.ifft(X * band_mask).real

    return x_band, freq, X, band_mask


def calc_power_spectrum(freq, X):
    power = np.abs(X) ** 2
    pos = freq > 0
    return freq[pos], power[pos]

# =========================================================
# 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)

soi_band, freq_soi_band, X_soi_band, mask_soi_band = fft_bandpass_filter(
    soi, dt=1.0, min_period=band_min_period, max_period=band_max_period
)
nino_band, freq_nino_band, X_nino_band, mask_nino_band = fft_bandpass_filter(
    nino34, dt=1.0, min_period=band_min_period, max_period=band_max_period
)

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 / freq_soi_pos
period_nino = 1.0 / freq_nino_pos

# =========================================================
# 4. ローパス・ハイパス図
# =========================================================
fig, axes = plt.subplots(2, 1, figsize=(14, 8), constrained_layout=True)
axes[0].plot(soi_time, soi, color="0.65", lw=1.0, label="Original SOI")
axes[0].plot(soi_time, soi_low, color="red", lw=2.0, label=f"Low-pass (period >= {cutoff_period_months:.0f} months)")
axes[0].plot(soi_time, soi_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("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("fft_filter_soi.png", dpi=300, bbox_inches="tight", facecolor="white")
plt.show()

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("fft_filter_nino34.png", dpi=300, bbox_inches="tight", facecolor="white")
plt.show()

# =========================================================
# 5. バンドパス図
# =========================================================
fig, axes = plt.subplots(2, 1, figsize=(14, 8), constrained_layout=True)
axes[0].plot(soi_time, soi, 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].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].set_xlabel("Time")
axes[1].legend(loc="upper right")
fig.savefig("fft_bandpass_nino34_soi.png", dpi=300, bbox_inches="tight", facecolor="white")
plt.show()

# =========================================================
# 6. 相関・有意性
# =========================================================
def effective_dof(x, y):
    r1 = np.corrcoef(x[:-1], x[1:])[0, 1]
    r2 = np.corrcoef(y[:-1], y[1:])[0, 1]
    N = len(x)
    Neff = N * (1 - r1 * r2) / (1 + r1 * r2)
    Neff = max(3, Neff)
    return Neff, r1, r2


def corr_significance(r, N):
    t = r * np.sqrt((N - 2) / (1 - r**2))
    p = 2 * (1 - stats.t.cdf(np.abs(t), df=N - 2))
    return t, p


def corr_regression_stats(x, y, use_neff=True):
    x = np.asarray(x, 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(x, y)[0, 1]
    R2 = r**2
    slope, intercept = np.polyfit(x, y, 1)
    if use_neff:
        Neff, r1, r2 = effective_dof(x, y)
        t, p = corr_significance(r, Neff)
    else:
        Neff, r1, r2 = N, np.nan, np.nan
        t, p = corr_significance(r, N)
    return dict(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)

# オリジナル：共通期間でそろえる
df_org = pd.DataFrame({"time": nino_time, "nino": nino34}).merge(
    pd.DataFrame({"time": soi_time, "soi": soi}),
    on="time", how="inner"
).dropna()

org_stats = corr_regression_stats(df_org["nino"].values, df_org["soi"].values, use_neff=True)

# バンドパス後：共通期間でそろえる
df_bp = pd.DataFrame({"time": nino_time, "nino": nino_band}).merge(
    pd.DataFrame({"time": soi_time, "soi": soi_band}),
    on="time", how="inner"
).dropna()

bp_stats = corr_regression_stats(df_bp["nino"].values, df_bp["soi"].values, use_neff=True)

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}")

# 散布図比較
fig, axes = plt.subplots(1, 2, figsize=(12, 5), constrained_layout=True)

for ax, stats_dict, title, xlabel, ylabel in [
    (axes[0], org_stats, "Original", "Niño3.4 SST anomaly", "SOI"),
    (axes[1], bp_stats, "Band-pass filtered 24–84 months", "Niño3.4 anomaly (band-pass)", "SOI (band-pass)")
]:
    x = stats_dict["x"]
    y = stats_dict["y"]
    xline = np.linspace(np.min(x), np.max(x), 200)
    yline = stats_dict["slope"] * xline + stats_dict["intercept"]
    ax.scatter(x, y, s=25, alpha=0.55)
    ax.plot(xline, yline, color="red", lw=2)
    ax.axhline(0, color="k", lw=0.8)
    ax.axvline(0, color="k", lw=0.8)
    ax.grid(True, alpha=0.4)
    ax.set_title(title)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    ax.text(
        0.05, 0.95,
        f"r = {stats_dict['r']:.2f}\nR² = {stats_dict['R2']:.2f}\nNeff = {stats_dict['Neff']:.1f}\np = {stats_dict['p']:.1e}",
        transform=ax.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()