# Source code for curveball.likelihood

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# This file is part of curveball.
# https://github.com/yoavram/curveball

# Copyright (c) 2015, Yoav Ram <yoav@yoavram.com>
from __future__ import print_function
from __future__ import division
from builtins import str
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
sns.set_style("ticks")
import curveball
import curveball.baranyi_roberts_model

[docs]def loglik(t, y, y_sig, f, penalty=None, **params):
r"""Computes the log-likelihood of seeing the data given a model
assuming normal distributed observation/measurement errors.

.. math::

\log{L(y | \theta)} = -\frac{1}{2} \sum_i { \log{(2 \pi \sigma_{i}^{2})} + \frac{(y - f(t_i; \theta))^2}{\sigma_{i}^{2}} }

which is the log-likelihood of seeing the data points :math:t_i, y_i
with measurement error :math:\sigma_i
given the model function :math:f, the model parameters :math:\theta,
and that the measurement error at time :math:t_i has a normal distribution with mean 0.

Parameters
----------
t : np.ndarray
one dimensional array of time
y : np.ndarray
one dimensional array of the means of the observations
y_sig : np.ndarray
one dimensional array of standrad deviations of the observations
f : callable
a function the calculates the expected observations (f(t)) from t and any parameters in params
penalty : callable
a function that calculates a scalar penalty from the parameters in params to be substracted from the log-likelihood
params : floats, optional
model parameters

Returns
-------
float
the log-likelihood result
"""
yhat = f(t, **params)
val = -0.5 * np.sum(np.log(2 * np.pi * y_sig ** 2) + (y - yhat) ** 2 / y_sig ** 2)
if penalty:
val -= penalty(**params)
return val

[docs]def ridge_regularization(lam, **center):
r"""Create a penaly function that employs the ridge regularization method:

.. math::

P = \lambda ||\theta - \theta_0||_2

where :math:\lambda is the regularization scale,
:math:\theta is the model parameters vector,
and :math:\theta_0 is the model parameters guess vector.
This is similar to using a multivariate Gaussian prior distribution on the model parameters
with the Gaussian centerd at :math:\theta_0 and scaled by :math:\lambda.

Parameters
----------
lam : float
the penalty factor or regularization scale
center : floats, optional
guesses of model parameters

Returns
-------
callable
the penalty function, accepts model parameters as float keyword arguments and returns a float penalty to the log-likelihood

Examples
--------
>>> penalty = ridge_regularization(1, y=0.1, K=1, r=1)
>>> loglik(t, y, y_sig, logistic, penalty=penalty, y0=0.12, K=0.98, r=1.1)
"""
def _ridge_regularization(**params):
return lam * np.linalg.norm([v - center.get(k, 0) for k, v in params.items() if np.isfinite(v)])
return _ridge_regularization

[docs]def loglik_r_nu(r_range, nu_range, df, f=curveball.baranyi_roberts_model.baranyi_roberts_function,
penalty=None, **params):
r"""Estimates the log-likelihood surface for :math:r and :math:\nu given data and a model function.

Parameters
----------
r_range, nu_range : numpy.ndarray
vectors of floats of :math:r and :math:\nu values on which to compute the log-likelihood
df : pandas.DataFrame
data frame with Time and OD columns
f : callable, optional
model function, defaults to :py:func:curveball.baranyi_roberts_model.baranyi_roberts_function
penalty : callable, optional
if given, the result of penalty will be substracted from the log-likelihood for each parameter set
params : floats
values for the model model parameters used by f

Returns
-------
np.ndarray
two-dimensional array of log-likelihood calculations;
value at index i, j will have the value for r_range[i] and nu_range[j]

--------
loglik
loglik_r_q0
"""
if not params:
params = dict()
t = df.Time.unique()
y = df.groupby('Time').OD.mean().as_matrix()
y_sig = df.groupby('Time').OD.std().as_matrix()

output = np.empty((len(r_range), len(nu_range)))
for i, r in enumerate(r_range):
params['r'] = r
for j, nu in enumerate(nu_range):
params['nu'] = nu
output[i,j] = loglik(t, y, y_sig, f, penalty, **params)

return output

[docs]def loglik_r_q0(r_range, q0_range, df, f=curveball.baranyi_roberts_model.baranyi_roberts_function,
penalty=None, **params):
r"""Estimates the log-likelihood surface for :math:r and :math:\nu given data and a model function.

Parameters
----------
r_range, q0_range : numpy.ndarray
vectors of floats of :math:r and :math:q_0 values on which to compute the log-likelihood
df : pandas.DataFrame
data frame with Time and OD columns
f : callable, optional
model function, defaults to :py:func:curveball.baranyi_roberts_model.baranyi_roberts_function
penalty : callable, optional
if given, the result of penalty will be substracted from the log-likelihood for each parameter set
params : floats
values for the model model parameters used by f

Returns
-------
np.ndarray
two-dimensional array of log-likelihood calculations;
value at index i, j will have the value for r_range[i] and q0_range[j]

--------
loglik
loglik_r_nu
"""
if not params:
params = dict()
t = df.Time.unique()
y = df.groupby('Time').OD.mean().as_matrix()
y_sig = df.groupby('Time').OD.std().as_matrix()

output = np.empty((len(r_range), len(q0_range)))
for i, r in enumerate(r_range):
params['r'] = r
for j, q0 in enumerate(q0_range):
params['q0'] = q0
output[i,j] = loglik(t, y, y_sig, f, penalty, **params)

return output

[docs]def plot_loglik(Ls, xrange, yrange, xlabel=None, ylabel=None, columns=4, fig_title=None, normalize=True,
ax_titles=None, cmap='viridis', colorbar=True, ax_width=4, ax_height=4, ax=None):
r"""Plots one or more log-likelihood surfaces.

Parameters
----------
Ls : sequence of numpy.ndarray
list or tuple of log-likelihood two-dimensional arrays; if one array is given it will be converted to a size 1 list
xrange, yrange : np.ndarray
values on x-axis and y-axis of the plot (rows and columns of Ls, respectively)
xlabel, ylabel : str, optional
strings for x and y labels
columns : int, optional
number of columns in case that Ls has more than one matrice
fig_title : str, optional
a title for the whole figure
normalize : bool, optional
if :py:const:True, all matrices will be plotted using a single color scale
ax_titles : list or tuple of str, optional
titles corresponding to the different matrices in Ls
cmap : str. optional
name of a matplotlib colormap (to see list, call :py:func:matplotlib.pyplot.colormaps), defaults to viridis
colorbar : bool, optional
if :py:const:True a colorbar will be added to the plot
ax_width, ax_height : int
width and height of each panel (one for each matrice in Ls)
ax : matplotlib axes or numpy.ndarray of axes
if given, will plot into ax, otherwise will create a new figure

Returns
-------
fig : matplotlib.figure.Figure
figure object
ax : numpy.ndarray
array of axis objects

Examples
--------
>>> L = loglik_r_nu(rs, nus, df, y0=y0, K=K, q0=q0, v=v)
>>> plot_loglik(L0, rs, nus, normalize=False, fig_title=fig_title, xlabel=r'$r$', ylabel=r'$\nu$', colorbar=False)
"""
if not isinstance(Ls, (list, tuple)):
Ls = [Ls]
columns = min(columns, len(Ls))
rows = int(np.ceil(len(Ls) / columns))
if ax is None:
fig, ax = plt.subplots(rows, columns, sharex=True, sharey=True, figsize=(ax_width*columns, ax_height*rows))
else:
fig = ax.figure
if not hasattr(ax, '__iter__'):
ax = np.array(ax, ndmin=2)
if ax.ndim == 1:
ax.resize((rows, columns))
vmin = np.nanmin(Ls)
vmax = np.nanmax(Ls)
norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)

for i, L in enumerate(Ls):
row = i // columns
col = i % columns
_ax = ax[row, col]

im = _ax.imshow(L.T, cmap=cmap, aspect=1, origin=(0,0))
if normalize:
im.set_norm(norm)

