Fitting a single dataset using Inverse Solver¶

In addition to using a parametric model or a control-points approach, we can also infer the production rates and its uncertainties using an inverse solver.

In [1]:

import numpy as np
import ticktack
from ticktack import fitting
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.style.use('seaborn-colorblind')

import numpy as np import ticktack from ticktack import fitting import matplotlib as mpl import matplotlib.pyplot as plt mpl.style.use('seaborn-colorblind')

In [2]:

sf = fitting.SingleFitter('Guttler15', 'Guttler15')
sf.load_data('miyake12.csv')
sf.compile_production_model("inverse_solver")

sf = fitting.SingleFitter('Guttler15', 'Guttler15') sf.load_data('miyake12.csv') sf.compile_production_model("inverse_solver")

Let's infer production rates from Monte Carlo samples of the d14c data in miyake12.csv using an inverse solver,

In [3]:

chain = sf.MC_reconstruct(iters=1000)
print("a sample from chain: \n", chain[0])
print("chain shape: ", chain.shape)

chain = sf.MC_reconstruct(iters=1000) print("a sample from chain: \n", chain[0]) print("chain shape: ", chain.shape)

100%|███████████████████████████████████████| 1000/1000 [00:11<00:00, 86.18it/s]

a sample from chain: 
 [ 3.30259166 -0.58035717  3.66990759 -1.22151059  3.84375996  2.44877767
  0.76190239  3.35334906  0.45990536 -1.69032714  5.3466338   1.61066658
 -0.98638563  4.77175444 -1.72507769 11.72212446 -0.24609793  1.83787816
 -1.52245896  4.8587668   0.79885272  3.56613786 -0.81719213  6.34347604
 -1.23213181  1.89105894  2.17688813  1.42199137]
chain shape:  (1000, 28)

Visualise the mean and the standard deviation of the chain,

In [4]:

fig, ax = plt.subplots(figsize=(6, 4), dpi=100)
mean = np.mean(chain, axis=0)
std = np.std(chain, axis=0)
ax.errorbar(sf.time_data, mean, color='k', drawstyle="steps")
ax.fill_between(sf.time_data, mean - std, mean + std,
                step='pre', alpha=0.6, facecolor=(0, 0, 0, .1),
                edgecolor=(0, 0, 0, 0.8), lw=1.5)
ax.set_xlabel("Year (CE)")
ax.set_ylabel("production rate (atoms/cm$^2$/s)");

fig, ax = plt.subplots(figsize=(6, 4), dpi=100) mean = np.mean(chain, axis=0) std = np.std(chain, axis=0) ax.errorbar(sf.time_data, mean, color='k', drawstyle="steps") ax.fill_between(sf.time_data, mean - std, mean + std, step='pre', alpha=0.6, facecolor=(0, 0, 0, .1), edgecolor=(0, 0, 0, 0.8), lw=1.5) ax.set_xlabel("Year (CE)") ax.set_ylabel("production rate (atoms/cm$^2$/s)");

Unlike our other approaches, we cannot directly call sf.dc14 on the inverse solver since this production rate model does not require burn-in. The reconstruction function below addresses this issue, but it is only appropriate for inverse solver,

In [5]:

def d14c(params):
    event = sf.run_event(y0=sf.steady_state_y0, params=(params,))
    binned_d14c = sf.cbm.bin_data(event[:, sf.box_idx], sf.oversample, sf.annual, growth=sf.growth)
    return (binned_d14c - 
            sf.steady_state_y0[sf.box_idx]) / sf.steady_state_y0[sf.box_idx] * 1000 + sf.d14c_data[0]

def d14c(params): event = sf.run_event(y0=sf.steady_state_y0, params=(params,)) binned_d14c = sf.cbm.bin_data(event[:, sf.box_idx], sf.oversample, sf.annual, growth=sf.growth) return (binned_d14c - sf.steady_state_y0[sf.box_idx]) / sf.steady_state_y0[sf.box_idx] * 1000 + sf.d14c_data[0]

Visualise the reconstruction,

In [6]:

fig, ax = plt.subplots(figsize=(6, 4), dpi=100)
colors = mpl.rcParams['axes.prop_cycle'].by_key()['color']
ax.errorbar(sf.time_data, sf.d14c_data, sf.d14c_data_error, fmt='o', 
            capsize=3, color="k", alpha=0.7, label="data")
ax.plot(sf.time_data, sf.d14c_data, color="k", alpha=0.7)
ax.errorbar(sf.time_data, d14c(mean), fmt='o', color=colors[0], label="reconstruction")
ax.legend(frameon=False)
ax.set_ylabel("$\Delta^{14}$C (‰)")
ax.set_xlabel("Year (CE)");

fig, ax = plt.subplots(figsize=(6, 4), dpi=100) colors = mpl.rcParams['axes.prop_cycle'].by_key()['color'] ax.errorbar(sf.time_data, sf.d14c_data, sf.d14c_data_error, fmt='o', capsize=3, color="k", alpha=0.7, label="data") ax.plot(sf.time_data, sf.d14c_data, color="k", alpha=0.7) ax.errorbar(sf.time_data, d14c(mean), fmt='o', color=colors[0], label="reconstruction") ax.legend(frameon=False) ax.set_ylabel("$\Delta^{14}$C (‰)") ax.set_xlabel("Year (CE)");