Romont site, MLE with 6 free parameters (grid search method).
For more info about the method used, see the notebook Inference_Notes.
This notebook has the following external dependencies:
In [1]:
import math
import csv
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
import yaml
%matplotlib inline
The mathematical model is available in the models Python module (see the notebook Models)
In [2]:
import models
In [3]:
profile_data = pd.read_csv('profiles_data/romont_10Be_26Al_profile_data.csv',
index_col='sample',
delim_whitespace=True,
quoting=csv.QUOTE_NONNUMERIC, quotechar='\"',
dtype={'depth': 'f', 'depth_g-cm-2': 'f',
'C': 'f', 'std': 'f'})
profile_data
Out[3]:
In [4]:
with open('profiles_data/romont_10Be_26Al_settings.yaml') as f:
romont_settings = yaml.load(f)
romont_settings
Out[4]:
The dataset is stored as a :class:pandas.DataFrame object.
The grid search method is implemented in the gridsearch module (see the notebook Grid-Search for more info).
In [5]:
import gridsearch
Create a new object for setup and results
In [6]:
gstest = gridsearch.CosmogenicInferenceGC(description='Romont 6 parameters')
Set the data
In [7]:
gstest.set_profile_measured(
profile_data['depth'].values,
profile_data['C'].values,
profile_data['std'].values,
profile_data['nuclide'].values
)
Set the model. The fit is based on both 10Be and 26Al concentrations vs. depth. The model also takes into account an unknown period exposure_burial, which started when the measured depth profile has been buried by other deposits (instantaneous burial is assumed, and we also assume that no erosion took place since the burial). The total age of the deposits is therefore given by exposure + exposure_burial.
In [8]:
P0_26Al_10Be_ratio = 6.75 # ratio of 26Al / 10Be production rates
P_0_10Be = romont_settings['P_0']
P_0_26Al = P0_26Al_10Be_ratio * P_0_10Be
burial_depth = 300 # burial depth [cm]
# build slices to distinguish between 10Be and 26Al
i10Be = gstest.oprofile['nuclide'] == "10Be"
i26Al = gstest.oprofile['nuclide'] == "26Al"
s10Be = np.s_[i10Be,...]
s26Al = np.s_[i26Al,...]
def C_10Be_26Al_romont(depth, erosion,
exposure, exposure_burial,
density, inheritance_10Be,
inheritance_26Al):
"""
10Be and 26Al for Romont with burial.
"""
C = np.empty_like(depth)
Cb = np.empty_like(depth)
# "dummy" broadcast to create full arrays
C = depth * erosion * exposure * exposure_burial \
* density * inheritance_10Be * inheritance_26Al \
* 0.
Cb = np.zeros_like(C)
C[s10Be] += models.C_10Be(depth[s10Be] - burial_depth,
erosion, exposure,
density, inheritance_10Be,
P_0=P_0_10Be)
C[s26Al] += models.C_26Al(depth[s26Al] - burial_depth,
erosion, exposure,
density, inheritance_26Al,
P_0=P_0_26Al)
Cb[s10Be] += models.C_10Be(depth[s10Be],
0., exposure_burial,
density, C[s10Be],
P_0=P_0_10Be)
Cb[s26Al] += models.C_26Al(depth[s26Al],
0., exposure_burial,
density, C[s26Al],
P_0=P_0_26Al)
return Cb
gstest.set_profile_model(C_10Be_26Al_romont)
Define the parameters to fit and their search ranges / steps. The order must be the same than the order of the arguments of the function used for the model!
In [9]:
gstest.set_parameter(
'erosion_rate',
[0., 0.5e-3, 20j],
stats.uniform(loc=0, scale=0.5e-3).pdf
)
gstest.set_parameter(
'exposure_time',
[1e5, 2.5e5, 20j],
stats.uniform(loc=1e5, scale=2.5e5).pdf
)
gstest.set_parameter(
'exposure_time_burial',
[0., 7e5, 20j],
stats.uniform(loc=0., scale=7e5).pdf
)
gstest.set_parameter(
'soil_density',
[2.0, 6.0, 20j],
stats.uniform(loc=2.0, scale=6.0).pdf
)
gstest.set_parameter(
'inheritance_10Be',
[2e5, 2.6e5, 20j],
stats.uniform(loc=2e5, scale=2.6e5).pdf
)
gstest.set_parameter(
'inheritance_26Al',
[0.6e6, 1.5e6, 20j],
stats.uniform(loc=0.6e6, scale=1e6).pdf
)
Grid search setup summary
In [10]:
print gstest.setup_summary()
Perform Maximum likelihood estimation on the search grid
In [11]:
gstest.compute_mle()
Get the MLE (i.e., the parameter values at the maximum likelihood), in the same order than the definition of the parameters
In [12]:
gstest.mle
Out[12]:
Plot the profile log-likelihood for each parameter. The blue lines represent the difference between the profile log-likelihood and the maximum log-likelihood, The intersections between the blue line and the black lines define the confidence intervals at the given confidence levels (based on the likelihood ratio test). The red lines indicate the true values.
In [13]:
%matplotlib inline
def plot_proflike1d(cobj, pname, clevels=[0.68, 0.95, 0.997],
true_val=None, ax=None):
p = cobj.parameters[pname]
pindex = cobj.parameters.keys().index(pname)
x = cobj.grid[pindex].flatten()
proflike = cobj.proflike1d[pindex]
if ax is None:
ax = plt.subplot(111)
difflike = proflike - cobj.maxlike
ax.plot(x, difflike, label='profile loglike')
ccrit = gridsearch.profile_likelihood_crit(
cobj.proflike1d[pindex],
cobj.maxlike,
clevels=clevels
)
ccrit -= cobj.maxlike
for lev, cc in zip(clevels, ccrit):
l = ax.axhline(cc, color='k')
hpos = x.min() + (x.max() + x.min()) * 0.05
ax.text(hpos, cc, str(lev * 100),
size=9, color = l.get_color(),
ha="center", va="center",
bbox=dict(ec='1',fc='1'))
if true_val is not None:
ax.axvline(true_val, color='r')
plt.setp(ax, xlabel=pname,
ylabel='profile log-like - max log-like',
xlim=p['range'][0:2],
ylim=[ccrit[-1], 0.])
def plot_proflike1d_all(cobj, n_subplot_cols=2, **kwargs):
n_subplots = len(cobj.parameters)
n_subplot_rows = int(math.ceil(1. *
n_subplots /
n_subplot_cols))
fig, aax = plt.subplots(nrows=n_subplot_rows,
ncols=n_subplot_cols,
**kwargs)
axes = aax.flatten()
fig.text(0.5, 0.975,
"Profile log-like: " + cobj.description,
horizontalalignment='center',
verticalalignment='top')
for i, pname in enumerate(cobj.parameters.keys()):
ax = axes[i]
plot_proflike1d(cobj, pname,
ax=ax)
plt.tight_layout()
plt.subplots_adjust(top=0.93)
plot_proflike1d_all(gstest, figsize=(11, 10))
Show the profile log-likelihood for couples of parameters. Confidence regions are also shown (also based on the likelihood ratio test).
In [14]:
def plot_proflike2d(cobj, p1p2, ax=None,
cmap='Blues', show_colorbar=True):
pname1, pname2 = p1p2
idim = cobj.parameters.keys().index(pname2)
jdim = cobj.parameters.keys().index(pname1)
if ax is None:
ax = plt.subplot(111)
X, Y = np.meshgrid(cobj.grid[idim].flatten(),
cobj.grid[jdim].flatten())
difflike = cobj.proflike2d[idim][jdim] - cobj.maxlike
ccrit = gridsearch.profile_likelihood_crit(
cobj.proflike2d[idim][jdim],
cobj.maxlike,
clevels=[0.68, 0.95]
)
ccrit -= cobj.maxlike
contours = np.linspace(np.median(difflike),
0,
10)
P2D = ax.contourf(Y, X, difflike,
contours,
cmap=plt.get_cmap(cmap))
ci68 = ax.contour(Y, X, difflike,
[ccrit[0]], colors='w',
linestyles='solid')
plt.clabel(ci68, fontsize=8, inline=True,
fmt='68')
ci95 = ax.contour(Y, X, difflike,
[ccrit[1]], colors=['k'],
linestyles='solid')
plt.clabel(ci95, fontsize=8, inline=True,
fmt='95')
ax.scatter(cobj.mle[jdim], cobj.mle[idim],
marker='*', s=60, c='w')
plt.setp(ax, xlabel=pname1, ylabel=pname2)
if show_colorbar:
plt.colorbar(P2D, ax=ax)
#ax.axhline(true_exposure, color='r')
#ax.axvline(true_erosion, color='r')
fig = plt.figure(figsize=(11, 11))
ax = plt.subplot(421)
plot_proflike2d(gstest, ('erosion_rate', 'exposure_time'), ax=ax)
ax2 = plt.subplot(422)
plot_proflike2d(gstest, ('erosion_rate', 'exposure_time_burial'), ax=ax2)
ax3 = plt.subplot(423)
plot_proflike2d(gstest, ('exposure_time', 'inheritance_10Be'), ax=ax3)
ax4 = plt.subplot(424)
plot_proflike2d(gstest, ('exposure_time', 'inheritance_26Al'), ax=ax4)
ax5 = plt.subplot(425)
plot_proflike2d(gstest, ('exposure_time', 'soil_density'), ax=ax5)
ax6 = plt.subplot(426)
plot_proflike2d(gstest, ('soil_density', 'inheritance_26Al'), ax=ax6)
ax7 = plt.subplot(427)
plot_proflike2d(gstest, ('exposure_time_burial', 'inheritance_26Al'), ax=ax7)
ax8 = plt.subplot(428)
plot_proflike2d(gstest, ('exposure_time_burial', 'inheritance_10Be'), ax=ax8)
plt.tight_layout()
Plot the measured concentrations and the predicted profile corresponding to the best fitted data model
In [15]:
sns.set_context('notebook')
depths = np.linspace(profile_data['depth'].min(),
profile_data['depth'].max(),
100)
print
Cm_fitted = C_10Be_26Al_romont(gstest.oprofile['depth'],
*gstest.mle)
plt.figure()
plt.plot(Cm_fitted[s10Be],
-gstest.oprofile['depth'][s10Be],
label='10Be best-fitted model')
plt.plot(Cm_fitted[s26Al],
-gstest.oprofile['depth'][s26Al],
label='26Al best-fitted model')
plt.errorbar(gstest.oprofile['C'][s10Be],
-gstest.oprofile['depth'][s10Be],
xerr=gstest.oprofile['std'][s10Be],
fmt='o', markersize=4,
label='10Be data')
plt.errorbar(gstest.oprofile['C'][s26Al],
-gstest.oprofile['depth'][s26Al],
xerr=gstest.oprofile['std'][s26Al],
fmt='o', markersize=4,
label='26Al data')
plt.setp(plt.gca(),
xlabel='10Be or 26Al concentration [atoms g-1]',
ylabel='-1 * depth [cm]',
xlim=[0, None], ylim=[None, 0])
plt.legend(loc='lower right')
Out[15]:
In [15]: