Profile Likelihood

Profile likelihood toolbox

Installation

Using pip

pip install profile_likelihood

From source

git clone https://git.physics.byu.edu/yonatank/profile_likelihood.git
pip install ./profile_likelihood  # Install the package

Theory

In the problem of fitting a theoretical model $f\left(t_m, \vec{\theta}\right)$ to the $M$ experimentally determined data points $y_m$ at times $t_m$, by assuming that the experimental errors for the data points are independent and Gaussian distributed with standard deviation of $\sigma$, the probability that a given model produced the observed data points is

    P\left(\vec{y}\middle|\vec{\theta}\right) =
		\prod_{m=1}^N \frac{1}{\sqrt{2\pi}\sigma}
        e^{-\left(f(t_m,\vec{\theta}) - y_m\right)^2/2\sigma^2}.

The likelihood function of this model, $L\left(\vec{\theta}\middle|\vec{y}\right)$, is the probability of the occurrence of the outcomes $\vec{y}$ given a set of parameters $\vec{\theta}$ of the model, $P\left(\vec{y}\middle|\vec{\theta}\right)$. Using the equation above, we can write

    L\left(\vec{\theta}\middle|\vec{y}\right) \propto
		\exp\left[ -C\left(\vec{\theta}\right) \right],

where $C\left(\vec{\theta}\right)$ is the cost function, given by

    C\left(\vec{\theta}\right) = \frac{1}{2} \sum_m
        \left(\frac{y_m - f\left(t_m, \vec{\theta}\right)}
			{\sigma_m} \right)^2

Suppose the model $f\left(t_m, \vec{\theta}\right)$ has $N$ parameters, written as $\{ \theta_0, \cdots, \theta_{N-1} \}$. The profile likelihood of the model for parameter $\theta_j$ is the possible maximum likelihood given the parameter $\theta_j$. The profile likelihood for parameter $\theta_j$ is calculated by setting $\theta_j$ to a fixed value, then maximizing the likelihood function (by minimizing the cost function) over the other parameters of the model. We repeat this computation across a range of $\theta_j$, $\left(\theta_j^{\min}, \theta_j^{\max}\right)$.

Basic usage

As an example, we want to run profile likelihood calculation to sum of exponential model

    f(\vec{\phi}; t) = \exp(-\phi_0*t) + \exp(-\phi_1*t),

with t = [0.5, 1.1, 3.0] and data y = [1.089, 0.717, 0.332]. We further want to constrain the parameters to be positive, which suggest to work with $\theta = \log{\phi}$ instead.

import numpy as np
from profile_likelihood import profile_likelihood

tdata = np.array([0.5, 1.1, 3.0]).reshape((-1, 1))  # Sampling time
ydata = np.array([1.089, 0.717, 0.332])  # Mock data

def residuals(theta):
	"""The function, representing the model, that will be used in the
	optimization process. This function takes an array of parameter and return
	an array of the residuals.
	"""
	phi = np.exp(theta)  # Undo the log transformation
	pred = np.sum(np.exp(-phi * tdata), axis=1)  # Model that gives predictions
	return ydata - pred  # Return the residual vector.

best_fit = np.array([-1.0, 1.0])  # Best-fit parameter
nparams = len(best_fit)  # Number of parameters
npred = 3  # Number of predictions

# Create likelihood object
pl = profile_likelihood(residuals, nparams, npred)

# Run computation
pl.compute(best_fit)

# Access results
pl.results

For more examples, see here.

List of frequently used methods

• profile_likelihood.compute
• profile_likelihood.save_results
• profile_likelihood.save_best_results
• profile_likelihood.load_results
• profile_likelihood.plot_likelihoods
• profile_likelihood.plot_paths
• profile_likelihood.plot_likelihoods_and_paths

Contact

Yonatan Kurniawan (kurniawanyo@outlook.com)