01.5 Python計算トレーニング

1. なぜこの練習が必要か

01 では Python 環境の構築を行い、02 ではすぐに SOI の計算に進む。しかし、SOI の計算では実は多くの操作をまとめて行っている。

• データを読む
• 平均を求める
• 平均からのずれを作る
• 標準偏差を求める
• 標準化する
• 2系列の差を取る
• 移動平均をかける

そこで、まずは小さなサンプルデータで、これらの計算を1つずつ確認する。

Pythonでよく出る書き方

2. 練習用データ

ここでは 12 か月分の簡単な配列を、仮の時系列データとして使う。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

x = np.array([10, 11, 13, 12, 14, 15, 16, 15, 14, 13, 12, 11], dtype=float)
y = np.array([ 8,  9, 10, 11, 12, 14, 15, 14, 13, 12, 10,  9], dtype=float)
month = np.arange(1, 13)

print(x)
print(y)

3. 配列と平均

まずは平均を計算する。

x_mean = np.mean(x)
y_mean = np.mean(y)

print("x mean =", x_mean)
print("y mean =", y_mean)

確認

• 平均とは何か
• 平均を取ると、時系列全体の代表値が1つ得られることを確認する

4. 偏差（anomaly）

平均を引くと、各値が「平均からどれだけずれているか」がわかる。これを anomaly と考える。

x_anom = x - x_mean
y_anom = y - y_mean

print(x_anom)
print(y_anom)
print("mean of x_anom =", np.mean(x_anom))

確認

• anomaly の平均は 0 に近くなる
• 正の値は平均より大きい、負の値は平均より小さいことを表す

5. 標準偏差と標準化

データのばらつきを表す量が標準偏差である。標準偏差で割ると、スケールをそろえた「標準化データ」になる。

x_std = np.sqrt(np.sum(x_anom**2) / x_anom.size)
y_std = np.sqrt(np.sum(y_anom**2) / y_anom.size)

x_norm = x_anom / x_std
y_norm = y_anom / y_std

print("x std =", x_std)
print("y std =", y_std)
print(x_norm)
print(y_norm)

確認

• 標準偏差が大きいほど、データのばらつきが大きい
• 標準化すると、単位の違うデータでも比較しやすくなる

6. 2系列の差

SOIでは、標準化した Tahiti と Darwin の差を用いる。ここでも同じことを試す。

diff = y_norm - x_norm
print(diff)
print("mean of diff =", np.mean(diff))

この差のばらつきを求めると、SOI計算に出てくる MSD と同じ考え方になる。

msd = np.sqrt(np.sum(diff**2) / diff.size)
print("MSD =", msd)

確認

• 2系列の差が大きい月はどこか
• MSD は「差のばらつき」を表していることを確認する

7. 移動平均

次に、時系列のノイズを少しなめらかにするため、移動平均を行う。

s = pd.Series(x)
x_ma3 = s.rolling(window=3, center=True).mean()

print(x_ma3)

図で確認する。

plt.figure(figsize=(8,4))
plt.plot(month, x, marker="o", label="original")
plt.plot(month, x_ma3, marker="o", label="3-point moving average")
plt.xlabel("Month")
plt.ylabel("Value")
plt.title("Moving Average")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

確認

• 移動平均をかけるとギザギザがなめらかになる
• 前後の端では値が欠けることがある

8. 自己相関の考え方

自己相関とは、「同じ時系列を少しずらして比べたとき、どれくらい似ているか」を見る量である。

def autocorrelation_simple(a, lag):
    x1 = a[:-lag]
    x2 = a[lag:]
    return np.corrcoef(x1, x2)[0,1]

for lag in [1, 2, 3]:
    print(lag, autocorrelation_simple(x, lag))

ラグを増やしながら自己相関を描く。

lags = np.arange(1, 6)
acf = [autocorrelation_simple(x, lag) for lag in lags]

plt.figure(figsize=(7,4))
plt.bar(lags, acf)
plt.axhline(0, color="k")
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.title("Autocorrelation (sample data)")
plt.tight_layout()
plt.show()

9. SOI計算との対応

付録：演習の解答例

STEP1 演習

z = np.array([5, 6, 7, 8, 9], dtype=float)
print("z mean =", np.mean(z))

STEP3 演習（標準化の確認）

print("mean(x_norm) =", np.mean(x_norm))
print("std(x_norm)  =", np.sqrt(np.sum((x_norm - np.mean(x_norm))**2) / x_norm.size))

STEP5 演習（window変更）

x_ma5 = pd.Series(x).rolling(window=5, center=True).mean()

plt.figure(figsize=(8,4))
plt.plot(month, x, marker="o", label="original")
plt.plot(month, x_ma5, marker="o", label="5-point moving average")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

STEP6 自己相関

lags = np.arange(1, 6)
acf = [autocorrelation_simple(x, lag) for lag in lags]
print(list(zip(lags, acf)))