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/130¥011-011/011-011.p4p'
ser_unknown = load_p4p_series(filepath, "unknown", 0)
print(ser_unknown[:5])
print(ser_unknown[-5:])
Year 0 251.0 1 354.0 2 390.0 3 291.0 4 322.0 Name: unknown, dtype: float64 Year 347 76.0 348 66.0 349 71.0 350 45.0 351 60.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, start_year_known)
display(stats_df.head(3))
display(stats_df.tail(3))
| end_year | r | tBP | GLK | |
|---|---|---|---|---|
| 0 | 314 | 0.004680 | 0.087551 | 0.541311 |
| 1 | 315 | -0.000319 | -0.005966 | 0.455840 |
| 2 | 316 | 0.010253 | 0.191827 | 0.532764 |
| end_year | r | tBP | GLK | |
|---|---|---|---|---|
| 521 | 835 | -0.033827 | -0.633210 | 0.472934 |
| 522 | 836 | 0.052526 | 0.984022 | 0.515670 |
| 523 | 837 | -0.013114 | -0.245366 | 0.467236 |
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 のピーク年: [320, 325, 336, 339, 348, 365, 374, 388, 391, 396, 404, 420, 428, 463, 465, 479, 485, 488, 517, 522, 528, 533, 552, 566, 572, 577, 583, 588, 591, 605, 616, 620, 631, 633, 636, 645, 648, 658, 660, 664, 666, 675, 680, 684, 692, 700, 706, 712, 727, 731, 737, 751, 759, 762, 769, 771, 773, 775, 779, 783, 797, 805, 808, 815, 820, 822, 828] tBP のピーク年: [320, 323, 325, 328, 331, 334, 336, 339, 341, 345, 348, 350, 356, 360, 365, 369, 372, 374, 377, 379, 382, 384, 388, 391, 393, 396, 400, 402, 404, 407, 410, 416, 418, 420, 425, 428, 432, 434, 436, 440, 443, 445, 448, 454, 456, 459, 463, 465, 467, 470, 472, 474, 476, 479, 482, 485, 488, 492, 495, 501, 503, 511, 517, 522, 524, 526, 528, 533, 540, 544, 547, 552, 557, 561, 563, 566, 568, 572, 577, 580, 583, 586, 588, 591, 593, 595, 597, 599, 605, 610, 612, 616, 620, 625, 627, 631, 633, 636, 640, 642, 645, 648, 650, 655, 658, 660, 662, 664, 666, 668, 670, 672, 675, 680, 684, 689, 692, 695, 700, 704, 706, 709, 712, 714, 720, 727, 729, 731, 733, 737, 739, 746, 751, 755, 757, 759, 762, 765, 767, 769, 771, 773, 775, 779, 783, 786, 790, 793, 795, 797, 800, 803, 805, 808, 815, 820, 822, 824, 828, 830, 834, 836] GLK のピーク年: [316, 323, 325, 336, 338, 341, 343, 350, 354, 361, 365, 370, 376, 382, 384, 387, 391, 396, 402, 407, 410, 420, 424, 427, 450, 452, 456, 463, 467, 474, 482, 485, 488, 498, 515, 517, 519, 526, 534, 540, 544, 549, 551, 553, 557, 566, 572, 583, 588, 591, 595, 597, 602, 612, 620, 627, 631, 633, 636, 638, 640, 642, 653, 658, 660, 662, 664, 666, 668, 679, 684, 686, 700, 702, 704, 709, 711, 717, 719, 721, 729, 739, 741, 744, 748, 755, 757, 759, 762, 766, 769, 773, 780, 786, 788, 795, 803, 808, 810, 812, 820, 822, 824, 828, 830]
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()
- “本当に一致” と判断するには、同じ一致終了年で r 高・tBP ≧3.5・GLK ≧0.60 が揃うのが理想
- 591 年にて全ての検定を有意に満たしており、法隆寺五重塔心柱の最外周は 591 年であると推定される