Skip to main content Skills Marketplace Discover and explore AI skills built by the community.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
Show prompt details A direct command skips the review prompt. Inspect the source before running it.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill radiocarbon-datingThe command stays on one line. Scroll horizontally to inspect it before copying.
Prefer a local copy? Download the files currently available to SkillsMP.
Download Zip Downloading... More from this repository
Related occupations SOC
Based on SOC occupation classification
name radiocarbon-dating description Radiocarbon (14C) dating calibration and Bayesian age modeling using IntCal20/SHCal20/Marine20 curves, sequence modeling with stratigraphic constraints, and uncertainty reporting.
tags ["archaeology","radiocarbon","bayesian","geochronology","calibration","dating"] version 1.0.0 authors [{"name":"awesome-rosetta-skills contributors","github":"@awesome-rosetta-skills"}] license MIT platforms ["claude-code","codex","gemini-cli","cursor"] dependencies ["numpy>=1.24.0","scipy>=1.10.0","matplotlib>=3.7.0","pandas>=2.0.0","requests>=2.28.0"] last_updated 2026-03-17
Radiocarbon Dating Analysis
Radiocarbon (14C) dating is the most widely used chronometric technique in archaeology,
providing calendar age estimates for organic materials up to ~50,000 years old. This skill
covers calibration against internationally agreed reference curves, Bayesian sequence
modeling, wiggle matching, and proper reporting conventions.
Background
Raw radiocarbon measurements are reported as "conventional radiocarbon ages" in years BP
(Before Present, where Present = 1950 CE). Because atmospheric 14C concentration has varied
over time, a calibration step is required to convert these measurements into calendar ages.
The IntCal20, SHCal20, and Marine20 calibration curves (Reimer et al. 2020) provide the
relationship between radiocarbon years and calendar years.
Key Concepts
BP (Before Present) : Radiocarbon years before 1950 CEcal BP : Calibrated (calendar) years before 1950 CEBCE/CE : Calendar years; BCE = cal BP + 1950 (for dates before 1950)Delta-14C (Δ14C) : Per-mil deviation of 14C/12C ratio from a standardF14C : Fraction modern — ratio of sample 14C to modern standard 14CReservoir effect : Marine/freshwater organisms appear older due to depleted 14C reservoirsMarine reservoir correction (ΔR) : Local deviation from global marine average
Installation
pip install numpy scipy matplotlib pandas requests
python -c "import urllib.request; urllib.request.urlretrieve(
'https://intcal.org/curves/intcal20.14c',
'intcal20.14c'
)"
Loading Calibration Curves
import numpy as np
import pandas as pd
import requests
from io import StringIO
from pathlib import Path
def load_calibration_curve (curve_name: str = "intcal20" ) -> pd.DataFrame:
"""
Load a calibration curve from file or download from intcal.org.
Parameters
----------
curve_name : str
One of 'intcal20', 'shcal20', 'marine20'
Returns
-------
pd.DataFrame
Columns: cal_BP, c14_age, c14_error, delta14C, delta14C_error
"""
curve_urls = {
"intcal20"
"https://intcal.org/curves/intcal20.14c"
"shcal20"
"https://intcal.org/curves/shcal20.14c"
"marine20"
"https://intcal.org/curves/marine20.14c"
"intcal20"
"intcal20.14c"
"shcal20"
"shcal20.14c"
"marine20"
"marine20.14c"
if
not
in
raise
f"Unknown curve: {curve_name} . Choose from {list (curve_urls)} "
if
else
print
f"Downloading {curve_name} from intcal.org ..."
30
for
in
if
not
"#"
and
"\n"
None
"cal_BP"
"c14_age"
"c14_error"
"delta14C"
"delta14C_error"
r"\s*,\s*"
"python"
"cal_BP"
True
return
def
_interpolate_curve
curve_df: pd.DataFrame, cal_years: np.ndarray
"""Interpolate curve c14_age and c14_error at given cal BP values."""
"cal_BP"
"c14_age"
"cal_BP"
"c14_error"
return
Single-Date Calibration import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from scipy.stats import norm
def calibrate_date (
c14_age: float ,
c14_error: float ,
curve_df: pd.DataFrame,
cal_range: tuple [float , float ] = (0 , 55000 float = 1.0 ,
) -> pd.DataFrame:
"""
Calibrate a single radiocarbon date against a calibration curve.
Parameters
----------
c14_age : float
Conventional radiocarbon age in BP.
c14_error : float
1-sigma laboratory measurement error.
curve_df : pd.DataFrame
Calibration curve loaded by load_calibration_curve().
cal_range : tuple
(min_cal_BP, max_cal_BP) search window.
resolution : float
Step size in calendar years (default 1 year).
Returns
-------
pd.DataFrame
Columns: cal_BP, probability (unnormalized then normalized), cal_CE.
"""
cal_years = np.arange(cal_range[0 ], cal_range[1 ] + resolution, resolution)
curve_mu, curve_sigma = _interpolate_curve(curve_df, cal_years)
combined_sigma = np.sqrt(c14_error**2 + curve_sigma**2 )
prob = norm.pdf(c14_age, loc=curve_mu, scale=combined_sigma)
prob /= prob.sum ()
cal_CE = 1950 - cal_years
return pd.DataFrame({
"cal_BP" : cal_years,
"cal_CE" : cal_CE,
"probability" : prob,
})
def hpd_intervals (cal_dist: pd.DataFrame, confidence: float = 0.95 ) -> list [tuple ]:
"""
Compute Highest Posterior Density (HPD) intervals.
Parameters
----------
cal_dist : pd.DataFrame
Output of calibrate_date().
confidence : float
Credible mass (e.g. 0.68 or 0.95).
Returns
-------
list of (start_cal_BP, end_cal_BP, probability_mass) tuples
"""
sorted_idx = np.argsort(cal_dist["probability" ].values)[::-1 ]
cum_prob = np.cumsum(cal_dist["probability" ].values[sorted_idx])
in_hpd = np.zeros(len (cal_dist), dtype=bool )
in_hpd[sorted_idx[cum_prob <= confidence]] = True
if cum_prob[np.sum (cum_prob <= confidence)] <= confidence + 1e-9 :
idx = np.sum (cum_prob <= confidence)
if idx < len (sorted_idx):
in_hpd[sorted_idx[idx]] = True
cal_years = cal_dist["cal_BP" ].values
intervals = []
in_block = False
start = None
block_prob = 0.0
for i, (yr, flag, p) in enumerate (
zip (cal_years, in_hpd, cal_dist["probability" ].values)
):
if flag and not in_block:
in_block = True
start = yr
block_prob = 0.0
if flag:
block_prob += p
if not flag and in_block:
intervals.append((start, cal_years[i - 1 ], block_prob))
in_block = False
if in_block:
intervals.append((start, cal_years[-1 ], block_prob))
return intervals
def summary_statistics (cal_dist: pd.DataFrame ) -> dict :
"""
Compute summary statistics for a calibrated date distribution.
Returns dict with keys: mean_cal_BP, median_cal_BP, mode_cal_BP,
sigma_cal_BP, hpd68, hpd95 (each as list of interval tuples).
"""
probs = cal_dist["probability" ].values
years = cal_dist["cal_BP" ].values
mean_bp = np.sum (years * probs)
variance = np.sum (probs * (years - mean_bp) ** 2 )
sigma = np.sqrt(variance)
cdf = np.cumsum(probs)
median_bp = years[np.searchsorted(cdf, 0.5 )]
mode_bp = years[np.argmax(probs)]
return {
"mean_cal_BP" : float (mean_bp),
"median_cal_BP" : float (median_bp),
"mode_cal_BP" : float (mode_bp),
"sigma_cal_BP" : float (sigma),
"mean_cal_CE" : float (1950 - mean_bp),
"hpd68" : hpd_intervals(cal_dist, 0.68 ),
"hpd95" : hpd_intervals(cal_dist, 0.95 ),
}
def plot_calibrated_date (
cal_dist: pd.DataFrame,
title: str = "Calibrated Radiocarbon Date" ,
c14_age: float = None ,
c14_error: float = None ,
curve_df: pd.DataFrame = None ,
show_hpd: bool = True ,
output_path: str = None ,
) -> None :
"""
Plot a calibrated radiocarbon date with optional calibration curve overlay.
Parameters
----------
cal_dist : pd.DataFrame
Output of calibrate_date().
title : str
Plot title.
c14_age, c14_error : float, optional
Original measurement for Gaussian overlay on 14C axis.
curve_df : pd.DataFrame, optional
Calibration curve to display.
show_hpd : bool
Shade 95% HPD intervals.
output_path : str, optional
Save figure to this path if provided.
"""
fig = plt.figure(figsize=(12 , 8 ))
gs = gridspec.GridSpec(2 , 2 , width_ratios=[3 , 1 ], height_ratios=[1 , 3 ],
hspace=0.05 , wspace=0.05 )
ax_curve = fig.add_subplot(gs[0 , 0 ])
ax_cal = fig.add_subplot(gs[1 , 0 ])
ax_c14 = fig.add_subplot(gs[1 , 1 ])
cal_years = cal_dist["cal_BP" ].values
probs = cal_dist["probability" ].values
if curve_df is not None :
mask = (curve_df["cal_BP" ] >= cal_years.min ()) & \
(curve_df["cal_BP" ] <= cal_years.max ())
sub = curve_df[mask]
ax_curve.plot(sub["cal_BP" ], sub["c14_age" ], "k-" , lw=1.5 , label="IntCal20" )
ax_curve.fill_between(
sub["cal_BP" ],
sub["c14_age" ] - 2 * sub["c14_error" ],
sub["c14_age" ] + 2 * sub["c14_error" ],
alpha=0.3 , color="gray" ,
)
if c14_age is not None :
ax_curve.axhline(c14_age, color="red" , lw=1 , ls="--" )
ax_curve.axhspan(c14_age - c14_error, c14_age + c14_error,
color="red" , alpha=0.15 )
ax_curve.set_ylabel("14C Age (BP)" )
ax_curve.xaxis.set_visible(False )
ax_curve.invert_xaxis()
ax_cal.fill_between(cal_years, probs, alpha=0.5 , color="steelblue" )
ax_cal.plot(cal_years, probs, color="steelblue" , lw=1 )
if show_hpd:
intervals_95 = hpd_intervals(cal_dist, 0.95 )
for start, end, mass in intervals_95:
mask = (cal_years >= start) & (cal_years <= end)
ax_cal.fill_between(cal_years[mask], probs[mask],
alpha=0.6 , color="steelblue" )
ax_cal.set_xlabel("Calibrated Age (cal BP)" )
ax_cal.set_ylabel("Probability Density" )
ax_cal.invert_xaxis()
if c14_age is not None and c14_error is not None :
y_range = np.linspace(c14_age - 4 * c14_error, c14_age + 4 * c14_error, 300 )
gauss = norm.pdf(y_range, c14_age, c14_error)
ax_c14.plot(gauss, y_range, color="red" , lw=1.5 )
ax_c14.fill_betweenx(y_range, gauss, alpha=0.3 , color="red" )
ax_c14.yaxis.set_visible(False )
ax_c14.set_xlabel("Probability" )
fig.suptitle(title, fontsize=14 , fontweight="bold" )
if output_path:
plt.savefig(output_path, dpi=150 , bbox_inches="tight" )
plt.tight_layout()
plt.show()
Bayesian Sequence Modeling from scipy.optimize import minimize
from typing import Optional
def sequence_model (
dates: list [dict ],
constraints: list [tuple [int , int ]] = None ,
n_iterations: int = 50000 ,
burn_in: int = 10000 ,
seed: int = 42 ,
) -> list [dict ]:
"""
Bayesian sequence model using MCMC for stratigraphically ordered dates.
Each date dict must have keys:
label: str
c14_age: float (BP)
c14_error: float
curve_df: pd.DataFrame (calibration curve)
constraints: list of (i, j) meaning date[i] is OLDER than date[j]
(i.e., date[i].cal_BP > date[j].cal_BP)
Returns list of dicts with posterior samples and statistics.
Parameters
----------
dates : list of dict
Each entry: {label, c14_age, c14_error, curve_df}
constraints : list of (int, int)
Ordering constraints: (older_index, younger_index)
n_iterations : int
Total MCMC iterations.
burn_in : int
Discarded warm-up iterations.
seed : int
Random seed for reproducibility.
"""
rng = np.random.default_rng(seed)
n = len (dates)
theta = np.zeros(n)
cal_dists = []
for i, d in enumerate (dates):
cd = calibrate_date(d["c14_age" ], d["c14_error" ], d["curve_df" ])
cal_dists.append(cd)
theta[i] = float (summary_statistics(cd)["median_cal_BP" ])
if constraints is None :
constraints = []
def log_likelihood (th ):
ll = 0.0
for i, d in enumerate (dates):
cd = cal_dists[i]
prob = np.interp(th[i], cd["cal_BP" ].values, cd["probability" ].values,
left=0.0 , right=0.0 )
if prob <= 0 :
return -np.inf
ll += np.log(prob)
return ll
def satisfies_constraints (th ):
for (older, younger) in constraints:
if th[older] <= th[younger]:
return False
return True
proposal_sigma = 50.0
current_ll = log_likelihood(theta)
samples = np.zeros((n_iterations - burn_in, n))
n_accepted = 0
for it in range (n_iterations):
idx = rng.integers(0 , n)
theta_prop = theta.copy()
theta_prop[idx] += rng.normal(0 , proposal_sigma)
if not satisfies_constraints(theta_prop):
continue
prop_ll = log_likelihood(theta_prop)
log_alpha = prop_ll - current_ll
if np.log(rng.uniform()) < log_alpha:
theta = theta_prop
current_ll = prop_ll
n_accepted += 1
if it >= burn_in:
samples[it - burn_in] = theta
acceptance_rate = n_accepted / n_iterations
print (f"MCMC acceptance rate: {acceptance_rate:.3 f} " )
results = []
for i, d in enumerate (dates):
samp = samples[:, i]
results.append({
"label" : d["label" ],
"c14_age" : d["c14_age" ],
"c14_error" : d["c14_error" ],
"posterior_samples" : samp,
"mean_cal_BP" : float (np.mean(samp)),
"median_cal_BP" : float (np.median(samp)),
"std_cal_BP" : float (np.std(samp)),
"hpd95_lower" : float (np.percentile(samp, 2.5 )),
"hpd95_upper" : float (np.percentile(samp, 97.5 )),
"mean_cal_CE" : float (1950 - np.mean(samp)),
})
return results
def wiggle_match (
radiocarbon_ages: list [float ],
errors: list [float ],
time_gaps: list [float ],
curve_df: pd.DataFrame,
search_range: tuple = (0 , 10000 ) -> dict :
"""
Wiggle matching: fit a floating tree-ring sequence to the calibration curve.
Parameters
----------
radiocarbon_ages : list of float
14C ages (BP) for each ring, outermost first.
errors : list of float
Measurement errors for each ring.
time_gaps : list of float
Known gaps in years between rings (len = len(radiocarbon_ages) - 1).
curve_df : pd.DataFrame
Calibration curve.
search_range : tuple
(min_cal_BP, max_cal_BP) for the outermost ring.
Returns
-------
dict with best_fit_cal_BP (outermost ring), chi2, and probability array.
"""
n_rings = len (radiocarbon_ages)
candidate_years = np.arange(search_range[0 ], search_range[1 ])
chi2_values = np.zeros(len (candidate_years))
for idx, outer_cal_BP in enumerate (candidate_years):
chi2 = 0.0
for ring_i in range (n_rings):
gap = sum (time_gaps[:ring_i]) if ring_i > 0 else 0
cal_bp_ring = outer_cal_BP + gap
curve_mu, curve_sigma = _interpolate_curve(curve_df, np.array([cal_bp_ring]))
combined_var = errors[ring_i]**2 + curve_sigma[0 ]**2
chi2 += (radiocarbon_ages[ring_i] - curve_mu[0 ])**2 / combined_var
chi2_values[idx] = chi2
prob = np.exp(-0.5 * chi2_values)
prob /= prob.sum ()
best_idx = np.argmin(chi2_values)
return {
"best_fit_cal_BP" : float (candidate_years[best_idx]),
"best_fit_cal_CE" : float (1950 - candidate_years[best_idx]),
"chi2_min" : float (chi2_values[best_idx]),
"candidate_years" : candidate_years,
"probability" : prob,
}
Example 1: Calibrate a Single Radiocarbon Date with IntCal20
def example_single_date ():
print ("=== Example 1: Single Date Calibration ===\n" )
curve = load_calibration_curve("intcal20" )
print (f"Loaded IntCal20: {len (curve)} data points" )
c14_age = 2850
c14_error = 35
cal_dist = calibrate_date(c14_age, c14_error, curve_df=curve)
stats = summary_statistics(cal_dist)
print (f"\nRadiocarbon age: {c14_age} ± {c14_error} BP" )
print (f"Calibrated mean: {stats['mean_cal_BP' ]:.0 f} cal BP "
f"({stats['mean_cal_CE' ]:.0 f} CE/BCE)" )
print (f"Calibrated median: {stats['median_cal_BP' ]:.0 f} cal BP" )
print (f"Sigma (1σ): ± {stats['sigma_cal_BP' ]:.0 f} years" )
print ("\n68.3% HPD intervals (cal BP):" )
for start, end, mass in stats["hpd68" ]:
ce_start = 1950 - end
ce_end = 1950 - start
era_s = "BCE" if ce_start < 0 else "CE"
era_e = "BCE" if ce_end < 0 else "CE"
print (f" {start:.0 f} – {end:.0 f} cal BP "
f"({abs (ce_start):.0 f} {era_s} – {abs (ce_end):.0 f} {era_e} ) "
f"[{mass*100 :.1 f} %]" )
print ("\n95.4% HPD intervals (cal BP):" )
for start, end, mass in stats["hpd95" ]:
ce_start = 1950 - end
ce_end = 1950 - start
era_s = "BCE" if ce_start < 0 else "CE"
era_e = "BCE" if ce_end < 0 else "CE"
print (f" {start:.0 f} – {end:.0 f} cal BP "
f"({abs (ce_start):.0 f} {era_s} – {abs (ce_end):.0 f} {era_e} ) "
f"[{mass*100 :.1 f} %]" )
plot_calibrated_date(
cal_dist,
title=f"Charcoal from Hearth — {c14_age} ± {c14_error} BP" ,
c14_age=c14_age,
c14_error=c14_error,
curve_df=curve,
output_path="calibrated_hearth.png" ,
)
print ("\nPlot saved to calibrated_hearth.png" )
if __name__ == "__main__" :
example_single_date()
Example 2: Bayesian Sequence Model for a Burial Site
def example_burial_sequence ():
print ("=== Example 2: Bayesian Sequence Model — Bronze Age Burial Site ===\n" )
curve = load_calibration_curve("intcal20" )
burial_dates = [
{"label" : "Layer 1 (youngest)" , "c14_age" : 2980 , "c14_error" : 40 , "curve_df" : curve},
{"label" : "Layer 2" , "c14_age" : 3060 , "c14_error" : 35 , "curve_df" : curve},
{"label" : "Layer 3" , "c14_age" : 3150 , "c14_error" : 45 , "curve_df" : curve},
{"label" : "Layer 4" , "c14_age" : 3280 , "c14_error" : 50 , "curve_df" : curve},
{"label" : "Layer 5 (oldest)" , "c14_age" : 3390 , "c14_error" : 40 , "curve_df" : curve},
]
constraints = [
(1 , 0 ),
(2 , 1 ),
(3 , 2 ),
(4 , 3 ),
]
print ("Running MCMC sequence model (this may take 10–30 seconds)..." )
results = sequence_model(
burial_dates,
constraints=constraints,
n_iterations=30000 ,
burn_in=5000 ,
seed=2024 ,
)
print ("\nPosterior estimates (sequence-constrained):\n" )
print (f"{'Label' :<25 } {'Mean (cal BP)' :>14 } {'95% HPD (cal BP)' :>22 } {'Mean (CE/BCE)' :>14 } " )
print ("-" * 80 )
for r in results:
ce = r["mean_cal_CE" ]
era = "BCE" if ce < 0 else "CE"
print (
f"{r['label' ]:<25 } "
f"{r['mean_cal_BP' ]:>14.0 f} "
f"{r['hpd95_lower' ]:>10.0 f} – {r['hpd95_upper' ]:<10.0 f} "
f"{abs (ce):>10.0 f} {era} "
)
print ("\nSequence consistency check:" )
for i in range (len (results) - 1 ):
older = results[i + 1 ]["mean_cal_BP" ]
younger = results[i]["mean_cal_BP" ]
ok = "OK" if older > younger else "VIOLATION"
print (f" {results[i+1 ]['label' ]} > {results[i]['label' ]} : "
f"{older:.0 f} > {younger:.0 f} [{ok} ]" )
fig, axes = plt.subplots(len (results), 1 , figsize=(12 , 10 ), sharex=True )
colors = plt.cm.viridis(np.linspace(0.2 , 0.8 , len (results)))
for i, (r, ax) in enumerate (zip (results, axes)):
samp = r["posterior_samples" ]
ax.hist(samp, bins=80 , density=True , color=colors[i], alpha=0.7 , edgecolor="none" )
ax.axvline(r["mean_cal_BP" ], color="red" , lw=1.5 , ls="--" , label="Mean" )
ax.axvspan(r["hpd95_lower" ], r["hpd95_upper" ], alpha=0.15 , color="red" )
ax.set_ylabel(r["label" ][:20 ], fontsize=8 )
ax.yaxis.set_ticklabels([])
axes[-1 ].set_xlabel("Calibrated Age (cal BP)" )
axes[-1 ].invert_xaxis()
fig.suptitle("Bayesian Sequence Model — Bronze Age Burial Site" , fontweight="bold" )
plt.tight_layout()
plt.savefig("burial_sequence.png" , dpi=150 , bbox_inches="tight" )
print ("\nSequence plot saved to burial_sequence.png" )
if __name__ == "__main__" :
example_burial_sequence()
Reporting Conventions When reporting calibrated dates in publications, follow OxCal/IntCal conventions:
Always specify the calibration curve and software version
Report both 68.3% (1σ) and 95.4% (2σ) HPD intervals
State the original 14C measurement in BP with ± 1σ laboratory error
For example: "2850 ± 35 BP (OxCal 4.4, IntCal20); 2σ: 3090–2870 cal BP (1060–920 BCE)"
Delta-14C and F14C Conversions def c14_age_to_F14C (c14_age_BP: float ) -> float :
"""Convert conventional radiocarbon age (BP) to Fraction Modern (F14C)."""
return np.exp(-c14_age_BP / 8267.0 )
def F14C_to_c14_age (f14c: float ) -> float :
"""Convert Fraction Modern back to conventional radiocarbon age (BP)."""
if f14c <= 0 :
raise ValueError("F14C must be > 0" )
return -8267.0 * np.log(f14c)
def delta14C_to_F14C (delta14C_permil: float ) -> float :
"""Convert Δ14C (per mil) to F14C."""
return 1.0 + delta14C_permil / 1000.0
if __name__ == "__main__" :
age = 3000
f = c14_age_to_F14C(age)
recovered = F14C_to_c14_age(f)
print (f"3000 BP → F14C={f:.4 f} → {recovered:.1 f} BP (round-trip check)" )
References
Reimer, P.J., et al. (2020). The IntCal20 Northern Hemisphere Radiocarbon Age Calibration
Curve (0–55 cal kBP). Radiocarbon , 62(4), 725–757.
Hogg, A.G., et al. (2020). SHCal20 Southern Hemisphere Calibration, 0–55,000 Years cal BP.
Radiocarbon , 62(4), 759–778.
Bronk Ramsey, C. (2009). Bayesian Analysis of Radiocarbon Dates. Radiocarbon , 51(1), 337–360.
Buck, C.E., Kenworthy, J.B., Litton, C.D., & Smith, A.F.M. (1991). Combining archaeological
and radiocarbon information: a Bayesian approach to calibration. Antiquity , 65, 808–821.
