D.1 Pythonコード全文
ファイル名: analyze_ancient_japan.py
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
analyze_ancient_japan.py
────────────────────────────────────────────────────────────────────────────
物部系神社(n=17)・百名山神社（n=18）・一之宮（n=51）と古代鉱山（n=27）の
空間統計分析スクリプト
【分析内容】
1. Haversine式による大圏距離計算
2. 記述統計（平均・中央値・標準偏差・30km圏内率など）
3. Mann-Whitney U検定（ノンパラメトリック）
4. Cohen's d 効果量
5. Bonferroni補正による多重検定評価
6. 検出力分析（Post-hoc power analysis）
7. 距離分布ヒストグラム、箱ひげ図、CDF、ランキング図、閾値分析図の生成
【対照群の選定方針】
百名山神社群（n=18）: 深田久弥『日本百名山』（1964年）に選定された山岳のうち
福島以南の本州に所在し、山岳信仰を示す神社が現存するものを機械的全件抽出。
箱根神社（箱根山）・戸隠神社（戸隠山）はこの基準を満たさないため含めない。
【入力ファイル】
data/mononobe_shrines_honshu.csv : 物部系神社 n=17
data/hyakumeizan_shrines.csv : 百名山神社 n=18
data/ichinomiya_shrines.csv : 一之宮 n=51
data/mines_ancient_honshu.csv : 古代鉱山 n=27
【出力】
images/0_histograms_distribution.png : 距離分布ヒストグラム（三群）
images/1_threshold_analysis.png : 距離閾値別カバレッジ
images/2_boxplot_detailed.png : 詳細箱ひげ図（三群）
images/3_mononobe_ranking.png : 物部系神社17社 距離ランキング
images/4_cdf_analysis.png : 累積分布関数（CDF）
────────────────────────────────────────────────────────────────────────────
"""
import csv
import math
import os
from typing import List, Tuple
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm
import matplotlib.patches as mpatches
import numpy as np
from scipy import stats
# ── フォント設定（日本語対応） ────────────────────────────────────────────
_jp_font = None
for name in ["Noto Sans CJK JP", "Noto Serif CJK JP", "IPAGothic", "IPAPGothic"]:
try:
fm.findfont(fm.FontProperties(family=name), fallback_to_default=False)
_jp_font = name
break
except Exception:
pass
if _jp_font:
plt.rcParams["font.family"] = _jp_font
else:
plt.rcParams["font.family"] = "DejaVu Sans"
plt.rcParams["axes.unicode_minus"] = False
plt.rcParams["figure.dpi"] = 150
# ── 定数 ─────────────────────────────────────────────────────────────────
EARTH_R = 6371.0 # km
THRESHOLD_KM = 30 # 実効支配距離の閾値
COLOR_M = "#c0392b" # 物部系神社（赤）
COLOR_H = "#2980b9" # 百名山神社（青）
COLOR_I = "#27ae60" # 一之宮（緑）
ALPHA = 0.01 # 有意水準
ALPHA_FWER = 0.05 # 族全体有意水準（Bonferroni補正の出発点）
BONFERRONI_N = 2 # 検定回数
ALPHA_BONF = ALPHA_FWER / BONFERRONI_N # Bonferroni補正後 α' = 0.025
# ════════════════════════════════════════════════════════════════════════════
# 1. ユーティリティ関数
# ════════════════════════════════════════════════════════════════════════════
def haversine(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
"""Haversine式による大圏距離（km）"""
lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])
dlat = lat2 - lat1
dlon = lon2 - lon1
a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
return 2 * math.asin(math.sqrt(a)) * EARTH_R
def nearest_mine_distance(shrine_lat: float, shrine_lon: float,
mines: List[Tuple]) -> Tuple[float, str]:
"""最近接鉱山までの距離と鉱山名を返す"""
min_d = float("inf")
min_name = ""
for m in mines:
d = haversine(shrine_lat, shrine_lon, m["lat"], m["lon"])
if d < min_d:
min_d = d
min_name = m["name"]
return min_d, min_name
def load_csv(path: str) -> List[dict]:
with open(path, newline="", encoding="utf-8") as f:
return list(csv.DictReader(f))
def cohen_d(x: np.ndarray, y: np.ndarray) -> float:
"""Pooled SD による Cohen's d"""
nx, ny = len(x), len(y)
pooled_std = math.sqrt(((nx - 1) * np.std(x, ddof=1) ** 2 +
(ny - 1) * np.std(y, ddof=1) ** 2) / (nx + ny - 2))
return abs(np.mean(x) - np.mean(y)) / pooled_std if pooled_std > 0 else 0.0
def effect_label(d: float) -> str:
if d >= 0.8:
return "Large"
elif d >= 0.5:
return "Medium"
elif d >= 0.2:
return "Small"
return "Negligible"
def post_hoc_power(d: float, n1: int, n2: int, alpha: float = 0.01) -> float:
"""
Two-sample t-test の事後検出力（Mann-Whitney の近似として利用）
Cohen's d と両サンプルサイズから計算。
"""
se = math.sqrt(1 / n1 + 1 / n2)
ncp = d / se # non-centrality parameter
z_alpha = stats.norm.ppf(1 - alpha / 2)
return float(stats.norm.cdf(ncp - z_alpha) + stats.norm.cdf(-ncp - z_alpha))
# ════════════════════════════════════════════════════════════════════════════
# 2. データ読み込みと距離計算
# ════════════════════════════════════════════════════════════════════════════
BASE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(BASE, "data")
IMG = os.path.join(BASE, "images")
os.makedirs(IMG, exist_ok=True)
print("=" * 60)
print(" 古代日本 空間統計分析")
print(" 物部神社・百名山神社・一之宮 vs 古代鉱山")
print("=" * 60)
# --- 読み込み
raw_mononobe = load_csv(os.path.join(DATA, "mononobe_shrines_honshu.csv"))
raw_hyakumeizan = load_csv(os.path.join(DATA, "hyakumeizan_shrines.csv"))
raw_ichinomiya = load_csv(os.path.join(DATA, "ichinomiya_shrines_full.csv"))
raw_mines = load_csv(os.path.join(DATA, "mines_ancient_honshu.csv"))
mines = [{"name": m["name"], "lat": float(m["lat"]), "lon": float(m["lon"])}
for m in raw_mines]
# --- 距離計算
def calc_distances(shrines, mines, label):
results = []
for s in shrines:
d, mine_name = nearest_mine_distance(float(s["lat"]), float(s["lon"]), mines)
results.append({
"shrine": s["name"],
"lat": float(s["lat"]),
"lon": float(s["lon"]),
"dist_km": round(d, 2),
"nearest_mine": mine_name,
})
dists = np.array([r["dist_km"] for r in results])
return results, dists
mononobe_data, d_mono = calc_distances(raw_mononobe, mines, "物部系神社")
hyakumeizan_data, d_hyaku = calc_distances(raw_hyakumeizan, mines, "百名山神社")
ichinomiya_data, d_ichi = calc_distances(raw_ichinomiya, mines, "一之宮")
print(f"\nデータ確認:")
print(f" 物部系神社 : {len(d_mono)}社")
print(f" 百名山神社 : {len(d_hyaku)}社")
print(f" 一之宮 : {len(d_ichi)}社")
print(f" 古代鉱山 : {len(mines)}鉱山")
# ════════════════════════════════════════════════════════════════════════════
# 3. 記述統計
# ════════════════════════════════════════════════════════════════════════════
def describe(arr: np.ndarray, label: str) -> dict:
n = len(arr)
w30 = (arr <= THRESHOLD_KM).sum()
w10 = (arr <= 10).sum()
return {
"label": label,
"n": n,
"mean": round(arr.mean(), 2),
"median": round(float(np.median(arr)), 2),
"std": round(arr.std(ddof=1), 2),
"min": round(arr.min(), 2),
"max": round(arr.max(), 2),
"within10_n": int(w10),
"within10_pct": round(w10 / n * 100, 1),
"within30_n": int(w30),
"within30_pct": round(w30 / n * 100, 1),
}
stats_mono = describe(d_mono, "物部系神社 (n=17)")
stats_hyaku = describe(d_hyaku, "百名山神社 (n=18)")
stats_ichi = describe(d_ichi, "一之宮 (n=51)")
print("\n" + "─" * 60)
print("【記述統計】")
print("─" * 60)
header = f"{'統計量':<16} {'物部系':>10} {'百名山':>10} {'一之宮':>10}"
print(header)
print("─" * 50)
rows = [
("n", stats_mono["n"], stats_hyaku["n"], stats_ichi["n"]),
("平均 (km)", stats_mono["mean"], stats_hyaku["mean"], stats_ichi["mean"]),
("中央値 (km)",stats_mono["median"], stats_hyaku["median"], stats_ichi["median"]),
("標準偏差", stats_mono["std"], stats_hyaku["std"], stats_ichi["std"]),
("最小値", stats_mono["min"], stats_hyaku["min"], stats_ichi["min"]),
("最大値", stats_mono["max"], stats_hyaku["max"], stats_ichi["max"]),
("10km以内数", stats_mono["within10_n"], stats_hyaku["within10_n"], stats_ichi["within10_n"]),
("10km以内%", stats_mono["within10_pct"], stats_hyaku["within10_pct"], stats_ichi["within10_pct"]),
("30km以内数", stats_mono["within30_n"], stats_hyaku["within30_n"], stats_ichi["within30_n"]),
("30km以内%", stats_mono["within30_pct"], stats_hyaku["within30_pct"], stats_ichi["within30_pct"]),
]
for name, v1, v2, v3 in rows:
print(f" {name:<14} {str(v1):>10} {str(v2):>10} {str(v3):>10}")
# ════════════════════════════════════════════════════════════════════════════
# 4. 推測統計 — Mann-Whitney U 検定 + Cohen's d
# ════════════════════════════════════════════════════════════════════════════
print("\n" + "─" * 60)
print("【推測統計】")
print("─" * 60)
# 分析1: 物部 vs 百名山神社
u1, p1 = stats.mannwhitneyu(d_mono, d_hyaku, alternative="two-sided")
d1 = cohen_d(d_mono, d_hyaku)
pw1 = post_hoc_power(d1, len(d_mono), len(d_hyaku), alpha=ALPHA)
print(f"\n ■ 物部系神社 (n={len(d_mono)}) vs 百名山神社 (n={len(d_hyaku)})")
print(f" Mann-Whitney U統計量 : {u1:.1f}")
print(f" p値 (両側検定): {p1:.6f}")
print(f" Cohen's d : {d1:.4f} [{effect_label(d1)}]")
print(f" 有意水準 α={ALPHA} : {'有意 ✓' if p1 < ALPHA else '非有意'}")
print(f" Bonferroni α'={ALPHA_BONF}: {'有意 ✓' if p1 < ALPHA_BONF else '非有意'}")
print(f" 事後検出力 : {pw1:.3f} ({pw1*100:.1f}%)")
# 分析2: 物部 vs 一之宮
u2, p2 = stats.mannwhitneyu(d_mono, d_ichi, alternative="two-sided")
d2 = cohen_d(d_mono, d_ichi)
pw2 = post_hoc_power(d2, len(d_mono), len(d_ichi), alpha=ALPHA)
print(f"\n ■ 物部系神社 (n={len(d_mono)}) vs 一之宮 (n={len(d_ichi)})")
print(f" Mann-Whitney U統計量 : {u2:.1f}")
print(f" p値 (両側検定) : {p2:.6f}")
print(f" Cohen's d : {d2:.4f} [{effect_label(d2)}]")
print(f" 有意水準 α={ALPHA} : {'有意 ✓' if p2 < ALPHA else '非有意'}")
print(f" Bonferroni α'={ALPHA_BONF}: {'有意 ✓' if p2 < ALPHA_BONF else '非有意'}")
print(f" 事後検出力 : {pw2:.3f} ({pw2*100:.1f}%)")
# 感度分析: n=14 (物部社名のみ) vs 百名山神社 / 一之宮
exclude_extended = ["大津神社（飛騨市）", "彌彦神社（弥彦村）", "石上神宮"]
d_mono14 = np.array([
mononobe_data[i]["dist_km"] for i, s in enumerate(raw_mononobe)
if s["name"] not in exclude_extended
])
u1b, p1b = stats.mannwhitneyu(d_mono14, d_hyaku, alternative="two-sided")
d1b = cohen_d(d_mono14, d_hyaku)
u2b, p2b = stats.mannwhitneyu(d_mono14, d_ichi, alternative="two-sided")
d2b = cohen_d(d_mono14, d_ichi)
print(f"\n ■ 感度分析: 物部社名のみ (n={len(d_mono14)}) vs 百名山神社 (n={len(d_hyaku)})")
print(f" p値: {p1b:.6f} Cohen's d: {d1b:.4f} [{effect_label(d1b)}]")
print(f" 有意 α={ALPHA}: {'✓' if p1b < ALPHA else '✗ (α=0.05: ' + ('✓' if p1b < 0.05 else '✗') + ')'}")
print(f"\n ■ 感度分析: 物部社名のみ (n={len(d_mono14)}) vs 一之宮 (n={len(d_ichi)})")
print(f" p値: {p2b:.6f} Cohen's d: {d2b:.4f} [{effect_label(d2b)}]")
print(f" 有意 α={ALPHA}: {'✓' if p2b < ALPHA else '✗'}")
# 正規性検定
_, p_shapiro_mono = stats.shapiro(d_mono)
_, p_shapiro_hyaku = stats.shapiro(d_hyaku)
_, p_shapiro_ichi = stats.shapiro(d_ichi)
print(f"\n ■ Shapiro-Wilk正規性検定")
print(f" 物部系: p={p_shapiro_mono:.4f} {'正規' if p_shapiro_mono > 0.05 else '非正規'}")
print(f" 百名山: p={p_shapiro_hyaku:.4f} {'正規' if p_shapiro_hyaku > 0.05 else '非正規'}")
print(f" 一之宮: p={p_shapiro_ichi:.4f} {'正規' if p_shapiro_ichi > 0.05 else '非正規'}")
# 等分散性検定
lev_stat, lev_p = stats.levene(d_mono, d_hyaku)
print(f"\n ■ Levene等分散性検定 (物部 vs 百名山)")
print(f" F={lev_stat:.3f}, p={lev_p:.4f} {'等分散' if lev_p > 0.05 else '非等分散'}")
# ════════════════════════════════════════════════════════════════════════════
# 5. グラフ描画
# ════════════════════════════════════════════════════════════════════════════
LABEL_M = f"物部系神社 (n={len(d_mono)})"
LABEL_H = f"百名山神社 (n={len(d_hyaku)})"
LABEL_I = f"一之宮 (n={len(d_ichi)})"
# ─── 図0：距離分布ヒストグラム（三群） ─────────────────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharey=False)
fig.suptitle("三群の最近接鉱山距離分布", fontsize=14, fontweight="bold", y=1.01)
bins = np.arange(0, 145, 10)
configs = [
(d_mono, COLOR_M, LABEL_M),
(d_hyaku, COLOR_H, LABEL_H),
(d_ichi, COLOR_I, LABEL_I),
]
for ax, (data, color, label) in zip(axes, configs):
ax.hist(data, bins=bins, color=color, alpha=0.75, edgecolor="white", linewidth=0.8)
med = np.median(data)
ax.axvline(med, color="#333", linestyle="--", linewidth=1.5,
label=f"中央値: {med:.1f} km")
ax.axvline(THRESHOLD_KM, color="gray", linestyle=":", linewidth=1.2, alpha=0.8)
ax.set_title(label, fontsize=11, pad=8)
ax.set_xlabel("最近接鉱山距離 (km)", fontsize=10)
ax.set_ylabel("神社数", fontsize=10)
ax.legend(fontsize=9)
ax.set_xlim(0, 140)
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
path0 = os.path.join(IMG, "0_histograms_distribution.png")
plt.savefig(path0, dpi=150, bbox_inches="tight")
plt.close()
print(f"\n → 保存: {path0}")
# ─── 図1：距離閾値別カバレッジ ───────────────────────────────────────────
thresholds = [10, 20, 30, 50, 100]
cov_mono = [(d_mono <= t).mean() * 100 for t in thresholds]
cov_hyaku = [(d_hyaku <= t).mean() * 100 for t in thresholds]
cov_ichi = [(d_ichi <= t).mean() * 100 for t in thresholds]
fig, ax = plt.subplots(figsize=(9, 6))
x = np.arange(len(thresholds))
w = 0.26
rects1 = ax.bar(x - w, cov_mono, w, label=LABEL_M, color=COLOR_M, alpha=0.85)
rects2 = ax.bar(x, cov_hyaku, w, label=LABEL_H, color=COLOR_H, alpha=0.85)
rects3 = ax.bar(x + w, cov_ichi, w, label=LABEL_I, color=COLOR_I, alpha=0.85)
for rects in [rects1, rects2, rects3]:
for rect in rects:
h = rect.get_height()
ax.annotate(f"{h:.0f}%",
xy=(rect.get_x() + rect.get_width() / 2, h),
xytext=(0, 3), textcoords="offset points",
ha="center", va="bottom", fontsize=8)
ax.set_xticks(x)
ax.set_xticklabels([f"{t} km" for t in thresholds])
ax.set_xlabel("距離閾値", fontsize=11)
ax.set_ylabel("閾値以内の神社の割合 (%)", fontsize=11)
ax.set_title("距離閾値別・鉱山圏内神社の割合（三群比較）", fontsize=13, fontweight="bold")
ax.legend(fontsize=10)
ax.set_ylim(0, 110)
ax.grid(axis="y", alpha=0.3)
fig.tight_layout()
path1 = os.path.join(IMG, "1_threshold_analysis.png")
plt.savefig(path1, dpi=150, bbox_inches="tight")
plt.close()
print(f" → 保存: {path1}")
# ─── 図2：詳細箱ひげ図 ───────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(10, 7))
bp = ax.boxplot(
[d_mono, d_hyaku, d_ichi],
patch_artist=True,
notch=False,
widths=0.5,
medianprops=dict(color="#1a3a6b", linewidth=2.5),
whiskerprops=dict(linewidth=1.5),
capprops=dict(linewidth=1.5),
flierprops=dict(marker="o", markersize=5, alpha=0.6),
)
colors = [COLOR_M, COLOR_H, COLOR_I]
for patch, color in zip(bp["boxes"], colors):
patch.set_facecolor(color)
patch.set_alpha(0.65)
# 平均値プロット
means = [d_mono.mean(), d_hyaku.mean(), d_ichi.mean()]
ax.scatter([1, 2, 3], means, color="red", zorder=5, s=60, label="平均値", marker="D")
ax.axhline(THRESHOLD_KM, color="gray", linestyle="--", linewidth=1.3, alpha=0.7,
label=f"{THRESHOLD_KM} km 閾値")
ax.set_xticks([1, 2, 3])
ax.set_xticklabels([LABEL_M, LABEL_H, LABEL_I], fontsize=10)
ax.set_ylabel("最近接鉱山距離 (km)", fontsize=11)
ax.set_title("三群の鉱山距離分布（箱ひげ図）", fontsize=13, fontweight="bold")
ax.legend(fontsize=10)
ax.grid(axis="y", alpha=0.3)
# p値注記
y_max = max(d_mono.max(), d_hyaku.max(), d_ichi.max()) * 1.05
ax.annotate(
f"p={p1:.4f} (Mann-Whitney)\nd={d1:.3f} [{effect_label(d1)}]",
xy=(1.5, y_max * 0.88), ha="center", fontsize=9,
bbox=dict(boxstyle="round,pad=0.3", facecolor="white", edgecolor=COLOR_M, alpha=0.8),
)
fig.tight_layout()
path2 = os.path.join(IMG, "2_boxplot_detailed.png")
plt.savefig(path2, dpi=150, bbox_inches="tight")
plt.close()
print(f" → 保存: {path2}")
# ─── 図3：物部系神社17社 ランキング ─────────────────────────────────────
sorted_mono = sorted(mononobe_data, key=lambda x: x["dist_km"])
names = [s["shrine"] for s in sorted_mono]
dists = [s["dist_km"] for s in sorted_mono]
mines_ = [s["nearest_mine"] for s in sorted_mono]
colors_bar = [COLOR_M if d <= THRESHOLD_KM else "#e8a0a0" for d in dists]
fig, ax = plt.subplots(figsize=(13, 7))
bars = ax.barh(range(len(names)), dists, color=colors_bar, alpha=0.85, edgecolor="white")
ax.axvline(THRESHOLD_KM, color="#333", linestyle="--", linewidth=1.8,
label=f"{THRESHOLD_KM} km 閾値")
for i, (d, mine) in enumerate(zip(dists, mines_)):
ax.text(d + 0.8, i, f"{d:.1f} km [{mine}]", va="center", fontsize=8.5, color="#222")
ax.set_yticks(range(len(names)))
ax.set_yticklabels(names, fontsize=9)
ax.set_xlabel("最近接鉱山距離 (km)", fontsize=11)
ax.set_title("物部系神社17社 — 最近接鉱山距離（昇順）", fontsize=13, fontweight="bold")
within30 = sum(1 for d in dists if d <= THRESHOLD_KM)
pct30 = within30 / len(dists) * 100
patch1 = mpatches.Patch(color=COLOR_M, alpha=0.85, label=f"30km以内 ({within30}社, {pct30:.1f}%)")
patch2 = mpatches.Patch(color="#e8a0a0", alpha=0.85, label=f"30km超 ({len(dists)-within30}社)")
ax.legend(handles=[patch1, patch2], fontsize=10, loc="lower right")
ax.set_xlim(0, max(dists) * 1.35)
ax.grid(axis="x", alpha=0.3)
ax.invert_yaxis()
fig.tight_layout()
path3 = os.path.join(IMG, "3_mononobe_ranking.png")
plt.savefig(path3, dpi=150, bbox_inches="tight")
plt.close()
print(f" → 保存: {path3}")
# ─── 図4：累積分布関数（CDF） ────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(10, 7))
for data, color, label in [(d_mono, COLOR_M, LABEL_M),
(d_hyaku, COLOR_H, LABEL_H),
(d_ichi, COLOR_I, LABEL_I)]:
xs = np.sort(data)
ys = np.arange(1, len(xs) + 1) / len(xs) * 100
ax.step(xs, ys, where="post", color=color, linewidth=2.2, label=label)
ax.axvline(THRESHOLD_KM, color="#555", linestyle="--", linewidth=1.5, alpha=0.8,
label=f"{THRESHOLD_KM} km 閾値")
# 30km地点の読み取り値注記
for data, color in [(d_mono, COLOR_M), (d_hyaku, COLOR_H), (d_ichi, COLOR_I)]:
y30 = (data <= THRESHOLD_KM).mean() * 100
ax.scatter([THRESHOLD_KM], [y30], color=color, zorder=5, s=55)
ax.annotate(f"{y30:.1f}%", xy=(THRESHOLD_KM, y30),
xytext=(2.5, 2.5), textcoords="offset points",
color=color, fontsize=9, fontweight="bold")
ax.set_xlabel("最近接鉱山距離 (km)", fontsize=11)
ax.set_ylabel("神社の累積割合 (%)", fontsize=11)
ax.set_title("三群の鉱山距離 累積分布関数（CDF）", fontsize=13, fontweight="bold")
ax.legend(fontsize=10)
ax.set_xlim(0, 140)
ax.set_ylim(0, 105)
ax.grid(alpha=0.3)
fig.tight_layout()
path4 = os.path.join(IMG, "4_cdf_analysis.png")
plt.savefig(path4, dpi=150, bbox_inches="tight")
plt.close()
print(f" → 保存: {path4}")
# ════════════════════════════════════════════════════════════════════════════
# 6. 最終サマリー
# ════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 60)
print("【最終統計サマリー】")
print("=" * 60)
print(f"\n 物部系神社 (n=17): 平均 {d_mono.mean():.2f}km / 中央値 {np.median(d_mono):.2f}km")
print(f"\n 対 百名山神社 (n={len(d_hyaku)}):")
print(f" Mann-Whitney U検定 p値: {p1:.6f}")
print(f" 効果量 (Cohen's d) : {d1:.4f} [{effect_label(d1)}]")
print(f" U統計量 : {u1:.1f}")
print(f"\n 対 一之宮 (n={len(d_ichi)}):")
print(f" Mann-Whitney U検定 p値: {p2:.6f}")
print(f" 効果量 (Cohen's d) : {d2:.4f} [{effect_label(d2)}]")
print(f" U統計量 : {u2:.1f}")
print(f"\n 30km圏内到達率:")
print(f" 物部系神社: {stats_mono['within30_pct']}% ({stats_mono['within30_n']}/{stats_mono['n']})")
print(f" 百名山神社: {stats_hyaku['within30_pct']}% ({stats_hyaku['within30_n']}/{stats_hyaku['n']})")
print(f" 一之宮 : {stats_ichi['within30_pct']}% ({stats_ichi['within30_n']}/{stats_ichi['n']})")
print()
print(" 生成画像:")
for p in [path0, path1, path2, path3, path4]:
print(f" {os.path.basename(p)}")
print("\n 完了。")
付録E 春成モデルのベイズ統計適合度（Amodel）の独立検証
証拠レベル1：OxCalによる客観的再現結果
【方法論的前提：なぜ春成論文のサンプルデータを使って春成モデルを批判できるのか】
本稿は春成モデルを批判しながら、その春成ら（2011）論文に由来するC14実測値（サンプルデータ）を再検証に使用している。これは矛盾ではない——批判対象と使用データは、以下のように明確に区別される。
✘ 批判の対象
①Prior（事前制約）——春成の編年順序の設定
②較正曲線の選択——検証不能なJCAL
③サンプル操作——ホケノ山小枝の除外
✔ 信頼して使用するもの
C14実測値そのもの——加速器質量分析計（AMS）が出力した物理的測定値。測定プロセスは国際標準に従っており、この数値を批判する根拠はない。
科学的反証の正しい形式は「同じ生データを、より適切なモデルで再計算して異なる結論を示すこと」である。異なるデータを使えば「別の話」になり反証にならない。本付録が春成ら（2011）のC14実測値をそのまま使い、PriorとIntCal20だけを変えて再計算するのは、反証方法論として最も論理的に正当な手続きである。
本稿の第8章において、坂本（2022）が報告した春成モデルの統計的適合度の低さ（Amodel=16%等）を引用して批判を行った。この統計的棄却が特定の研究者による恣意的な操作ではなく、公開データに基づく数学的・客観的な帰結であることを証明するため、著者自身が最新の国際標準較正曲線（IntCal20）とオックスフォード大学のOxCal（v4.4.4）を用いて独立した再検証を実施した。
E.1 春成モデル使用試料の全一覧（56点・属性付き）
春成ら（2011）付表および坂本（2022）に基づき著者が整理。試料種別・属性は原文付表の記載を優先。リスク区分は春成モデルの¹⁴C年代指標としての信頼性評価（高=古木効果・リザーバー効果の危険あり、中=一定の対策あり、低=単年性・高信頼性）。⚠は本稿追試でA値低下またはPoor agreementを確認した試料。
試料の信頼性分布サマリー（春成モデル採用56点）
23.2%
低リスク高信頼性
（低+低* 13/56点）
低リスク高信頼性試料は23.2%（13点）存在するが、その2倍超の高リスク試料（33点）に統計的に支配されている。さらに唯一の最高信頼性アンカー（ホケノ山小枝）が除外され、モデルに方向性が生じない構造となっている。
高：古木効果・リザーバー効果の危険
中：古木効果の可能性（辺材採取等で軽減）・吹（リザーバー効果の留保あり）
低：単年性・高信頼性
低*：材質は高信頼性だが層位撹乱疑い
| 試料ID |
フェーズ |
遺跡・古墳 |
試料種別 |
属性 |
リスク |
¹⁴C（BP） |
± |
採否 |
| Phase 1 大和VI-4様式（大福遺跡） |
| NRSK-65 | Phase 1 | 大福遺跡 | 土器付着物 | 煤（甕胴中外面） | 高 | 1910 | ±20 | 採用 |
| NRSK-66 | Phase 1 | 大福遺跡 | 土器付着物 | 煤（甕胴上外面） | 高 | 1920 | ±15 | 採用 |
| NRSK-C63 | Phase 1 | 大福遺跡 | 種実 | コメ（炭化米） | 低 | 1925 | ±20 | 採用 |
| Phase 2 庄内0式期（纒向遺跡） |
| NRSK-21 | Phase 2 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面） | 高 | 1920 | ±25 | 採用 |
| NRSK-22_b | Phase 2 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面） | 高 | 1910 | ±25 | 採用 |
| NRSK-23_b | Phase 2 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面） | 高 | 1915 | ±25 | 採用 |
| NRSK-33_a | Phase 2 | 大福遺跡26次 | 土器付着物 | 焦（甕胴内面） | 高 | 1920 | ±25 | 採用 |
| NRSK-33_b | Phase 2 | 大福遺跡26次 | 土器付着物 | 煤（甕胴外面） | 高 | 1965 | ±25 | 採用 |
| NRSK-34_b | Phase 2 | 大福遺跡26次 | 土器付着物 | 煤（甕胴外面） | 高 | 1905 | ±25 | 採用 |
| Phase 3 庄内1式期（纒向遺跡） |
| NRSK-24 | Phase 3 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面） | 高 | 1930 | ±25 | 採用 |
| NRSK-25_a | Phase 3 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面） | 高 | 1925 | ±25 | 採用 |
| NRSK-25_b | Phase 3 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面） | 高 | 1910 | ±25 | 採用 |
| NRSK-26 | Phase 3 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面） | 高 | 1910 | ±25 | 採用 |
| Phase 4 庄内3式期（纒向遺跡） |
| NRSK-27 | Phase 4 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面付着） | 高 | 1890 | ±25 | 採用 |
| NRSK-29_b | Phase 4 | 纒向遺跡 | 土器付着物 | 煤（甕胴外面付着） | 高 | 1865 | ±25 | 採用 |
| Phase 5 纒向石塚墳丘墓期（纒向石塚・矢塚古墳） |
| NRSK-C14 | Phase 5 | 纒向石塚古墳 | 木材 | 小枝 | 低 | 1830 | ±30 | 採用 |
| NRSK-C16 | Phase 5 | 纒向石塚古墳 | 種実 | 植物遺体（種実類） | 低 | 1910 | ±30 | 採用 |
| NRSK-C17 | Phase 5 | 纒向石塚古墳 | 木材 | ヒノキ（植物層） | 中 | 1890 | ±30 | 採用 |
| NRSK-C18 | Phase 5 | 纒向石塚古墳 | 木材 | 不明環孔材（墳丘下湿地層） | 中 | 1880 | ±30 | 採用 |
| NRSK-6 | Phase 5 | 矢塚古墳 | 土器付着物 | 煤（甕胴外面） | 高 | 1820 | ±30 | 採用 |
| Phase 6a 東田大塚墳丘墓築造前 |
| NRSK-2 | Phase 6a | 東田大塚古墳 | 土器付着物 | 煤（甕胴外面） | 高 | 1860 | ±30 | 採用 |
| NRSK-3_a | Phase 6a | 東田大塚古墳 | 土器付着物 | 焦（甕胴内面） | 高 | 1820 | ±30 | 採用 |
| NRSK-3_b | Phase 6a | 東田大塚古墳 | 土器付着物 | 煤（甕胴外面） | 高 | 1780 | ±30 | 採用 |
| NRSK-C1 | Phase 6a | 東田大塚古墳 | 種実 | ウリ科 種子 | 低 | 1850 | ±30 | 採用 |
| NRSK-C2 | Phase 6a | 東田大塚古墳 | 種実 | モモ 核 | 低* | 1730 | ±30 | 採用⚠ A=5.7% |
| Phase 6b 東田大塚築造中 |
| NRSK-C7 | Phase 6b | 東田大塚古墳 | その他 | タケ亜科イネ科 カゴ材 | 高 | 1760 | ±30 | 採用⚠ A=37.5% |
| NRSK-C8 | Phase 6b | 東田大塚古墳 | その他 | タケ亜科イネ科 カゴ材 | 高 | 1730 | ±30 | 採用⚠ A=6.9% |
| Phase 6c 箸墓古墳築造前（箸墓古墳7次調査・県調査SX01） ※このフェーズのIDプレフィックスはNRKS |
| NRKS-02 | Phase 6c | 箸墓古墳周辺 | 土器付着物 | 煤（口縁外面） | 高 | 1840 | ±30 | 採用 |
| NRKS-03 | Phase 6c | 箸墓古墳周辺 | 土器付着物 | 煤（口縁外面） | 高 | 1780 | ±30 | 採用 |
| NRKS-04 | Phase 6c | 箸墓古墳周辺 | 土器付着物 | 吹（口縁外面） | 中 | 1820 | ±30 | 採用 |
| NRKS-06 | Phase 6c | 箸墓古墳周辺 | 土器付着物 | 吹（甕F口縁外面） | 中 | 1830 | ±40 | 採用 |
| NRKS-09 | Phase 6c | 箸墓古墳周辺 | 土器付着物 | 煤（甕SY口縁外面） | 高 | 1910 | ±40 | 採用⚠ Poor |
| NRKS-10 | Phase 6c | 箸墓古墳周辺 | 土器付着物 | 煤（甕S口縁外面） | 高 | 1840 | ±40 | 採用 |
| NRKS-11 | Phase 6c | 箸墓古墳周辺 | 土器付着物 | 煤（甕F口縁外面） | 高 | 1740 | ±40 | 採用⚠ Poor |
| NRKS-13 | Phase 6c | 箸墓古墳周辺 | 土器付着物 | 煤（甕口縁外面） | 高 | 1870 | ±40 | 採用⚠ Poor |
| Phase 7 布留0式新相（箸墓古墳期） |
| NRSK-1 | Phase 7 | 箸墓古墳周辺 | 土器付着物 | 煤（甕胴外面） | 高 | 1710 | ±30 | 採用 |
| NRSK-4_a | Phase 7 | 箸墓古墳周辺 | 土器付着物 | 焦（甕胴内面） | 高 | 1840 | ±30 | 採用⚠ Poor |
| NRSK-4_b | Phase 7 | 箸墓古墳周辺 | 土器付着物 | 煤（甕胴外面） | 高 | 1750 | ±30 | 採用⚠ Poor |
| NRSK-C4 | Phase 7 | 箸墓古墳周辺 | 木材 | エノキ属 加工木（樹皮直下） | 中 | 1650 | ±30 | 採用⚠ A=54.3% |
| NRSK-C5 | Phase 7 | 箸墓古墳周辺 | 木材 | ヤナギ属 自然木、樹幹（最外縁） | 中 | 1670 | ±30 | 採用 |
| NRSK-C6 | Phase 7 | 箸墓古墳周辺 | 木材 | ヤナギ属 自然木、枝（最外縁） | 低 | 1760 | ±30 | 採用 |
| NRSK-7 | Phase 7 | 箸墓古墳周辺 | 土器付着物 | 煤（甕胴外面） | 高 | 1790 | ±30 | 採用 |
| NRSK-C21 | Phase 7 | 箸墓古墳周辺 | 木材 | 乾燥木材（枝状、周濠下層） | 中 | 1820 | ±30 | 採用 |
| NRSK-C23 | Phase 7 | 箸墓古墳周辺 | 木材 | 材（周濠上層腐植土） | 中 | 1800 | ±30 | 採用 |
| NRSK-C24 | Phase 7 | 箸墓古墳周辺 | 木材 | 小枝（周濠上層腐植土） | 低 | 1720 | ±30 | 採用 |
| NRSK-C25 | Phase 7 | 箸墓古墳周辺 | 木材 | 枝（周濠上層腐植土） | 低 | 1800 | ±30 | 採用 |
| Phase 8 布留1式期（布留遺跡） |
| NRSK-46 | Phase 8 | 布留遺跡 | 土器付着物 | 煤（布留1式甕胴外面） | 高 | 1710 | ±20 | 採用 |
| NRSK-C49 | Phase 8 | 布留遺跡 | 種実 | ウリ 種子 | 低* | 1680 | ±20 | 採用⚠ A=41.7% |
| NRSK-C50_a | Phase 8 | 布留遺跡 | 種実 | モモ 核 | 低 | 1740 | ±20 | 採用 |
| NRSK-C50_b | Phase 8 | 布留遺跡 | 種実 | ウリ 種子 | 低 | 1730 | ±20 | 採用 |
| NRSK-C50_c | Phase 8 | 布留遺跡 | 種実 | ヘチマ 種子（推定） | 低 | 1710 | ±20 | 採用 |
| NRSK-55_a | Phase 8 | 布留遺跡 | 木材 | 木材（詳細は坂本2022付表参照） | 中 | 1750 | ±20 | 採用 |
| NRSK-55_b | Phase 8 | 布留遺跡 | 木材 | 木材（詳細は坂本2022付表参照） | 中 | 1800 | ±15 | 採用 |
| Phase 9 布留2式期（上之庄遺跡） |
| NRSK-59_a | Phase 9 | 上之庄遺跡 | 土器付着物 | 焦（甕胴下内面） | 高 | 1770 | ±20 | 採用 |
| NRSK-59_b | Phase 9 | 上之庄遺跡 | 土器付着物 | 煤（甕胴下外面） | 高 | 1775 | ±15 | 採用 |
| NRSK-C60 | Phase 9 | 上之庄遺跡 | 種実 | モモ 核 | 低* | 1710 | ±15 | 採用⚠ A=52.7% |
| ホケノ山古墳試料——最高信頼性・春成モデルが除外（本稿代替モデルで復権） |
小枝 No.1
（奥山2008） | ホケノ山 | ホケノ山古墳 | 木材 | 小枝（最外年輪含む） | 低 | 1710 | ±20 | 除外 |
小枝 No.2
（奥山2008） | ホケノ山 | ホケノ山古墳 | 木材 | 小枝（最外年輪含む） | 低 | 1690 | ±20 | 除外 |
【試料種別の凡例】煤＝土器外面付着炭化物（古木効果・使用木材の年代が反映される危険）、焦＝土器内面付着炭化物（リザーバー効果の危険）、吹＝土器外面付着物の特殊形態（原文記載の種別を採用）、カゴ材＝伝世・撹乱のリスク大。小枝・枝（最外縁）は最外年輪を含み古木効果なし。種実・炭化米は単年性（当年材）。いずれも理論的には最高信頼性。低*は材質は高信頼性だが当該モデルでA値が著しく低下しており層位撹乱の疑いがある試料。⚠は本稿追試でA値が基準（60%）を下回ることを確認した試料。試料種別の一次情報は春成ら（2011）付表、¹⁴C年代値はOxCalコード（付録E.2）記載値を使用。
E.2 検証条件と使用したOxCalコード
春成ら（2011）の論文に掲載された測定データ（異常値・混入とみなされたものを除く）を使用し、同論文が主張する土器型式および遺構の変遷順序（庄内0式→庄内1式→庄内3式→布留0式古相→布留0式新相→布留1式→布留2式）を、OxCalのSequenceおよびPhaseコマンドを用いて事前制約（Prior）として記述した。使用したコードの全文は以下の通りである。
Plot()
{
/* 較正曲線の指定：最新の国際標準 IntCal20 */
Curve("IntCal20","IntCal20.14c");
/* 春成モデルが想定する土器型式・遺構の変遷順序 */
Sequence("Harunari_Model_Sequence")
{
Boundary("Start");
Phase("1_Yamato_VI_4") /* 大和VI-4様式（大福28次） */
/* ※坂本（2022）は唐古・鍵遺跡（①〜③、NRTK系サンプル）を「参考資料」として */
/* 検討から除外しており、本モデルもこれに従いNRTK-32,-34,-35を含まない。 */
{
R_Date("NRSK-65", 1910, 20);
R_Date("NRSK-66", 1920, 15);
R_Date("NRSK-C63", 1925, 20);
};
Boundary("Trans_1_2");
Phase("2_Shonai_0") /* 庄内0式期 */
{
R_Date("NRSK-21", 1920, 25);
R_Date("NRSK-22_b", 1910, 25);
R_Date("NRSK-23_b", 1915, 25);
R_Date("NRSK-33_a", 1920, 25);
R_Date("NRSK-33_b", 1965, 25);
R_Date("NRSK-34_b", 1905, 25);
};
Boundary("Trans_2_3");
Phase("3_Shonai_1") /* 庄内1式期 */
{
R_Date("NRSK-24", 1930, 25);
R_Date("NRSK-25_a", 1925, 25);
R_Date("NRSK-25_b", 1910, 25);
R_Date("NRSK-26", 1910, 25);
};
Boundary("Trans_3_4");
Phase("4_Shonai_3") /* 庄内3式期 */
{
R_Date("NRSK-27", 1890, 25);
R_Date("NRSK-29_b", 1865, 25);
};
Boundary("Trans_4_5");
Phase("5_Makimuku_Ishizuka") /* ⑦ 纏向石塚墳丘墓期 */
/* ※坂本（2022）図7に基づく。NRSK-6（矢塚）は本文「除外」とあるが図7に[A:92]で */
/* 掲載されているため、図7の実装に準拠して採用。 */
{
R_Date("NRSK-C14", 1830, 30);
R_Date("NRSK-C16", 1910, 30);
R_Date("NRSK-C17", 1890, 30);
R_Date("NRSK-C18", 1880, 30);
R_Date("NRSK-6", 1820, 30);
};
Boundary("Trans_5_6");
/* ⑧ 坂本図7に従い築造前・築造中・箸墓築造前の3サブフェーズに分割 */
Sequence("Phase_8_Sub")
{
Phase("6a_Higashida_Before") /* ⑧-1 東田大塚墳丘墓築造前 */
{
R_Date("NRSK-2", 1860, 30);
R_Date("NRSK-3_a", 1820, 30);
R_Date("NRSK-3_b", 1780, 30);
R_Date("NRSK-C1", 1850, 30);
R_Date("NRSK-C2", 1730, 30);
};
Boundary("Trans_6a_6b");
Phase("6b_Higashida_During") /* ⑧-2 東田大塚築造中 */
{
R_Date("NRSK-C7", 1760, 30);
R_Date("NRSK-C8", 1730, 30);
};
Boundary("Trans_6b_6c");
Phase("6c_Hashihaka_Before") /* ⑧-3 箸墓古墳築造前 */
{
R_Date("NRKS-02", 1840, 30);
R_Date("NRKS-03", 1780, 30);
R_Date("NRKS-04", 1820, 30);
R_Date("NRKS-06", 1830, 40);
R_Date("NRKS-09", 1910, 40);
R_Date("NRKS-10", 1840, 40);
R_Date("NRKS-11", 1740, 40);
R_Date("NRKS-13", 1870, 40);
};
};
Boundary("Trans_6_7");
Phase("7_Furu_0_New") /* ⑨ 東田大塚壕埋没後・箸墓築造後 */
{
R_Date("NRSK-1", 1710, 30);
R_Date("NRSK-4_a", 1840, 30);
R_Date("NRSK-4_b", 1750, 30);
R_Date("NRSK-C4", 1650, 30);
R_Date("NRSK-C5", 1670, 30);
R_Date("NRSK-C6", 1760, 30);
R_Date("NRSK-7", 1790, 30);
R_Date("NRSK-C21", 1820, 30);
R_Date("NRSK-C23", 1800, 30);
R_Date("NRSK-C24", 1720, 30);
R_Date("NRSK-C25", 1800, 30);
};
Boundary("Trans_7_8");
Phase("8_Furu_1") /* ⑩ 布留1式期 */
{
R_Date("NRSK-46", 1710, 20);
R_Date("NRSK-C49", 1680, 20);
R_Date("NRSK-C50_a", 1740, 20);
R_Date("NRSK-C50_b", 1730, 20);
R_Date("NRSK-C50_c", 1710, 20);
R_Date("NRSK-55_a", 1750, 20);
R_Date("NRSK-55_b", 1800, 15);
};
Boundary("Trans_8_9");
Phase("9_Furu_2") /* ⑪ 布留2式期 */
{
R_Date("NRSK-59_a", 1770, 20);
R_Date("NRSK-59_b", 1775, 15);
R_Date("NRSK-C60", 1710, 15);
};
Boundary("End");
};
};
※本コードの実行結果（Amodel=19.3%）の詳細は、本文8.1.3節を参照のこと。
