年輪年代法による年代推定¶

COFECHA の既定曲線によるデトレンド¶

パブリックドメインに公開します、改造等はご自由に

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()
No description has been provided for this image
  • “本当に一致” と判断するには、同じ一致終了年で r 高・tBP ≧3.5・GLK ≧0.60 が揃うのが理想
  • 591 年にて全ての検定を有意に満たしており、法隆寺五重塔心柱の最外周は 591 年であると推定される