.. _computational_considerations:

Estimating Model Parameters
#############################
.. note::
    See :ref:`sum_link_dep_one_cmp_v_cl_stderr_tut` for the :ref:`tut_script` used here.

The previous chapters show how the observed concentrations are influenced by the parameters that define them such as the number of compartments and their function, how flows between compartments are specified, the shape of the dosing function, and the residual error model. We have explored these topics by taking models with given parameters, and using them to generate datasets of synthetic observations. This is often known as a *forward* problem.

We now turn to the more difficult *inverse* problem where we are given a dataset of measurements and a parameterized model, and want to estimate the most likely parameters for the given data. To solve this problem, we first need to make three design decisions.

First, we need a way to quantify the goodness-of-fit for a set of parameters, |ie| a single number that enables us to compare one set of parameters with another and determine which agrees better with the data. A popular choice for this is the :term:`likelihood` that quantifies how likely the data are, given the model parameters.

Second, we must specify an algorithm for finding the "best" model parameters from the set of possible options. In practice, finding *the* best set of parameters is difficult or impossible, so we instead look for a *locally* optimal set that is better than all other nearby sets [DennisSchnabel1987]_ [NocedalWright2006]_.

Finally, we need a set of criteria that determine whether we have reached the best local solution, |ie| whether the algorithm has converged.


Likelihood and the Objective Function
==========================================
Before we can find the "best" fit of a model to a set of observations, we first need a measure of what is "good". A common measure of goodness-of-fit is a model's likelihood [Millar2011]_: how likely are the observations to arise if the given set of model parameters is correct? Although this is *not* a probability (the integral over model parameters does not equal one), for a given set of observations and a given model structure a higher likelihood indicates a better fit.

In practice, we compute the Maximum Likelihood by minimizing the *negative* of the likelihood (because most optimization algorithms are written to find the minimum of an error or cost function). Moreover, where the likelihood is a member of the `exponential family of distributions <https://en.wikipedia.org/wiki/Exponential_family>`_ it is mathematically convenient to minimize *-2 log(likelihood)* which is commonly referred to as the :term:`objective function`, |objv|.

Using this measure of cost, we can take a "bird's eye view" of the surface in two dimensions by plotting it as a contour map for pairs of parameter values. In :ref:`dep_one_cmp_cl`, for example, we looked at a model with three parameters: *KA*, *CL* and *V*. Fixing *V* at its true value - which we know because these observations are synthetic - we can see how the cost surface varies with respect to *KA* and *CL* (:numref:`fig_dep_one_cmp_surface`).

.. figure:: /_autogen/surfaces/dep_one_cmp_v_cl_stderr.pyml_output/dep_one_cmp_v_cl_stderr_fit.pyml_output/objv_contour.*
    :name: fig_dep_one_cmp_surface
    :align: center
    :width: 75%
    
    *KA* |vs| *CL* surface plot of the Objective Function Value (|objv|) for a bolus dose, administered to a one compartment model with absorption

In this example, we see that the cost function has a single minimum in this region and is at its lowest close to the true values: *KA = 0.3* and *CL = 3*. Its lowest point is close to, not exactly at, the true value because of noise added to the observations.


Minimization Algorithm
=========================
Once we have specified the cost function we want to minimize, we need an algorithm that will find the minimum. Because some models have large datasets and complex models with many parameters, it is not practical to find the minimum via an exhaustive search over every possible combination of parameters.

We therefore adopt an incremental approach that takes an initial guess at the parameter values, then looks for a different set of values nearby that has a lower cost. Repeating this process takes us on a path across the likelihood surface until there are no nearby solutions that are better than the current one. This, by definition, is the *local* minimum.

.. figure:: /_autogen/surfaces/dep_one_cmp_v_cl_stderr.pyml_output/dep_one_cmp_v_cl_stderr_fit.pyml_output/objv_contour_path.*
    :name: fig_dep_one_cmp_optpath
    :align: center
    :width: 75%
    
    KA |vs| CL surface plot overlaid with the path taken by the optimization algorithm to the minimum
    
Depending on the shape of the cost surface and the particular algorithm used, finding the local minimum can take only a few steps or it can require many steps over a winding path. In the example given (:numref:`fig_dep_one_cmp_optpath`), a local minimum was found by |popy| in a single step. |popy|'s search algorithm is very good at finding nearby minima with very few steps, however it is impossible to guarantee that this is the global minimum. This is true for all search algorithms.

Details of optimization algorithms are beyond the scope of this manual, but there are many good references for the interested reader [DennisSchnabel1987]_ [NocedalWright2006]_.


Initial Values and Convergence
===================================
When using an local optimization algorithm, the current estimate heads "down the hill" to a nearby point where the cost is lower until there is nowhere lower to go. In some cases, there may be more than one point that is *locally* the lowest point. As a result, starting the algorithm from some points will lead to one local minimum whereas starting from other initial points will lead to other local minima.

The set of points that lead to a given local minimum are known as its :term:`basin of convergence`. In simple models, this is rarely a major problem. However, in complex models or if the fitted parameters are infeasible it is worthwhile to test alternative starting values, especially if the initial guess at a starting estimate was a long way from the fitted value.
