TUTORIAL_MCMC_sampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Markov Chain Monte Carlo Sampling (insoluble monolayer example)\n",
    "\n",
    "Markov Chain Monte Carlo (MCMC) sampling is a well-established method of using Bayes' theorem to determine estimates for the uncertainty in model input parameters. MultilayerPy employs the `emcee` python package, which is widely used for this purpose (https://doi.org/10.1086/670067). The full up-to-date documentation for emcee is found here: https://emcee.readthedocs.io/en/stable/index.html.\n",
    "\n",
    "There is an initial example which will demonstrate how to set up the system and it will run into a common problem. An improvement to the system is then presented along with a summary of what we can learn from an MCMC analysis.\n",
    "\n",
    "### Summary of MCMC\n",
    "MCMC sampling essentially involves initiating an ensemble of \"walkers\" (parameter sets) in a pre-defined parameter space. These walkers are also referred to as chains. At each iteration of the algorithm a new position is proposed for each walker. The probability that this next position is accepted is dependent on the probability of the previous step (i.e. the goodness of fit). Each iteration builds up a chain of parameter values which eventually tends towards the region of maximum likelihood. As a result, after a period of \"burn-in\" time where the chains find and converge around the region of maximum likelihood, the chains find an equilibrium position and walk around this region. A histogram can then be plotted for each model parameter, revealing each parameter's probability distribution (for that model-experimental data system) - this is not always gaussian.  \n",
    "\n",
    "The user is encouraged to find out more in the detailed documentation and paper referred to above. \n",
    "\n",
    "The model system we will use here is identical to the system studied in the \"insoluble monolayers\" tutorial notebook. In practice MCMC works best when initialising the ensemble of \"walkers\" around what is most likely the global minimum (hence the initial global optimisation). The focus will therefore be on the MCMC sampling procedure carried out after the initial global optimisation.\n",
    "\n",
    "The dataset used is from Woden et al. (2021) (https://doi.org/10.5194/acp-21-1325-2021).\n",
    "\n",
    "This dataset is from a monolayer of deuterated oleic acid on an aqueous subphase. We will assume that none of the products evaporate from the surface and keep them lumped together as a \"products\" component for simplicity. \n",
    "\n",
    "The reaction scheme is consistent with other MultilayerPy tutorials: `oleic acid + ozone --> products`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "# this stops numpy paralellising some processes\n",
    "# when it comes to parallel MCMC sampling we don't want to parallelise sub-processes (slows down otherwise)\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"1\"\n",
    "\n",
    "# importing the necessaries\n",
    "import numpy as np\n",
    "import multilayerpy.build as build \n",
    "import multilayerpy.simulate as simulate\n",
    "import multilayerpy.optimize as optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the ModelType class\n",
    "from multilayerpy.build import ModelType\n",
    "\n",
    "# import the ReactionScheme class\n",
    "from multilayerpy.build import ReactionScheme\n",
    "\n",
    "# define the model type (KM-SUB in this case) and geometry (film)\n",
    "mod_type = ModelType('km-sub','film')\n",
    "\n",
    "# build the reaction tuple list, in this case only 1 tuple in the list (for 1 reaction)\n",
    "# component 1 (oleic acid) reacts with component 2 (ozone)\n",
    "reaction_tuple_list = [(1,2)]\n",
    "\n",
    "# build the product tuple list, only component 3 (products) is a product\n",
    "# a tuple with a single value inside is defined (value,)\n",
    "product_tuple_list = [(3,)]\n",
    "\n",
    "# now construct the reaction scheme\n",
    "# we can give it a name and define the nuber of components as below\n",
    "reaction_scheme = ReactionScheme(mod_type,name='oleic acid ozonolysis',\n",
    "                                                   reactants=reaction_tuple_list,\n",
    "                                                products=product_tuple_list)\n",
    "\n",
    "# let's print out a representation of the reaction scheme\n",
    "reaction_scheme.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import ModelComponent class\n",
    "from multilayerpy.build import ModelComponent\n",
    "\n",
    "# making model components\n",
    "\n",
    "# oleic acid\n",
    "OA = ModelComponent(1,reaction_scheme,name='Oleic acid')\n",
    "\n",
    "# ozone, declare that it is a gas\n",
    "O3 = ModelComponent(2,reaction_scheme,gas=True,name='Ozone') \n",
    "\n",
    "# products\n",
    "prod = ModelComponent(3,reaction_scheme, name='Reaction products')\n",
    "\n",
    "# collect into a dictionary\n",
    "model_components_dict = {'1':OA,\n",
    "                        '2':O3,\n",
    "                        '3':prod}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import DiffusionRegime class\n",
    "from multilayerpy.build import DiffusionRegime\n",
    "\n",
    "# making the diffusion dictionary\n",
    "diff_dict = {'1' : None,\n",
    "             '2': None,\n",
    "             '3':None}  \n",
    "\n",
    "# make diffusion regime\n",
    "diff_regime = DiffusionRegime(mod_type,model_components_dict,diff_dict=diff_dict)\n",
    "\n",
    "# call it to build diffusion code ready for the builder\n",
    "diff_regime()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import ModelBuilder class\n",
    "from multilayerpy.build import ModelBuilder\n",
    "\n",
    "# create the model object, ignore [1,2,3] etc at the end\n",
    "model = ModelBuilder(reaction_scheme,model_components_dict,diff_regime)\n",
    "\n",
    "# build the model. Will save a file, don't include the date in the model filename\n",
    "model.build(date=False)\n",
    "\n",
    "# print out the parameters required for the model to run\n",
    "print(model.req_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making an insoluble monolayer\n",
    "\n",
    "An insoluble monolayer can be made by setting the bulk diffusion coefficient of that component = 0.0 cm2 s-1. This means that that component does not diffuse at all in the bulk. We will nominally create 5 model bulk layers and a bulk film thickness of 2 µm to satisfy the model building process. However, there is essentially no exchange of material between the bulk and surface layers. Henry's law coefficient is also set to 0.0.\n",
    "\n",
    "The surface reaction rate constant determined by Woden et al is 2.2e-10 cm2 s-1. Ozone concentration was 323 ppb (323 x 2.46e10 cm-3). We will vary `alpha_s_0_2`, `Td_2` and `k_1_2_surf` in the model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# import the Simulate class\n",
    "from multilayerpy.simulate import Simulate\n",
    "\n",
    "# import the Parameter class\n",
    "from multilayerpy.build import Parameter\n",
    "\n",
    "# make the parameter dictionary\n",
    "# SETTING BULK DIFFUSION AND HENRY'S LAW PARAMETERS TO 0.0\n",
    "param_dict = {'delta_3':Parameter(1e-7),\n",
    "              'alpha_s_0_2':Parameter(2e-3,vary=True,bounds=(1e-4,1e-1)),\n",
    "              'delta_2':Parameter(0.4e-7),\n",
    "              'Db_2':Parameter(0.0),\n",
    "              'delta_1':Parameter(0.8e-7),\n",
    "              'Db_1':Parameter(0.0),\n",
    "              'Db_3':Parameter(0.0),\n",
    "              'k_1_2':Parameter(0.0),\n",
    "              'H_2':Parameter(0.0),\n",
    "              'Xgs_2': Parameter(323.0 * 2.46e10),\n",
    "              'Td_2': Parameter(1e-7,vary=True,bounds=(1e-9,1e-6)),\n",
    "              'w_2':Parameter(3.6e4),\n",
    "              'T':Parameter(294.0),\n",
    "              'k_1_2_surf':Parameter(2.2e-10)}\n",
    "# 'Td_2': Parameter(1e-7,vary=True,bounds=(1e-8,1e-6))\n",
    "\n",
    "# import the data\n",
    "from multilayerpy.simulate import Data\n",
    "\n",
    "raw_data = np.genfromtxt('woden_etal_acp_data_uncert.csv',delimiter=',')\n",
    "\n",
    "data = Data(raw_data)\n",
    "\n",
    "# make the simulate object with the model and parameter dictionary\n",
    "sim = Simulate(model,param_dict,data=data)\n",
    "\n",
    "# define required parameters\n",
    "n_layers = 5\n",
    "rp = 2e-4 # radius in cm\n",
    "time_span = [0,800] # in s\n",
    "n_time = 999 # number of timepoints to save to output\n",
    "\n",
    "#spherical V and A\n",
    "# use simulate.make_layers function\n",
    "V, A, layer_thick = simulate.make_layers(mod_type,n_layers,rp)\n",
    "\n",
    "# initial conc. of everything\n",
    "bulk_conc_dict = {'1':0.0,'2':0.0,'3':0.0} # key=model component number, value=bulk conc, units: cm-3\n",
    "\n",
    "# initial surf conc of oleic acid calculated from neutron reflectometry model fit\n",
    "surf_conc_dict = {'1':262571428571428.56,'2':0.0,'3':0.0} # key=model component number, value=surf conc, units: cm-2\n",
    "\n",
    "y0 = simulate.initial_concentrations(mod_type,bulk_conc_dict,surf_conc_dict,n_layers) \n",
    "    \n",
    "# now run the model\n",
    "output = sim.run(n_layers,rp,time_span,n_time,V,A,layer_thick,Y0=y0)\n",
    "\n",
    "%matplotlib inline\n",
    "# plot the model\n",
    "sim.plot(norm=True)\n",
    "\n",
    "# There may be some runtime warnings because we are forcing some divisions by 0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now optimise the model\n",
    "\n",
    "# import the optimize module and Optimizer object\n",
    "import multilayerpy.optimize\n",
    "from multilayerpy.optimize import Optimizer\n",
    "\n",
    "fitter = Optimizer(sim)\n",
    "\n",
    "res = fitter.fit(method='least_squares',weighted=False);\n",
    "\n",
    "\n",
    "sim.plot(norm=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC sampling\n",
    "\n",
    "Now we can set up MCMC sampling. We need to define the number of walkers (chains) we want to use. The more the merrier. Normally a multiple of the number of CPUs available in your system. Use `os.cpu_count()` to find this out if you don't know already. My machine has 12 cores available. \n",
    "\n",
    "Here, I have selected 120 walkers and 1000 samples. We will decide how many of the initial steps to discard (or burn) after the sampling procedure. \n",
    "\n",
    "In order to paralellise the algorithm, the `multiprocessing` package is used and the `Pool` object is imported. Using the `with Pool() as pool:` syntax will ensure parallelisation happens without unnecessary problems. \n",
    "\n",
    "The `sampler` defined below is an `emcee.EnsembleSampler` object and can be manipulated as described in the emcee documentation (https://emcee.readthedocs.io/en/stable/user/sampler/). \n",
    "\n",
    "**NOTE** if you are carrying this out in an IDE such as Spyder on a Windows system, you may need to invoke the `if __name__ == \"__main__\":` syntax before typing the code outlined below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MCMC sampling\n",
    "walkers = 120\n",
    "samples = 1000\n",
    "\n",
    "# set the numpy random seed so that the analysis is reproducible\n",
    "np.random.seed(1)\n",
    "\n",
    "from multiprocessing import Pool\n",
    "\n",
    "with Pool() as pool:\n",
    "\n",
    "    sampler = fitter.sample(n_walkers=walkers,samples=samples,pool=pool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to look at the chains that have been generated and plot what they look like. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the chains (see emcee documentation for details)\n",
    "chains = sampler.get_chain()\n",
    "\n",
    "# printing the shape of the chains\n",
    "print('chains = (number of samples, number of walkers , number of dimensions)',chains.shape)\n",
    "\n",
    "# we can plot the chains with the Optimizer object (fitter in this case)\n",
    "fitter.plot_chains()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's decide how many initial steps to burn\n",
    "# the chains seem to have expanded to their equilibrium distributions after ~ 800 steps\n",
    "\n",
    "burn_no = 800\n",
    "thin = 1\n",
    "\n",
    "# redefine chains now discarding the burn-in number of samples\n",
    "chains = sampler.get_chain(discard=burn_no,thin=thin)\n",
    "print(chains.shape)\n",
    "\n",
    "# plot the chains again \n",
    "fitter.plot_chains(discard=burn_no)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the sampler for more steps after the burn-in stage will allow us to get a more accurate estimation of the parameter probability distribution.\n",
    "\n",
    "The idea of autocorrelation and thinning will not be discussed here but MCMC sampling studies should ideally be accompanied by an analysis of autocorrelation (see references at the start of this tutorial). There should be at least the details regarding the number of samples, walkers, burn-in steps and parameter bounds for someone to reproduce the analysis. Better yet, a Jupyter notebook going through the analysis! \n",
    "\n",
    "Below is one method of how to present the distributions as a corner plot using the `corner` module. If you do not already have this module `$ pip install corner` will install the package into your python environment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the flattened chains (chains from all walkers joined together)\n",
    "flat_chains = sampler.get_chain(flat=True,discard=burn_no,thin=thin)\n",
    "\n",
    "import corner\n",
    "\n",
    "# order of labels is the same as the order of varying parameters printed during MCMC sampling\n",
    "labels = ['alpha_s_0_2', 'Td_2']\n",
    "\n",
    "# corner.corner is the best method for this \n",
    "# we can plot 25th and 75th quantiles\n",
    "fig = corner.corner(\n",
    "    flat_chains, quantiles=[0.75,0.25],\n",
    "    bins=20,\n",
    "    labels=labels\n",
    "    \n",
    ");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What could this mean?\n",
    "\n",
    "Ideally there would be a nice gaussian blob in the middle of this corner plot. But there isn't. These two parameters are highly correlated (with an interesting banana shape). We could still quote the mean of these parameter probability distributions but it would be wise to present this plot along with any analysis. Most uncertainties are assumed to be the standard deviation of some gaussian distribution. The shape of these distributions is not gaussian. This is an advantage of this kind of analysis.\n",
    "\n",
    "Take-home message from this:\n",
    "1. Either `alpha_s_0_2` or `Td_2` should be held constant and/or experimentally constrained in some way. They are both associated with the adsorption and desorption of ozone to the surface. \n",
    "2. Could `alpha_s_0_2` and/or `Td_2`be dependent on the composition of the monolayer? See the \"composition-dependent surface adsorption\" tutorial for how to construct such a model. \n",
    "\n",
    "**NOTE** as the number of model components and layers increases, so does the time it takes for a single model run. This would significantly increase the time taken to optimise and sample the model-data system. Paralellising the global optimisation and sampling procedure is recommended for this (possible with MutilayerPy). At some point, however, a lot of CPUs would be required for bulky models. A high-performance computing resource (such as a computer cluster) would be recommended for this. \n",
    "\n",
    "## Plotting the outcome\n",
    "\n",
    "We can take a random `n_samples` number of samples from the MCMC run to plot. Firstly, because the optimiser didn't save each model output as the MCMC sampling procedure progressed (see note below), we need to run the model with this sub-sample of model parameters. \n",
    "\n",
    "MultilayerPy has the utility to do this. Let's grab 50 random samples from the MCMC sampling run, run and save those model outputs and plot them in 3 lines. \n",
    "\n",
    "**Note** the optimiser did not save each model output after each iteration due to memory concerns. For example, running 120 walkers for 1000 steps produces 120000 samples. The amount of data to store quickly ramps up. Therefore it is appropriate to store only the runs we want for presentation. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples_plot = 200\n",
    "fitter.get_chain_outputs(n_samples=n_samples_plot,n_burn=burn_no,override_stop_run=True)\n",
    "fitter.plot_chain_outputs()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is not a good outcome. The MCMC chains have not converged fully and more sampling would be required. You can see that some of the sampled runs are nowhere near the data. An analysis of the autocorrelation time would also help. This is not the focus of this notebook but interested readers can consult the references at the start of this notebook. Holding one of these highly correlated parameters constant may be a good way to proceed...\n",
    "\n",
    "## An improvement\n",
    "\n",
    "Let's hold `Td_2` at 1e-7 s. Which is a reasonable guess for ozone at room temperature. We'll only allow `alpha_s_0_2` to vary. We can do the same optimisation + MCMC sampling procedure as before (to save time, we'll take fewer samples). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make the parameter dictionary\n",
    "# SETTING BULK DIFFUSION AND HENRY'S LAW PARAMETERS TO 0.0\n",
    "param_dict = {'delta_3':Parameter(1e-7),\n",
    "              'alpha_s_0_2':Parameter(1e-3,vary=True,bounds=(1e-4,1e-2)),\n",
    "              'delta_2':Parameter(0.4e-7),\n",
    "              'Db_2':Parameter(0.0),\n",
    "              'delta_1':Parameter(0.8e-7),\n",
    "              'Db_1':Parameter(0.0),\n",
    "              'Db_3':Parameter(0.0),\n",
    "              'k_1_2':Parameter(0.0),\n",
    "              'H_2':Parameter(0.0),\n",
    "              'Xgs_2': Parameter(323.0 * 2.46e10),\n",
    "              'Td_2': Parameter(1e-7),\n",
    "              'w_2':Parameter(3.6e4),\n",
    "              'T':Parameter(294.0),\n",
    "              'k_1_2_surf':Parameter(2.2e-10)}\n",
    "\n",
    "\n",
    "# import the data\n",
    "from multilayerpy.simulate import Data\n",
    "raw_data = np.genfromtxt('woden_etal_acp.csv',delimiter=',')\n",
    "raw_errors = np.genfromtxt('woden_etal_errors.csv',delimiter=',')\n",
    "\n",
    "actual_errors = (raw_errors[:,1] - raw_data[:,1]) * 2\n",
    "\n",
    "collected_data = np.column_stack((raw_data,actual_errors))\n",
    "\n",
    "data = Data(collected_data[:6,:])\n",
    "\n",
    "# make the simulate object with the model and parameter dictionary\n",
    "sim = Simulate(model,param_dict,data=data)\n",
    "\n",
    "# define required parameters\n",
    "n_layers = 5\n",
    "rp = 2e-4 # radius in cm\n",
    "time_span = [0,800] # in s\n",
    "n_time = 999 # number of timepoints to save to output\n",
    "\n",
    "#spherical V and A\n",
    "# use simulate.make_layers function\n",
    "V, A, layer_thick = simulate.make_layers(mod_type,n_layers,rp)\n",
    "\n",
    "# initial conc. of everything\n",
    "bulk_conc_dict = {'1':0.0,'2':0.0,'3':0.0} # key=model component number, value=bulk conc, units: cm-3\n",
    "\n",
    "# initial surf conc of oleic acid calculated from neutron reflectometry model fit\n",
    "surf_conc_dict = {'1':262571428571428.56,'2':0.0,'3':0.0} # key=model component number, value=surf conc, units: cm-2\n",
    "\n",
    "y0 = simulate.initial_concentrations(mod_type,bulk_conc_dict,surf_conc_dict,n_layers) \n",
    "    \n",
    "# now run the model\n",
    "output = sim.run(n_layers,rp,time_span,n_time,V,A,layer_thick,Y0=y0)\n",
    "\n",
    "%matplotlib inline\n",
    "# plot the model\n",
    "sim.plot(norm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now optimise the model\n",
    "\n",
    "# import the optimize module and Optimizer object\n",
    "import multilayerpy.optimize\n",
    "from multilayerpy.optimize import Optimizer\n",
    "\n",
    "fitter = Optimizer(sim)\n",
    "\n",
    "res = fitter.fit(method='differential_evolution',popsize=20);\n",
    "\n",
    "\n",
    "sim.plot(norm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MCMC sampling\n",
    "walkers = 120\n",
    "samples = 200\n",
    "\n",
    "# set the numpy random seed so that the analysis is reproducible\n",
    "np.random.seed(2)\n",
    "\n",
    "from multiprocessing import Pool\n",
    "\n",
    "with Pool() as pool:\n",
    "\n",
    "    sampler = fitter.sample(n_walkers=walkers,samples=samples,pool=pool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "burn_no = 100\n",
    "thin = 1\n",
    "\n",
    "# redefine chains now discarding the burn-in number of samples\n",
    "chains = sampler.get_chain(discard=burn_no,thin=thin)\n",
    "print(chains.shape)\n",
    "\n",
    "# plot the chains again \n",
    "fitter.plot_chains(discard=burn_no)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the flattened chains (chains from all walkers joined together)\n",
    "flat_chains = sampler.get_chain(flat=True,discard=burn_no,thin=thin)\n",
    "\n",
    "import corner\n",
    "\n",
    "# order of labels is the same as the order of varying parameters printed during MCMC sampling\n",
    "labels = ['alpha_s_0_2']\n",
    "\n",
    "# corner.corner is the best method for this \n",
    "# using quantiles allows you to define a line for a confidence interval for your data\n",
    "fig = corner.corner(\n",
    "    flat_chains, quantiles=[0.975,0.025],\n",
    "    bins=20,\n",
    "    labels=labels)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n_samples_plot = 50\n",
    "chain_outputs = fitter.get_chain_outputs(n_samples=n_samples_plot,n_burn=burn_no)\n",
    "fitter.plot_chain_outputs()\n",
    "\n",
    "# lets print the 2.5 and 97.5 % quantiles for this data (an uncertainty range to quote)\n",
    "upper, lower = np.quantile(flat_chains,[0.975,0.025])\n",
    "mean = np.mean(flat_chains)\n",
    "std = np.std(flat_chains)\n",
    "print(f'mean alpha_s_0: {mean}')\n",
    "print(f'std alpha_s_0: {std}')\n",
    "print(f'2.5% quantile: {lower}')\n",
    "print(f'97.5% quantile: {upper}')\n",
    "\n",
    "# let's save the optimised model output to a .csv file. A summary of what we have done here\n",
    "sim.save_params_csv()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ah, much better! The histogram for `alpha_s_0_2` is looking nice and gaussian and the sample of model outputs seem to agree well with the data. \n",
    "\n",
    "### What can we learn from this?\n",
    "* MCMC sampling can be a powerful tool in inferring parameters from your model-data system.\n",
    "* Care must be taken when selecting the parameters to vary. Some can be highly correlated (as we saw here). Can you constrain one of them?\n",
    "* The number of walkers, samples, burn-in steps and thinning steps matter. There is a large body of information available about the best way to carry out an MCMC sampling study. MultilayerPy uses the well-established `emcee` package to facilitate this but caution is needed when drawing conclusions from such an analysis. Quote what you did! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inter-quartile range\n",
    "print(np.quantile(flat_chains,[0.25,0.75]))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}