if colorbar:
plt.colorbar(im, ax=_ax, label='Log Likelihood', fraction=0.03, pad=0.25, format='%1.e')

_ax.xaxis.grid(color='k', ls='--', alpha=0.5)
_ax.yaxis.grid(color='k', ls='--', alpha=0.5)

xstep = len(xrange) // 5
xticks = list(range(0, len(xrange), xstep)) + [len(xrange)-1]
xticklabels = xrange[xticks]
_ax.set_xticks(xticks)
_ax.set_xticklabels(['{:.2g}'.format(x) for x in xticklabels], rotation='vertical')

ystep = len(xrange) // 5
yticks = list(range(0, len(yrange), ystep)) + [len(yrange)-1]
yticklabels = yrange[yticks]
_ax.set_yticks(yticks)
_ax.set_yticklabels(['{:.2g}'.format(y) for y in yticklabels])

imax, jmax = (L == np.nanmax(L)).nonzero()
_ax.scatter(imax, jmax, marker='o', s=50, color='r')

if row == rows-1 and xlabel:
_ax.set_xlabel(xlabel)
if col == 0 and ylabel:
_ax.set_ylabel(ylabel)
if ax_titles:
_ax.set_title(ax_titles[i])

if fig_title:
fig.text(0.5, 1, fig_title, fontsize=24, horizontalalignment='center')
fig.tight_layout()
sns.despine()
return fig, ax

[docs]def plot_model_loglik(m, df, fig_title=None):
r"""Plot the log-ikelihood surfaces for :math:\nu over :math:r and :math:q_0 over :math:r for given data and model fitting result.

Parameters
----------
m : lmfit.model.ModelResult
model for which to plot the log-likelihood surface
df : pandas.DataFrame
data frame with Time and OD columns used to fit the model
fig_title : str
title for the plot

Returns
-------
fig : matplotlib.figure.Figure
figure object
ax : numpy.ndarray
array of axis objects

Examples
--------
>>> m = curveball.models.fit_model(df)
>>> curveball.likelihood.plot_model_loglik(m, df)

"""
K = m.best_values['K']
y0 = m.best_values['y0']
r = m.best_values['r']
nu = m.best_values.get('nu', 1)
q0 = m.best_values.get('q0', np.inf)
v = m.best_values.get('v', np.inf)

rs = np.logspace(-2, 4, 100)
nus = np.logspace(-3, 2.5, 100)
q0s = np.logspace(-3, 2.5, 100)

ir = ( abs(rs - r).argmin() )
inu = ( abs(nus - nu).argmin() )
iq0 = ( abs(q0s - m.best_values.get('q0', np.inf)).argmin() )
ir_in = ( abs(rs - m.init_values['r']).argmin() )
inu_in = ( abs(nus - m.init_values.get('nu', 1)).argmin() )
iq0_in = ( abs(q0s - m.init_values.get('q0', np.inf)).argmin() )

L0 = loglik_r_nu(rs, nus, df, y0=y0, K=K, q0=q0, v=v)
L1 = loglik_r_q0(rs, q0s, df, y0=y0, K=K, nu=nu, v=v)

w, h= plt.rcParams['figure.figsize']
fig, ax = plt.subplots(1, 2, figsize=(w * 2, h))

plot_loglik(L0, rs, nus, normalize=False, fig_title=fig_title,
xlabel=r'$r$', ylabel=r'$\nu$', colorbar=False, ax=ax[0])
ax[0].scatter(ir, inu, marker='s', s=30, color='w', edgecolors='k')
ax[0].scatter(ir_in, inu_in, marker='^', s=30, color='w', edgecolors='k')

plot_loglik(L1, rs, q0s, normalize=False,
xlabel=r'$r$', ylabel=r'$q_0$', colorbar=True, ax=ax[1])
ax[1].scatter(ir, iq0, marker='s', s=30, color='w', edgecolors='k')
ax[1].scatter(ir_in, iq0_in, marker='^', s=30, color='w', edgecolors='k')

fig.tight_layout()
sns.despine()
return fig, ax