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

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/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()
No description has been provided for this image
In [ ]: