<jupyter_start><jupyter_text>Bayesian approach - Test case - 4 free parameters  An example of applying the Bayesian approach with 4 free parameters, using the [PyMC](https://github.com/pymc-devs/pymc) package.For more info about the method used, see the notebook [Inference_Notes](Inference_Notes.ipynb).This example (a test case) is based on a generic dataset of 10Be concentration vs. depth, which is drawn from a distribution with given "true" parameters.This notebook has the following external dependencies:<jupyter_code>import math

import numpy as np
import pandas as pd
import pymc
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline<jupyter_output><empty_output><jupyter_text>The mathematical (deterministic) model The mathematical model is available in the `models` Python module (see the notebook [Models](Models.ipynb)).<jupyter_code>import models<jupyter_output><empty_output><jupyter_text>The data The dataset is generated using the following parameter values<jupyter_code># the true parameters 
true_erosion = 1e-3
true_exposure = 1e5
true_density = 2.0
true_inheritance = 0.

# depths and sample size
depth_minmax = [50, 500]
N = 30

# perturbations
err_magnitude = 20.
err_variability = 5.<jupyter_output><empty_output><jupyter_text>The `gendata` Python module is used to generate the dataset (see the notebook [Datasets](Datasets.ipynb)).<jupyter_code>import gendata

profile_data = gendata.generate_dataset(
    models.C_10Be,
    [true_erosion, true_exposure, true_density, true_inheritance],
    zlimits=depth_minmax,
    n=N,
    err=[err_magnitude, err_variability]
)

profile_data<jupyter_output><empty_output><jupyter_text>The statistical model (used for computing the PPD values) PyMC allows to define the statistical model in a very nice way. The "nodes" of the model - i.e., the elements that contribute to the Posterior probability distribution - must be defined first.- Define the prior probability distribution for each free parameter (here the uniform distribution is used, with given bounds)<jupyter_code>erosion_rate = pymc.Uniform('erosion_rate', 0., 4e-3)
exposure_time = pymc.Uniform('exposure_time', 0., 4e5)
soil_density = pymc.Uniform('soil_density', 1.8, 2.3)
inheritance = pymc.Uniform('inheritance', 0, 3e5)<jupyter_output><empty_output><jupyter_text>- Define the deterministic part, i.e., the model used to predict nucleide concentration vs. depth<jupyter_code>@pymc.deterministic
def mprofile(erosion=erosion_rate, exposure=exposure_time,
             density=soil_density, inheritance=inheritance):
    """
    Modelled concentration profile.
    """
    return models.C_10Be(profile_data['depth'].values,
                         erosion, exposure,
                         density, inheritance)<jupyter_output><empty_output><jupyter_text>- Define the distribution of measured nucleide concentrations, which is here assumed to be a Gaussian with $\mu$ is predicted by the model above and $\sigma$ the known measurement errors. The likelihood will be derived from this distribution.<jupyter_code>d = pymc.Normal('data_distribution', mprofile, 1./profile_data['std']**2,
                value=profile_data['C'].values, observed=True)<jupyter_output><empty_output><jupyter_text>PyMC builds the model automatically from all nodes.<jupyter_code>smodel = pymc.Model({'erosion_rate': erosion_rate,
                     'exposure_time': exposure_time,
                     'soil_density': soil_density,
                     'inheritance': inheritance,
                     'mprofile': mprofile,
                     'data_distribution': d})<jupyter_output><empty_output><jupyter_text>PyMC use the prior distributions to generate a random initial value for each parameter (further needed)<jupyter_code>print erosion_rate.value
print exposure_time.value
print soil_density.value
print inheritance.value<jupyter_output>0.0030383483373
353740.615011
1.92284833369
58714.6972633<jupyter_text>The log-PPD value that correspond to the initial values is given by the `logp` attribute<jupyter_code>smodel.logp<jupyter_output><empty_output><jupyter_text>Fitting the statistical model :class:`pymc.MAP` allows to "clone" a model and set all parameters to their maximum a posteriori values. It uses the scipy's [optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html) package to maximize the log-PPD over the parameter space.<jupyter_code>maxppd_smodel = pymc.MAP(smodel)
maxppd_smodel.fit(method='fmin', iterlim=15000, tol=1e-6, verbose=0)<jupyter_output><empty_output><jupyter_text>See the resulting optimized values assigned to each parameter (and compare it to the true values defined above, it should be closer than the initial values)<jupyter_code>print erosion_rate.value
print exposure_time.value
print soil_density.value
print inheritance.value<jupyter_output>0.00399999999843
398135.161274
2.2999872895
0.00526119438799<jupyter_text>The log-PPD value at the maximum<jupyter_code>maxppd_smodel.logp_at_max<jupyter_output><empty_output><jupyter_text>Markov-Chain Monte-Carlo sampling :class:`pymc.MCMC` allows to sample the parameter space following the PPD with MCMC. It is a good idea to use the results of the log-PPD maximization as initial values for MCMC. Metropolis-Hastings is the default samplerused by MCMC.<jupyter_code>mcmc_smodel = pymc.MCMC(smodel)<jupyter_output><empty_output><jupyter_text>Generate samples, discard the first $n$ samples (burn-in, `burn` argument) and keep only one sample every $t$ samples (thinning, `thin` argument) in order to obtain *idependent* samples *from the PPD*.<jupyter_code>mcmc_smodel.sample(50000, burn=5000, thin=40)<jupyter_output>[-----------------100%-----------------] 50000 of 50000 complete in 84.6 sec<jupyter_text>Plot some characteristics of the samples for each parameter in order to evaluate the quality of sampling.- trace: plot of the values of the - burnt and thinned - samples (in the order of which it has been generated), useful to check if the burn-in setting was chosen appropriately (no trend visible for the firsts samples)- acorr: autocorrelation of the values of the - burnt and thinned - samples, useful to check if the thinning setting was chosen appropriately (rapid decrease of autocorrelation)- histogram of the samples<jupyter_code>pymc.Matplot.plot(mcmc_smodel)<jupyter_output>Plotting exposure_time
Plotting soil_density
Plotting erosion_rate
Plotting inheritance<jupyter_text>Print a summary about the distribution of the samples for each parameter. It shows the values of various moments and percentiles of the PPD, the 95% Highest Posterior Density (HPD) interval (see [here](http://www.bayesian-inference.com/credible) for an explanation) and an estimation of the Monte-Carlo integration error.<jupyter_code>mcmc_smodel.summary()<jupyter_output>exposure_time:
 
	Mean             SD               MC Error        95% HPD interval
	------------------------------------------------------------------
	78405.582        12864.894        1195.936[  57536.206  102054.876]
	
	
	Posterior quantiles:
	
	2.5             25              50              75             97.5
	 |---------------|===============|===============|---------------|
	59099.561        67750.8         76585.822      88069.717     105402.98
	

soil_density:
 
	Mean             SD               MC Error        95% HPD interval
	------------------------------------------------------------------
	1.983            0.07             0.005            [ 1.832  2.108]
	
	
	Posterior quantiles:
	
	2.5             25              50              75             97.5
	 |---------------|===============|===============|---------------|
	1.852            1.933           1.983          2.029         2.129
	

erosion_rate:
 
	Mean             SD               MC Error        95% HPD inte[...]<jupyter_text>Show the marginal distributions of the generated samples for each parameter. Continuous lines are obtained from the samples by Kernel Density Estimation (KDE). The result is an estimation of the marginal PPD of each parameter. The red lines indicate the true values.<jupyter_code>sns.set(style="darkgrid")
sns.set_context('paper')

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(8, 6))

sns.distplot(mcmc_smodel.trace('exposure_time')[:], 
             hist=False, kde_kws={"shade": True},
             ax=axes[0][0])
axes[0][0].axvline(true_exposure, color='r')
axes[0][0].set_xlabel('exposure time [yr]')

sns.distplot(mcmc_smodel.trace('erosion_rate')[:] * 1e4, 
             hist=False, kde_kws={"shade": True},
             ax=axes[0][1])
axes[0][1].axvline(true_erosion * 1e4, color='r')
axes[0][1].set_xlabel('erosion rate [m Myr-1]')

sns.distplot(mcmc_smodel.trace('soil_density')[:], 
             hist=False, kde_kws={"shade": True},
             ax=axes[1][0])
axes[1][0].axvline(true_density, color='r')
axes[1][0].set_xlabel('soil_density [g cm-3]')

sns.distplot(mcmc_smodel.trace('inheritance')[:], 
             hist=False, kde_kws={"shade": True},
             ax=axes[1][1])
axes[1][1].axvline(true_inheritance, color='r')
axes[1][1].set_xlabel('inheritance [atoms g-1]')

plt.setp(axes, yticks=[])<jupyter_output><empty_output><jupyter_text>Join plot of the generated samples with marginals, obtained by KDE. This represents an approximation of the 2D marginal PPDs. The red lines indicate the true values.<jupyter_code>sns.set(style="white")
sns.set_context('paper')

# exposure vs. erosion
p = sns.jointplot(mcmc_smodel.trace('erosion_rate')[:] * 1e4,
                  mcmc_smodel.trace('exposure_time')[:],
                  xlim=[0, 15], ylim=[4e4, 1.2e5],
                  kind="kde", size=5, space=0)
plt.setp(p.ax_joint,
         xlabel='erosion rate [m Myr-1]',
         ylabel='exposure time [yr]')
p.ax_joint.axhline(true_exposure, color='r')
p.ax_joint.axvline(true_erosion * 1e4, color='r')

# exposure vs. inheritance
p = sns.jointplot(mcmc_smodel.trace('inheritance')[:],
                  mcmc_smodel.trace('exposure_time')[:],
                  xlim=[0, 4000], ylim=[4e4, 1.2e5],
                  kind="kde", size=5, space=0)
plt.setp(p.ax_joint,
         xlabel='inheritance [atoms g-1]',
         ylabel='exposure time [yr]')
p.ax_joint.axhline(true_exposure, color='r')
p.ax_joint.axvline(true_inheritance, color='r')

# exposure vs. density
p = sns.jointplot(mcmc_smodel.trace('soil_density')[:],
                  mcmc_smodel.trace('exposure_time')[:],
                  xlim=[1.8, 2.3], ylim=[4e4, 1.2e5],
                  kind="kde", size=5, space=0)
plt.setp(p.ax_joint,
         xlabel='soil density [g cm-3]',
         ylabel='exposure time [yr]')
p.ax_joint.axhline(true_exposure, color='r')
p.ax_joint.axvline(true_density, color='r')

# density vs. erosion
p = sns.jointplot(mcmc_smodel.trace('erosion_rate')[:] * 1e4,
                  mcmc_smodel.trace('soil_density')[:],
                  xlim=[0, 15], ylim=[1.8, 2.3],
                  kind="kde", size=5, space=0)
plt.setp(p.ax_joint,
         xlabel='erosion rate [m Myr-1]',
         ylabel='soil density [g cm-3]')
p.ax_joint.axhline(true_density, color='r')
p.ax_joint.axvline(true_erosion * 1e4, color='r')

# inheritance vs. erosion
p = sns.jointplot(mcmc_smodel.trace('erosion_rate')[:] * 1e4,
                  mcmc_smodel.trace('inheritance')[:],
                  xlim=[0, 15], ylim=[0, 4000],
                  kind="kde", size=5, space=0)
plt.setp(p.ax_joint,
         xlabel='erosion rate [m Myr-1]',
         ylabel='inheritance [atoms g-1]')
p.ax_joint.axhline(true_inheritance, color='r')
p.ax_joint.axvline(true_erosion * 1e4, color='r')<jupyter_output><empty_output><jupyter_text>Plot the measured concentrations vs. depth and the predicted concentration profiles corresponding to the mean and specific quantiles of the PPD. The shaded grey area indicates the 95% HPD interval.<jupyter_code>sns.set(style="darkgrid")
sns.set_context('paper')

depths = np.linspace(profile_data['depth'].min(),
                     profile_data['depth'].max(),
                     100)

Cm_true = models.C_10Be(depths, true_erosion, true_exposure,
                        true_density, true_inheritance)

res_erosion_rate = erosion_rate.stats()
res_exposure_time = exposure_time.stats()
res_soil_density = soil_density.stats()
res_inheritance = inheritance.stats()
res_mprofile = mprofile.stats()

Cm_mean = models.C_10Be(depths,
                        res_erosion_rate['mean'],
                        res_exposure_time['mean'],
                        res_soil_density['mean'],
                        res_inheritance['mean'])

Cm_median = models.C_10Be(depths,
                          res_erosion_rate['quantiles'][50],
                          res_exposure_time['quantiles'][50],
                          res_soil_density['quantiles'][50],
                          res_inheritance['quantiles'][50])

plt.figure()
plt.errorbar(profile_data['C'],
             -profile_data['depth'],
             xerr=profile_data['std'],
             fmt='o', markersize=4,
             label='data')
plt.plot(Cm_mean, -depths, label='PPD mean')
plt.plot(Cm_median, -depths, label='PPD median')
plt.plot(Cm_true, -depths, color='r', label='true model')
plt.fill_betweenx(-profile_data['depth'],
                  res_mprofile['95% HPD interval'][:,0],
                  res_mprofile['95% HPD interval'][:,1],
                  color='k', alpha=0.3)
plt.setp(plt.gca(),
         xlabel='10Be concentration [atoms g-1]',
         ylabel='-1 * depth [cm]',
         xlim=[0, None], ylim=[None, 0])

plt.legend(loc='lower right')<jupyter_output><empty_output><jupyter_text>Show the correlation matrix of the free parameters (computed from the samples generated by MCMC).<jupyter_code>traces = np.column_stack((mcmc_smodel.trace('erosion_rate')[:],
                          mcmc_smodel.trace('exposure_time')[:],
                          mcmc_smodel.trace('soil_density')[:],
                          mcmc_smodel.trace('inheritance')[:]))

traces_df = pd.DataFrame(data=traces,
                         columns=['erosion rate', 'exposure time',
                                  'soil density', 'inheritance'])

f, ax = plt.subplots()
cmap = sns.blend_palette(["#00008B", "#6A5ACD", "#F0F8FF",
                          "#FFE6F8", "#C71585", "#8B0000"],
                         as_cmap=True)
sns.corrplot(traces_df, annot=True, sig_stars=False,
             diag_names=False, cmap=cmap, ax=ax)
f.tight_layout()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup