************
MCMC Fitting
************

Another fitting method available in :meth:`CRModel.fit() <mgrowthctrl.models.crm.model.CRModel.fit>` is Markov Chain Monte Carlo (MCMC). For the least-squares approach, see :doc:`lsq`.

The :meth:`CRModel.fit_mcmc() <mgrowthctrl.models.crm.model.CRModel.fit_mcmc>` method performs Bayesian parameter estimation using MCMC sampling via the ``pymcmcstat`` library.

Background
==========

To quantify uncertainty in the fitted parameters, we sample from the posterior distribution. As in :doc:`lsq`, the difference between model predictions and observations is measured by the sum of squared scaled errors (SSE),

.. math::

   \mathrm{SSE}(\theta)
   = \sum_i \sum_k
   \left(\frac{\hat{X}_i(t_k;\theta) - X_i(t_k)}{s_{X,i}}\right)^2
   + \sum_j \sum_k
   \left(\frac{\hat{S}_j(t_k;\theta) - S_j(t_k)}{s_{S,j}}\right)^2\,,

where :math:`\theta` denotes the vector of model parameters, :math:`t_k` are the measurement times, :math:`X_i(t_k)` and :math:`S_j(t_k)` are the observed concentrations, :math:`\hat{X}_i(t_k;\theta)` and :math:`\hat{S}_j(t_k;\theta)` are the corresponding model predictions, and :math:`s_{X,i}` and :math:`s_{S,j}` are scaling factors for species and metabolite residuals, respectively.

The sampler combines this likelihood with weakly informative priors on the parameters to generate posterior samples, and an initial burn-in period is discarded before analysis.

Example
=======

We load experimental time-series data from the `BT_WC_export.csv <https://gitlab.com/ComputationalScience/mgrowthctrl/-/raw/main/examples/datasets/BT_WC_export.csv>`_ file and fit a consumer-resource model (CRM) to the data.

.. code-block:: python

    from mgrowthctrl.utils.data import Dataloader
    from mgrowthctrl.models.crm.model import CRModel

    data = Dataloader()
    data.load_local_data(
        "examples/datasets/BT_WC_export.csv",
        x_selector=r"BT counts",
    )

    model = CRModel.from_single_species_data(
        df=data.df,
        time_col=data.time_col,
        x_col=data.X_names[0],
        s_cols=data.S_names,
    )


We then fit the model with least squares to obtain a good initial parameter set.

.. code-block:: python

    model.fit(
        df=data.df,
        time_col=data.time_col,
        x_cols=data.X_names,
        s_cols=data.S_names,
    )

Next, we use that fit as a starting point for MCMC sampling to quantify uncertainty.

.. code-block:: python

    result_mcmc = model.fit_mcmc(
        nsimu=10_000,
        sigma2=12.0,
        burn=0.8,
        seed=1234,
    )

    mcmc_sims = result_mcmc["mcmc_sims"]
    sschain = result_mcmc["sschain"]
    post = result_mcmc["post"]
    param_names = result_mcmc["param_names"]

Finally, we visualize the MCMC error evolution, parameter distributions, and model trajectories.

.. code-block:: python

    from mgrowthctrl.utils.plot import (
        plot_data_with_overlay,
        plot_mcmc_error_evolution,
        plot_parameter_distributions,
    )

    plot_mcmc_error_evolution(sschain, save_prefix="mcmc_error")

    plot_parameter_distributions(
        post,
        theta_mle=model.theta_opt_,
        names=param_names,
        save_prefix="param_dist",
    )

    plot_data_with_overlay(
        df=data.df,
        sim=model.sim_opt_,
        overlay="model",
        time_col=data.time_col,
        x_cols=data.X_names,
        s_cols=data.S_names,
        mcmc_sims=mcmc_sims,
        thin_alpha=0.1,
        thin_lw=0.5,
        thick_lw=2.5,
        shade_percentiles=True,
        percentiles=(25, 75),
        save_prefix="final_mcmc_fit",
    )

The MCMC error evolution should look similar to the following example:

.. image:: ../_static/mcmc_error.png
   :alt: MCMC error evolution
   :align: center
   :width: 50%

The parameter distributions should look similar to the following example:

.. image:: ../_static/param_dist.png
   :alt: MCMC parameter distributions
   :align: center
   :width: 80%

The fitted model trajectories should look similar to the following example:

.. image:: ../_static/crm_mcmc_example.png
   :alt: CRM MCMC example
   :align: center

Backend API
===========

.. autofunction:: mgrowthctrl.models.crm.fit.crm_fit_mcmc
    :noindex:

.. autofunction:: mgrowthctrl.models.crm.fit.crm_fit_mcmc_from_df