In [1]:
import xml.etree.ElementTree as ET
import pandas as pd
import numpy as np
from pathlib import Path
from scipy.optimize import curve_fit
from statsmodels.tsa.ar_model import AutoReg
import matplotlib.pyplot as plt
from matplotlib import rcParams
# 日本語フォント対応（Windows環境用：游ゴシック or MS Gothic）
rcParams['font.family'] = 'Yu Gothic' # or 'MS Gothic'
p4p 形式 XML データの読み込み
In [2]:
def load_p4p_series(path, core_name, first_year):
"""
.p4p の <DATA> を読み込み、Year を index にした Series を返す。
MissingRing==1 行は NaN。
Parameters
----------
path : str # .p4p ファイルパス
core_name : str # コア名
first_year : int# 先頭の暦年（手動入力）
Returns
-------
pd.Series (index=Year, value=RingWidth or NaN)
"""
xml = ET.parse(path).getroot()
# CDATA 配列を抽出
data_text = xml.find(".//RECORD/DATA").text
rows = [r.strip().split('\t') for r in data_text.strip().splitlines()]
df = pd.DataFrame(rows, columns=["Year","Flag1","Flag2","Missing","RingWidth"]).astype(float)
# Past の Year 列は “1,2,3…” の場合がある → first_year で上書き
df["Year"] = first_year + np.arange(len(df))
#df.loc[df["Missing"] == 1, "RingWidth"] = 0.001 # 欠損データは極小値で埋める
return df.set_index("Year")["RingWidth"].rename(core_name)
樹木の年輪幅 (Ring Width, RW) には 「樹齢に伴う長期減衰トレンド」 と 「気候など外部要因による年々の振動成分」 が重なって入っています。 年輪年代法（クロスデーティング）が狙うのは 後者＝気候シグナルの同期 であって、前者はむしろノイズになります。
そのため実務では、樹齢トレンドを取り除き、年々変動だけを取り出した「リング幅指標 (Ring-Width Index, RWI)」 を作った上で相関や t 値（Baillie & Pilcher tBP）・GLK (Gleichläufigkeit) を計算します。 この前処理を入れないと、たいてい相関は 0.3〜0.4 止まりで「合っていないように見える」ことが多いです。
In [3]:
# COFECHA 既定デトレンド & AR(1) 残差
from scipy.optimize import curve_fit
from statsmodels.tsa.ar_model import AutoReg
def _neg_exp(x,a,b,k): return a*np.exp(b*x)+k # y=a·e^{bx}+k
def detrend_cofecha(series):
y, x = series.values, np.arange(len(series))
try: # 1) 負の指数
p0 = (y[0], -0.01, y[-1])
popt,_ = curve_fit(_neg_exp, x, y, p0=p0, maxfev=10000)
trend = _neg_exp(x,*popt)
if (trend<=0).any(): raise RuntimeError
except Exception: # 2) 直線 or 3) 全平均
m,c = np.polyfit(x,y,1)
trend = m*x+c if (m!=0) else np.full_like(y, y.mean())
return pd.Series(y/trend, index=series.index)
def prewhiten_ar1(idx_series, lag=1):
"""
AR(1) で自己相関を除去し残差を返す。
先頭 lag 行の NaN を直後の値で埋める (bfill)。
"""
s = idx_series.dropna()
if len(s) < lag + 2:
# データ不足なら全 NaN
return pd.Series(np.nan, index=idx_series.index)
# モデル用に連番インデックス
tmp = s.copy()
tmp.index = pd.RangeIndex(len(tmp))
model = AutoReg(tmp, lags=lag, old_names=False).fit()
resid = model.resid # 長さ len(tmp)-lag
# 元の暦年インデックスに戻し
resid.index = s.index[lag:]
resid = resid.reindex(idx_series.index)
# 先頭 NaN を直後値で埋める
resid = resid.bfill()
return resid
In [4]:
def compute_yearly_stats(master_resid, core_resid, start_year):
n = len(core_resid)
rows = []
for lag in range(len(master_resid) - n + 1):
# 位置合わせで切り出し
seg_m = master_resid.iloc[lag:lag+n].reset_index(drop=True)
seg_r = core_resid.reset_index(drop=True)
# 相関と tBP
if seg_m.std() == 0 or seg_r.std() == 0:
r, t = 0.0, 0.0
else:
r = seg_m.corr(seg_r)
t = np.nan if abs(r) >= 1.0 else r * np.sqrt(n-2) / np.sqrt(1 - r*r)
# GLK
glk = (np.sign(np.diff(seg_m.values)) == np.sign(np.diff(seg_r.values))).mean()
end_year = start_year + lag + n - 1
rows.append({"end_year": end_year, "r": r, "tBP": t, "GLK": glk})
return pd.DataFrame(rows)
参照データの読み込み¶
In [5]:
# 暦年標準Eパターン（紀元前37年～838年）
start_year_known = -37
filepath = 'data/36A.p4p'
ser_known = load_p4p_series(filepath, "known", start_year_known)
print(ser_known[:5])
print(ser_known[-5:])
Year -37 47.0 -36 46.0 -35 48.0 -34 49.0 -33 42.0 Name: known, dtype: float64 Year 833 83.0 834 98.0 835 80.0 836 45.0 837 91.0 Name: known, dtype: float64
In [6]:
# デトレンド + 残差
resid_known = prewhiten_ar1(detrend_cofecha(ser_known))
比較データの読み込み¶
In [7]:
# 纏向石塚古墳板材
filepath = 'data/035￥54/54.p4p'
ser_unknown = load_p4p_series(filepath, "unknown", 0)
print(ser_unknown[:5])
print(ser_unknown[-5:])
Year 0 142.0 1 105.0 2 121.0 3 99.0 4 142.0 Name: unknown, dtype: float64 Year 241 45.0 242 62.0 243 30.0 244 68.0 245 85.0 Name: unknown, dtype: float64
In [8]:
# デトレンド + 残差
resid_unknown = prewhiten_ar1(detrend_cofecha(ser_unknown))
In [9]:
# known を“暫定マスター”とみなして unknown を検定
stats_df = compute_yearly_stats(resid_known, resid_unknown[-200:], start_year_known) # resid_unknown[-200:] として外周200年分を比較
display(stats_df.head(3))
display(stats_df.tail(3))
| end_year | r | tBP | GLK | |
|---|---|---|---|---|
| 0 | 162 | -0.078677 | -1.110532 | 0.457286 |
| 1 | 163 | -0.023420 | -0.329637 | 0.572864 |
| 2 | 164 | 0.020628 | 0.290326 | 0.457286 |
| end_year | r | tBP | GLK | |
|---|---|---|---|---|
| 673 | 835 | -0.049806 | -0.701710 | 0.467337 |
| 674 | 836 | -0.086190 | -1.217330 | 0.477387 |
| 675 | 837 | 0.118432 | 1.678292 | 0.527638 |
In [10]:
from scipy.signal import find_peaks
# --- 3) 局所ピーク検出 ---
peaks_r, _ = find_peaks(stats_df["r"], prominence=0.1)
peaks_t, _ = find_peaks(np.nan_to_num(stats_df["tBP"], nan=-np.inf), prominence=0.5)
peaks_g, _ = find_peaks(stats_df["GLK"], prominence=0.05)
years_r = stats_df["end_year"].iloc[peaks_r].tolist()
years_t = stats_df["end_year"].iloc[peaks_t].tolist()
years_g = stats_df["end_year"].iloc[peaks_g].tolist()
print("r のピーク年:", years_r)
print("tBP のピーク年:", years_t)
print("GLK のピーク年:", years_g)
r のピーク年: [169, 174, 180, 183, 194, 196, 204, 210, 212, 220, 223, 226, 234, 237, 240, 245, 250, 252, 257, 259, 266, 269, 272, 276, 283, 295, 297, 300, 304, 309, 321, 323, 329, 332, 337, 340, 346, 352, 355, 359, 371, 378, 383, 388, 390, 395, 404, 418, 423, 430, 438, 446, 453, 455, 464, 466, 469, 472, 475, 477, 489, 492, 497, 502, 505, 509, 521, 524, 526, 531, 534, 538, 540, 545, 549, 553, 556, 559, 561, 570, 572, 577, 584, 592, 595, 598, 600, 605, 609, 611, 616, 619, 625, 628, 630, 638, 641, 647, 649, 657, 661, 664, 670, 672, 680, 687, 689, 694, 698, 701, 704, 708, 710, 714, 716, 724, 730, 733, 738, 741, 745, 751, 754, 765, 768, 771, 779, 781, 783, 785, 788, 799, 802, 804, 812, 815, 818, 827, 829, 831, 834] tBP のピーク年: [164, 167, 169, 171, 174, 176, 180, 183, 187, 191, 194, 196, 201, 204, 207, 210, 212, 216, 218, 220, 223, 226, 229, 231, 234, 237, 240, 243, 245, 250, 252, 257, 259, 264, 266, 269, 272, 276, 281, 283, 285, 288, 290, 292, 295, 297, 300, 302, 304, 306, 309, 311, 317, 319, 321, 323, 326, 329, 332, 337, 340, 342, 346, 348, 352, 355, 357, 359, 361, 367, 371, 378, 381, 383, 386, 388, 390, 392, 395, 397, 404, 406, 411, 413, 418, 421, 423, 427, 430, 432, 436, 438, 440, 442, 446, 449, 453, 455, 458, 460, 464, 466, 469, 472, 475, 477, 484, 489, 492, 494, 497, 502, 505, 507, 509, 513, 518, 521, 524, 526, 529, 531, 534, 536, 538, 540, 542, 545, 549, 553, 556, 559, 561, 564, 570, 572, 575, 577, 580, 584, 586, 589, 592, 595, 598, 600, 603, 605, 609, 611, 614, 616, 619, 622, 625, 628, 630, 634, 638, 641, 644, 647, 649, 653, 657, 661, 664, 666, 668, 670, 672, 678, 680, 683, 685, 687, 689, 691, 694, 698, 701, 704, 708, 710, 714, 716, 718, 722, 724, 727, 730, 733, 735, 738, 741, 743, 745, 747, 751, 754, 758, 760, 765, 768, 771, 773, 775, 779, 781, 783, 785, 788, 791, 797, 799, 802, 804, 807, 809, 812, 815, 818, 821, 827, 829, 831, 834] GLK のピーク年: [163, 166, 174, 176, 180, 185, 187, 190, 193, 200, 203, 205, 210, 212, 214, 220, 223, 226, 228, 234, 237, 240, 243, 245, 252, 257, 259, 262, 264, 269, 272, 274, 276, 279, 287, 292, 295, 299, 301, 304, 306, 309, 312, 314, 317, 321, 324, 326, 329, 332, 337, 343, 346, 348, 350, 353, 355, 357, 359, 361, 364, 366, 368, 373, 375, 381, 383, 386, 388, 391, 393, 395, 404, 406, 411, 413, 418, 421, 423, 425, 428, 430, 433, 436, 438, 442, 450, 453, 455, 464, 466, 469, 472, 475, 477, 485, 487, 489, 492, 496, 499, 502, 504, 511, 519, 521, 523, 526, 529, 531, 534, 536, 538, 540, 543, 549, 554, 556, 559, 561, 563, 565, 580, 583, 588, 590, 592, 596, 598, 600, 603, 605, 607, 609, 611, 619, 622, 625, 628, 630, 634, 642, 645, 647, 649, 652, 655, 657, 659, 661, 664, 670, 672, 678, 680, 685, 687, 689, 694, 698, 701, 704, 707, 712, 714, 716, 718, 721, 724, 726, 728, 730, 732, 735, 738, 741, 743, 745, 751, 754, 758, 760, 763, 765, 768, 771, 773, 775, 779, 785, 787, 794, 796, 799, 802, 804, 807, 809, 814, 817, 827, 829, 831, 834]
In [11]:
fig, axs = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
# (a) r
axs[0].plot(stats_df["end_year"], stats_df["r"], label="r")
axs[0].scatter(years_r, stats_df["r"].iloc[peaks_r], color="red")
for yr in years_r:
axs[0].text(yr, stats_df.loc[stats_df.end_year==yr, "r"].values[0],
str(int(yr)), ha="center", va="bottom", fontsize=8)
axs[0].axhline(0.5, color="gray", ls="--")
axs[0].set_ylabel("r")
axs[0].set_title("Correlation (r) Peaks")
axs[0].grid(True)
# (b) tBP
axs[1].plot(stats_df["end_year"], stats_df["tBP"], color="green", label="tBP")
axs[1].scatter(years_t, stats_df["tBP"].iloc[peaks_t], color="red")
for yr in years_t:
axs[1].text(yr, stats_df.loc[stats_df.end_year==yr, "tBP"].values[0],
str(int(yr)), ha="center", va="bottom", fontsize=8)
axs[1].axhline(3.5, color="gray", ls="--")
axs[1].set_ylabel("tBP")
axs[1].set_title("tBP Peaks")
axs[1].grid(True)
# (c) GLK
axs[2].plot(stats_df["end_year"], stats_df["GLK"], color="magenta", label="GLK")
axs[2].scatter(years_g, stats_df["GLK"].iloc[peaks_g], color="red")
for yr in years_g:
axs[2].text(yr, stats_df.loc[stats_df.end_year==yr, "GLK"].values[0],
str(int(yr)), ha="center", va="bottom", fontsize=8)
axs[2].axhline(0.6, color="gray", ls="--")
axs[2].set_ylabel("GLK")
axs[2].set_xlabel("End Year")
axs[2].set_title("GLK Peaks")
axs[2].grid(True)
plt.tight_layout()
plt.show()
In [ ]: