loss

loss#

class pygom.loss.Gamma(y, weights=None, shape=2.0)[source]#

Gamma distribution loss object

Parameters

y (array like) – observation
Shape (float) – shape (a in latex equations)

diff2Loss(yhat, apply_weighting=True)[source]#

Twice derivative of the loss function with respect to yhat. See:

Jupiter notebook “Loss function Calculations.ipnyb”

Parameters

yhat (array like) – observation
apply_weighting (boolean) – If True multiplies array of residuals by weightings, else raw residuals are used.

Returns

scnd_deriv_yhat – \(\frac{a \left(- \hat{y} + 2 y\right)}{\hat{y}^{3}}\)

Return type

array like

diff_loss(yhat, apply_weighting=True)[source]#

Derivative of the loss function with respect to yhat which is See:

File “Loss function Calculations.ipnyb”

Parameters

yhat (array like) – prediction
apply_weighting (boolean) – If True multiplies array of residuals by weightings, else raw residuals are used.

Returns

first_deriv_yhat – \(\frac{a \left(\hat{y} - y\right)}{\hat{y}^{2}}\)

Return type

array like

loss(yhat, apply_weighting=True)[source]#

The loss under a gamma distribution. Defined as the negative log-likelihood of the gamma distirbution in terms of mean and shape. See: Bolker, B. M. (2008). Gamma. In Ecological Models in R (pp. 131–133).

Princeton University Press. File “Loss function Calculations.ipnyb”

Parameters

yhat (array like) – prediction
apply_weighting (boolean) – Residuals are not used in calculation, as such True or False argument makes no difference. Argument has been kept as the equivalent calculation is made with residuals for other loss functions, so without it missing arguments errors could occur.

Return type

negative log-likelihood, \(\mathcal{L}(\hat{y}; y,a)\)

class pygom.loss.GammaLoss(theta, ode, x0, t0, t, y, state_name, state_weight=None, shape=2, target_param=None, target_state=None)[source]#: Realizations from a Gamma distribution taking parameters mean and shape.

class pygom.loss.NegBinom(y, weights=None, k=1.0)[source]#

Negative Binomial distribution loss object

Parameters

y (array like) – observation
k (float) – Overdispersion parameter (k=mean+mean(mean/variance))

diff2Loss(yhat, apply_weighting=True)[source]#

Twice derivative of the loss function with respect to yhat. See:

Jupiter notebook “Loss function Calculations.ipnyb” Bolker, B. M. (2008). Negative Binomial. In Ecological Models in R (pp. 124–126). Princeton University Press

Parameters

yhat (array like) – observation
apply_weighting (boolean) – Residuals are not used in calculation, as such True or False argument makes no difference. Argument has been kept as the equivalent calculation is made with residuals for other loss functions, so without it missing arguments errors could occur.
apply_weighting – If True multiplies array of residuals by weightings, else raw residuals are used.

Returns

scnd_deriv_yhat –

math

`frac{k(hat{y}(k + hat{y}) + hat{y}(y -hat{y}) + (k +: hat{y})(y - hat{y})}{hat{y}^{2}(k + hat{y})^{2}}`

Return type

array like

diff_loss(yhat, apply_weighting=True)[source]#

Derivative of the loss function with respect to yhat which is See:

Jupiter notebook “Loss function Calculations.ipnyb” Bolker, B. M. (2008). Negative Binomial. In Ecological Models in R (pp. 124–126). Princeton University Press

Parameters

yhat (array like) – observation apply_weighting: boolean Residuals are not used in calculation, as such True or False argument makes no difference. Argument has been kept as the equivalent calculation is made with residuals for other loss functions, so without it missing arguments errors could occur
apply_weighting (boolean) – If True multiplies array of residuals by weightings, else raw residuals are used.

Returns

first_deriv_yhat – \(\frac{k(\hat{y}-y)}{\hat{y}(k + \hat{y})}\)

Return type

array like

loss(yhat, apply_weighting=True)[source]#

The loss under a Negative Binomial distribution. Defined as the negative log-likelihood of the Negative Binomial 2 distribution.

Parameters

yhat (array like) – observation
apply_weighting (boolean) – Residuals are not used in calculation, as such True or False argument makes no difference. Argument has been kept as the equivalent calculation is made with residuals for other loss functions, so without it missing arguments errors could occur.

Return type

negative log-likelihood, \(\mathcal{L}(\hat{y}; y,k)\)

class pygom.loss.NegBinomLoss(theta, ode, x0, t0, t, y, state_name, state_weight=None, k=1, target_param=None, target_state=None)[source]#: Realizations from a Negative Binomial distribution

class pygom.loss.Normal(y, weights=None, sigma=1.0)[source]#

Normal distribution loss object

Parameters

y (array like) – observation
sigma (float) – standard deviation

diff2Loss(yhat, apply_weighting=True)[source]#

Twice derivative of the normal loss.

Parameters

yhat (array like) – observations
apply_weighting (boolean) – Residuals are not used in calculation, as such True or False argument makes no difference. Argument has been kept as the equivalent calculation is made with residuals for other loss functions, so without it missing arguments errors could occur.

Returns

s – inverse of the variance with shape = yhat.shape

Return type

array like

diff_loss(yhat, apply_weighting=True)[source]#

Derivative of the loss function which is \(\sigma^{-1}(y - \hat{y})\)

Parameters

yhat (array like) – observation
apply_weighting (boolean) – If True multiplies array of residuals by weightings, else raw residuals are used.

Returns

r – \(\nabla \mathcal{L}(\hat{y}, y)\)

Return type

array like

loss(yhat, apply_weighting=True)[source]#

The loss under a normal distribution. Defined as the negative log-likelihood here.

Parameters

yhat (array like) – observation
apply_weighting (boolean) – If True multiplies array of residuals by weightings, else raw residuals are used.

Returns

negative log-likelihood, – mathcalfrac{1}{sqrt{2pi}sigma}e^{-frac{(y-hat{y})^{2}}{2sigma^{2}}

Return type

math:

class pygom.loss.NormalLoss(theta, ode, x0, t0, t, y, state_name, state_weight=None, sigma=1.0, target_param=None, target_state=None)[source]#: Realizations from a Normal distribution

class pygom.loss.Poisson(y, weights=None)[source]#

Poisson distribution loss object

Parameters: y (array like) – observation

diff2Loss(yhat, apply_weighting=True)[source]#

Twice derivative of the Poisson loss.

Parameters

yhat (array like) – observations
apply_weighting (boolean) – Residuals are not used in calculation, as such True or False argument makes no difference. Argument has been kept as the equivalent calculation is made with residuals for other loss functions, so without it missing arguments errors could occur.

Returns

s – \(\frac{y}{\hat{y}^{2}}\) with shape = yhat.shape

Return type

array like

diff_loss(yhat, apply_weighting=True)[source]#

Derivative of the loss function, \(1 - y\hat{y}^{-1}\)

Parameters

yhat (array like) – observation
apply_weighting (boolean) – If True multiplies array of residuals by weightings, else raw residuals are used.

Return type

\(\nabla \mathcal{L}(\hat{y},y)\)

loss(yhat, apply_weighting=True)[source]#

The loss under a Poisson distribution. Defined as the negative log-likelihood here.

Parameters

yhat (array like) – observation
apply_weighting (boolean) – Residuals are not used in calculation, as such True or False argument makes no difference. Argument has been kept as the equivalent calculation is made with residuals for other loss functions, so without it missing arguments errors could occur.

Return type

negative log-likelihood, \(\mathcal{L}(\hat{y}, y)\)

class pygom.loss.PoissonLoss(theta, ode, x0, t0, t, y, state_name, state_weight=None, target_param=None, target_state=None)[source]#: Realizations from a Poisson distribution

class pygom.loss.Square(y, weights=None)[source]#

Square loss object

Parameters: y (array like) – observations

diff2Loss(yhat, apply_weighting=True)[source]#

Twice derivative of the square loss. Which is simply 2.

Parameters

yhat (array like) – observations
apply_weighting (boolean) – Residuals are not used in calculation, as such True or False argument makes no difference. Argument has been kept as the equivalent calculation is made with residuals for other loss functions, so without it missing arguments errors could occur.

Returns

either a scalar, vector or matrix depending on the shape of of the input yhat

Return type

array with values of 2

diff_loss(yhat, apply_weighting=True)[source]#

Derivative under square loss. Assuming that we are solving the minimization problem i.e. our objective function is the negative of the log-likelihood

Parameters

yhat (array like) – observation
apply_weighting (boolean) – If True multiplies array of residuals by weightings, else raw residuals are used.

Return type

\(-2(y_{i} - \hat{y}_{i})\)

loss(yhat, apply_weighting=True)[source]#

Loss under square loss. Not really saying much here

Parameters

yhat (array like) – observation
apply_weighting (boolean) – If True multiplies array of residuals by weightings, else raw residuals are used.

Return type

\(\sum_{i=1}^{n} (y-\hat{y})^{2}\)

class pygom.loss.SquareLoss(theta, ode, x0, t0, t, y, state_name, state_weight=None, target_param=None, target_state=None)[source]#: The square loss function

class pygom.loss.UnivariateSpline(x, y, w=None, bbox=[None, None], k=3, s=None, ext=0, check_finite=False)[source]#

1-D smoothing spline fit to a given set of data points.

Fits a spline y = spl(x) of degree k to the provided x, y data. s specifies the number of knots by specifying a smoothing condition.

Parameters

x ((N,) array_like) – 1-D array of independent input data. Must be increasing; must be strictly increasing if s is 0.
y ((N,) array_like) – 1-D array of dependent input data, of the same length as x.
w ((N,) array_like, optional) – Weights for spline fitting. Must be positive. If w is None, weights are all 1. Default is None.
bbox ((2,) array_like, optional) – 2-sequence specifying the boundary of the approximation interval. If bbox is None, bbox=[x[0], x[-1]]. Default is None.
k (int, optional) – Degree of the smoothing spline. Must be 1 <= k <= 5. k = 3 is a cubic spline. Default is 3.
s (float or None, optional) –
Positive smoothing factor used to choose the number of knots. Number of knots will be increased until the smoothing condition is satisfied:
```
sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s
```
However, because of numerical issues, the actual condition is:
```
abs(sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) - s) < 0.001 * s
```
If s is None, s will be set as len(w) for a smoothing spline that uses all data points. If 0, spline will interpolate through all data points. This is equivalent to InterpolatedUnivariateSpline. Default is None. The user can use the s to control the tradeoff between closeness and smoothness of fit. Larger s means more smoothing while smaller values of s indicate less smoothing. Recommended values of s depend on the weights, w. If the weights represent the inverse of the standard-deviation of y, then a good s value should be found in the range (m-sqrt(2*m),m+sqrt(2*m)) where m is the number of datapoints in x, y, and w. This means s = len(w) should be a good value if 1/w[i] is an estimate of the standard deviation of y[i].
ext (int or str, optional) –
Controls the extrapolation mode for elements not in the interval defined by the knot sequence.
- if ext=0 or ‘extrapolate’, return the extrapolated value.
- if ext=1 or ‘zeros’, return 0
- if ext=2 or ‘raise’, raise a ValueError
- if ext=3 or ‘const’, return the boundary value.
Default is 0.
check_finite (bool, optional) – Whether to check that the input arrays contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination or non-sensical results) if the inputs do contain infinities or NaNs. Default is False.

See also

BivariateSpline: a base class for bivariate splines.
SmoothBivariateSpline: a smoothing bivariate spline through the given points
LSQBivariateSpline: a bivariate spline using weighted least-squares fitting
RectSphereBivariateSpline: a bivariate spline over a rectangular mesh on a sphere
SmoothSphereBivariateSpline: a smoothing bivariate spline in spherical coordinates
LSQSphereBivariateSpline: a bivariate spline in spherical coordinates using weighted least-squares fitting
RectBivariateSpline: a bivariate spline over a rectangular mesh
InterpolatedUnivariateSpline: a interpolating univariate spline for a given set of data points.
bisplrep: a function to find a bivariate B-spline representation of a surface
bisplev: a function to evaluate a bivariate B-spline and its derivatives
splrep: a function to find the B-spline representation of a 1-D curve
splev: a function to evaluate a B-spline or its derivatives
sproot: a function to find the roots of a cubic B-spline
splint: a function to evaluate the definite integral of a B-spline between two given points
spalde: a function to evaluate all derivatives of a B-spline

Notes

The number of data points must be larger than the spline degree k.

NaN handling: If the input arrays contain nan values, the result is not useful, since the underlying spline fitting routines cannot deal with nan. A workaround is to use zero weights for not-a-number data points:

>>> import numpy as np
>>> from scipy.interpolate import UnivariateSpline
>>> x, y = np.array([1, 2, 3, 4]), np.array([1, np.nan, 3, 4])
>>> w = np.isnan(y)
>>> y[w] = 0.
>>> spl = UnivariateSpline(x, y, w=~w)

Notice the need to replace a nan by a numerical value (precise value does not matter as long as the corresponding weight is zero.)

References

Based on algorithms described in [1]_, [2]_, [3]_, and 4:

1: P. Dierckx, “An algorithm for smoothing, differentiation and integration of experimental data using spline functions”, J.Comp.Appl.Maths 1 (1975) 165-184.
2: P. Dierckx, “A fast algorithm for smoothing data on a rectangular grid while using spline functions”, SIAM J.Numer.Anal. 19 (1982) 1286-1304.
3: P. Dierckx, “An improved algorithm for curve fitting with spline functions”, report tw54, Dept. Computer Science,K.U. Leuven, 1981.
4: P. Dierckx, “Curve and surface fitting with splines”, Monographs on Numerical Analysis, Oxford University Press, 1993.

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import UnivariateSpline
>>> rng = np.random.default_rng()
>>> x = np.linspace(-3, 3, 50)
>>> y = np.exp(-x**2) + 0.1 * rng.standard_normal(50)
>>> plt.plot(x, y, 'ro', ms=5)

Use the default value for the smoothing parameter:

>>> spl = UnivariateSpline(x, y)
>>> xs = np.linspace(-3, 3, 1000)
>>> plt.plot(xs, spl(xs), 'g', lw=3)

Manually change the amount of smoothing:

>>> spl.set_smoothing_factor(0.5)
>>> plt.plot(xs, spl(xs), 'b', lw=3)
>>> plt.show()

antiderivative(n=1)[source]#

Construct a new spline representing the antiderivative of this spline.

Parameters: n (int, optional) – Order of antiderivative to evaluate. Default: 1
Returns: spline – Spline of order k2=k+n representing the antiderivative of this spline.
Return type: UnivariateSpline

Notes

New in version 0.13.0.

See also

splantider, derivative

Examples

>>> import numpy as np
>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(0, np.pi/2, 70)
>>> y = 1 / np.sqrt(1 - 0.8*np.sin(x)**2)
>>> spl = UnivariateSpline(x, y, s=0)

The derivative is the inverse operation of the antiderivative, although some floating point error accumulates:

>>> spl(1.7), spl.antiderivative().derivative()(1.7)
(array(2.1565429877197317), array(2.1565429877201865))

Antiderivative can be used to evaluate definite integrals:

>>> ispl = spl.antiderivative()
>>> ispl(np.pi/2) - ispl(0)
2.2572053588768486

This is indeed an approximation to the complete elliptic integral \(K(m) = \int_0^{\pi/2} [1 - m\sin^2 x]^{-1/2} dx\):

>>> from scipy.special import ellipk
>>> ellipk(0.8)
2.2572053268208538

derivative(n=1)[source]#

Construct a new spline representing the derivative of this spline.

Parameters: n (int, optional) – Order of derivative to evaluate. Default: 1
Returns: spline – Spline of order k2=k-n representing the derivative of this spline.
Return type: UnivariateSpline

See also

splder, antiderivative

Notes

New in version 0.13.0.

Examples

This can be used for finding maxima of a curve:

>>> import numpy as np
>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(0, 10, 70)
>>> y = np.sin(x)
>>> spl = UnivariateSpline(x, y, k=4, s=0)

Now, differentiate the spline and find the zeros of the derivative. (NB: sproot only works for order 3 splines, so we fit an order 4 spline):

>>> spl.derivative().roots() / np.pi
array([ 0.50000001,  1.5       ,  2.49999998])

This agrees well with roots \(\pi/2 + n\pi\) of \(\cos(x) = \sin'(x)\).

derivatives(x)[source]#

Return all derivatives of the spline at the point x.

Parameters: x (float) – The point to evaluate the derivatives at.
Returns: der – Derivatives of the orders 0 to k.
Return type: ndarray, shape(k+1,)

Examples

>>> import numpy as np
>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(0, 3, 11)
>>> y = x**2
>>> spl = UnivariateSpline(x, y)
>>> spl.derivatives(1.5)
array([2.25, 3.0, 2.0, 0])

get_coeffs()[source]#: Return spline coefficients.

get_knots()[source]#

Return positions of interior knots of the spline.

Internally, the knot vector contains 2*k additional boundary knots.

get_residual()[source]#

Return weighted sum of squared residuals of the spline approximation.

This is equivalent to:

sum((w[i] * (y[i]-spl(x[i])))**2, axis=0)

integral(a, b)[source]#

Return definite integral of the spline between two given points.

Parameters

a (float) – Lower limit of integration.
b (float) – Upper limit of integration.

Returns

integral – The value of the definite integral of the spline between limits.

Return type

float

Examples

>>> import numpy as np
>>> from scipy.interpolate import UnivariateSpline
>>> x = np.linspace(0, 3, 11)
>>> y = x**2
>>> spl = UnivariateSpline(x, y)
>>> spl.integral(0, 3)
9.0

which agrees with \(\int x^2 dx = x^3 / 3\) between the limits of 0 and 3.

A caveat is that this routine assumes the spline to be zero outside of the data limits:

>>> spl.integral(-1, 4)
9.0
>>> spl.integral(-1, 0)
0.0

roots()[source]#

Return the zeros of the spline.

Notes

Restriction: only cubic splines are supported by FITPACK. For non-cubic splines, use PPoly.root (see below for an example).

Examples

For some data, this method may miss a root. This happens when one of the spline knots (which FITPACK places automatically) happens to coincide with the true root. A workaround is to convert to PPoly, which uses a different root-finding algorithm.

For example,

>>> x = [1.96, 1.97, 1.98, 1.99, 2.00, 2.01, 2.02, 2.03, 2.04, 2.05]
>>> y = [-6.365470e-03, -4.790580e-03, -3.204320e-03, -1.607270e-03,
...      4.440892e-16,  1.616930e-03,  3.243000e-03,  4.877670e-03,
...      6.520430e-03,  8.170770e-03]
>>> from scipy.interpolate import UnivariateSpline
>>> spl = UnivariateSpline(x, y, s=0)
>>> spl.roots()
array([], dtype=float64)

Converting to a PPoly object does find the roots at x=2:

>>> from scipy.interpolate import splrep, PPoly
>>> tck = splrep(x, y, s=0)
>>> ppoly = PPoly.from_spline(tck)
>>> ppoly.roots(extrapolate=False)
array([2.])

See also

sproot, PPoly.roots

set_smoothing_factor(s)[source]#

Continue spline computation with the given smoothing factor s and with the knots found at the last call.

This routine modifies the spline in place.

pygom.loss.asymptotic(obj, alpha=0.05, theta=None, lb=None, ub=None)[source]#

Finds the confidence interval at the \(\alpha\) level under the \(\mathcal{X}^{2}\) assumption for the likelihood

Parameters

obj (ode object) – an object initialized from BaseLoss
alpha (numeric, optional) – confidence level, \(0 < \alpha < 1\). Defaults to 0.05.
theta (array like, optional) – the MLE parameters. Defaults to None which then theta will be inferred from the input obj
lb (array like, optional) – expected lower bound
ub (array like, optional) – expected upper bound

Returns

l (array like) – lower confidence interval
u (array like) – upper confidence interval

pygom.loss.bootstrap(obj, alpha=0.05, theta=None, lb=None, ub=None, iteration=0, full_output=False)[source]#

Finds the confidence interval at the \(\alpha\) level via bootstrap

Parameters

obj (ode object) – an object initialized from BaseLoss
alpha (numeric, optional) – confidence level, \(0 < \alpha < 1\). Defaults to 0.05.
theta (array like, optional) – the MLE parameters
lb (array like, optional) – upper bound for the parameters
ub (array like, optional) – lower bound for the parameters
iteration (int, optional) – number of bootstrap samples, defaults to 0 which is interpreted as \(min(2n, 100)\) where \(n\) is the number of data points.
full_output (bool) – if the full set of estimates is required.

Returns

l (array like) – lower confidence interval
u (array like) – upper confidence interval

pygom.loss.cost_grad_interpolant(ode, spline_list, t, theta)[source]#

Returns the cost (sum of squared residuals) and the gradient between the first derivative of the interpolant and the function of the ode

Parameters

ode (DeterministicOde) – an ode object
spline_list (list) – list of scipy.interpolate.UnivariateSpline
t (array like) – time
theta (array list) – parameter value

Returns

cost (double) – sum of squared residuals
g – gradient of the squared residuals

pygom.loss.cost_interpolant(ode, spline_list, t, theta, vec=True, aggregate=True)[source]#

Returns the cost (sum of squared residuals) between the first derivative of the interpolant and the function of the ode

Parameters

ode (DeterministicOde) – an ode object
spline_list (list) – list of scipy.interpolate.UnivariateSpline
t (array like) – time
theta (array list) – paramter value
vec (bool, optional) – if the matrix should be flattened to be a vector
aggregate (bool, optional) – sum the vector/matrix

Returns

cost – sum of squared residuals

Return type

double

pygom.loss.cost_sample(ode, fxApprox, xApprox, t, theta, vec=True, aggregate=True)[source]#

Returns the cost (sum of squared residuals) between the first derivative of the interpolant and the function of the ode using samples at time points t.

Parameters

ode (DeterministicOde) – an ode object
fxApprox (list) – list of approximated values for the first derivative
xApprox (list) – list of approximated values for the states
t (array like) – time
theta (array list) – parameter value
vec (bool, optional) – if the matrix should be flattened to be a vector.
aggregate (bool/str, optional) – sum the vector/matrix. If this is equals to ‘int’ then the Simpsons rule is applied to the samples. Also changes the behaviour of vec, where True outputs a vector where the elements contain the values of the integrand on each of the dimensions of the ode. False returns the sum of this vector, a scalar.

Returns

r – the cost or the residuals if vec is True

Return type

array list

See also

residual_sample()

pygom.loss.geometric(obj, alpha=0.05, theta=None, method='jtj', geometry='o', full_output=False)[source]#

Finds the geometric confidence interval under profiling at the \(\alpha\) level the normal approximation

Parameters

obj (ode object) – an object initialized from BaseLoss
alpha (numeric) – confidence level, \(0 < \alpha < 1\)
theta (array like, optional) – the MLE parameters. When None given, it tries to estimate the optimal using methods provided by obj
method (string) – construction of the covariance matrix. jtj is the \(J^{\top}\) where \(J\) is the Jacobian of the ode. ‘hessian’ is the hessian of the ode while ‘fisher’ is the fisher information found by \(cov(\nabla_{\theta}\mathcal{L})\).
geometry (string) – the two types of geometry defined in [Moolgavkar1987]. c geometry uses the covariance at the maximum likelihood estimate \(\hat{\theta}\), while the ‘o’ geometry is the covariance defined at point \(\theta\).
full_output (bool, optional) – If True then both the l_path and u_path will be outputted, else only the point estimates of l and u

Returns

l (array like) – lower confidence interval
u (array like) – upper confidence interval
l_path (list) – path from \(\hat{\theta}\) to the lower \(1 - \alpha/2\) point for all parameters
u_path (list) – same as l_path but for the upper confidence interval

pygom.loss.get_init(y, t, ode, theta=None, full_output=False)[source]#

Get an initial guess of theta given the observations y and the corresponding time points t.

Parameters

y (:array like) – observed values
t (array like) – time
ode (DeterministicOde) – an ode object
theta (array like) – parameter value
full_output (bool, optional) – True if the optimization result should be returned. Defaults to False.

Returns

theta – a guess of the parameters

Return type

array like

pygom.loss.grad_sample(ode, fxApprox, xApprox, t, theta, vec=False, output_residual=False)[source]#

Returns the gradient of the objective value using the state values of the interpolant given samples at time points t. Note that the parameters taken here is chosen to be same as cost_sample() for convenience.

Parameters

ode (DeterministicOde) – an ode object
fxApprox (list) – list of approximated values for the first derivative
xApprox (list) – list of approximated values for the states
t (array like) – time
theta (array list) – parameter value
vec (bool, optional) – if the matrix should be flattened to be a vector
output_residual (bool, optional) – if True, then the residuals will be returned as an additional argument

Returns

g – gradient of the objective function

Return type

numpy.ndarray

See also

jac_sample()

pygom.loss.interpolate(solution, t, s=0)[source]#

Interpolate the solution of the ode given the time points and a suitable smoothing vector using univariate spline

Parameters

solution (numpy.ndarray) – f(t) of the ode with the rows correspond to time
t (array like) – time
s (smoothing scalar, optional) – greater or equal to zero

Returns

splineList – of scipy.interpolate.UnivariateSpline

Return type

list

pygom.loss.jac_sample(ode, fxApprox, xApprox, t, theta, vec=True)[source]#

Returns the Jacobian of the objective value using the state values of the interpolant given samples at time points t. Note that the parameters taken here is chosen to be same as cost_sample() for convenience.

Parameters

ode (DeterministicOde) – an ode object
fxApprox (list) – list of approximated values for the first derivative
xApprox (list) – list of approximated values for the states
t (array like) – time
theta (array list) – parameter value
vec (bool, optional) – if the matrix should be flattened to be a vector

Returns

r – the residuals

Return type

array list

See also

cost_sample()

pygom.loss.leastsq(func, x0, args=(), Dfun=None, full_output=False, col_deriv=False, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)[source]#

Minimize the sum of squares of a set of equations.

x = arg min(sum(func(y)**2,axis=0))
         y

Parameters

func (callable) – Should take at least one (possibly length N vector) argument and returns M floating point numbers. It must not return NaNs or fitting might fail. M must be greater than or equal to N.
x0 (ndarray) – The starting estimate for the minimization.
args (tuple, optional) – Any extra arguments to func are placed in this tuple.
Dfun (callable, optional) – A function or method to compute the Jacobian of func with derivatives across the rows. If this is None, the Jacobian will be estimated.
full_output (bool, optional) – If True, return all optional outputs (not just x and ier).
col_deriv (bool, optional) – If True, specify that the Jacobian function computes derivatives down the columns (faster, because there is no transpose operation).
ftol (float, optional) – Relative error desired in the sum of squares.
xtol (float, optional) – Relative error desired in the approximate solution.
gtol (float, optional) – Orthogonality desired between the function vector and the columns of the Jacobian.
maxfev (int, optional) – The maximum number of calls to the function. If Dfun is provided, then the default maxfev is 100*(N+1) where N is the number of elements in x0, otherwise the default maxfev is 200*(N+1).
epsfcn (float, optional) – A variable used in determining a suitable step length for the forward- difference approximation of the Jacobian (for Dfun=None). Normally the actual step length will be sqrt(epsfcn)*x If epsfcn is less than the machine precision, it is assumed that the relative errors are of the order of the machine precision.
factor (float, optional) – A parameter determining the initial step bound (factor * || diag * x||). Should be in interval (0.1, 100).
diag (sequence, optional) – N positive entries that serve as a scale factors for the variables.

Returns

x (ndarray) – The solution (or the result of the last iteration for an unsuccessful call).
cov_x (ndarray) – The inverse of the Hessian. fjac and ipvt are used to construct an estimate of the Hessian. A value of None indicates a singular matrix, which means the curvature in parameters x is numerically flat. To obtain the covariance matrix of the parameters x, cov_x must be multiplied by the variance of the residuals – see curve_fit. Only returned if full_output is True.
infodict (dict) – a dictionary of optional outputs with the keys:

nfev
The number of function calls

fvec
The function evaluated at the output

fjac
A permutation of the R matrix of a QR factorization of the final approximate Jacobian matrix, stored column wise. Together with ipvt, the covariance of the estimate can be approximated.

ipvt
An integer array of length N which defines a permutation matrix, p, such that fjac*p = q*r, where r is upper triangular with diagonal elements of nonincreasing magnitude. Column j of p is column ipvt(j) of the identity matrix.

qtf
The vector (transpose(q) * fvec).

Only returned if full_output is True.
mesg (str) – A string message giving information about the cause of failure. Only returned if full_output is True.
ier (int) – An integer flag. If it is equal to 1, 2, 3 or 4, the solution was found. Otherwise, the solution was not found. In either case, the optional output variable ‘mesg’ gives more information.

See also

least_squares: Newer interface to solve nonlinear least-squares problems with bounds on the variables. See method='lm' in particular.

Notes

“leastsq” is a wrapper around MINPACK’s lmdif and lmder algorithms.

cov_x is a Jacobian approximation to the Hessian of the least squares objective function. This approximation assumes that the objective function is based on the difference between some observed target data (ydata) and a (non-linear) function of the parameters f(xdata, params)

func(params) = ydata - f(xdata, params)

so that the objective function is

  min   sum((ydata - f(xdata, params))**2, axis=0)
params

The solution, x, is always a 1-D array, regardless of the shape of x0, or whether x0 is a scalar.

Examples

>>> from scipy.optimize import leastsq
>>> def func(x):
...     return 2*(x-3)**2+1
>>> leastsq(func, 0)
(array([2.99999999]), 1)

pygom.loss.minimize_scalar(fun, bracket=None, bounds=None, args=(), method=None, tol=None, options=None)[source]#

Local minimization of scalar function of one variable.

Parameters

fun (callable) –
Objective function. Scalar function, must return a scalar.

Suppose the callable has signature f0(x, *my_args, **my_kwargs), where my_args and my_kwargs are required positional and keyword arguments. Rather than passing f0 as the callable, wrap it to accept only x; e.g., pass fun=lambda x: f0(x, *my_args, **my_kwargs) as the callable, where my_args (tuple) and my_kwargs (dict) have been gathered before invoking this function.
bracket (sequence, optional) – For methods ‘brent’ and ‘golden’, bracket defines the bracketing interval and is required. Either a triple (xa, xb, xc) satisfying xa < xb < xc and func(xb) < func(xa) and func(xb) < func(xc), or a pair (xa, xb) to be used as initial points for a downhill bracket search (see scipy.optimize.bracket). The minimizer res.x will not necessarily satisfy xa <= res.x <= xb.
bounds (sequence, optional) – For method ‘bounded’, bounds is mandatory and must have two finite items corresponding to the optimization bounds.
args (tuple, optional) – Extra arguments passed to the objective function.
method (str or callable, optional) –
Type of solver. Should be one of:
- Brent
- Bounded
- Golden
- custom - a callable object (added in version 0.14.0), see below
Default is “Bounded” if bounds are provided and “Brent” otherwise. See the ‘Notes’ section for details of each solver.
tol (float, optional) – Tolerance for termination. For detailed control, use solver-specific options.
options (dict, optional) –
A dictionary of solver options.

maxiterint
Maximum number of iterations to perform.

dispbool
Set to True to print convergence messages.

See show_options() for solver-specific options.

Returns

res – The optimization result represented as a OptimizeResult object. Important attributes are: x the solution array, success a Boolean flag indicating if the optimizer exited successfully and message which describes the cause of the termination. See OptimizeResult for a description of other attributes.

Return type

OptimizeResult

See also

minimize: Interface to minimization algorithms for scalar multivariate functions
show_options: Additional options accepted by the solvers

Notes

This section describes the available solvers that can be selected by the ‘method’ parameter. The default method is the "Bounded" Brent method if bounds are passed and unbounded "Brent" otherwise.

Method Brent uses Brent’s algorithm [1]_ to find a local minimum. The algorithm uses inverse parabolic interpolation when possible to speed up convergence of the golden section method.

Method Golden uses the golden section search technique [1]_. It uses analog of the bisection method to decrease the bracketed interval. It is usually preferable to use the Brent method.

Method Bounded can perform bounded minimization [2]_ [3]_. It uses the Brent method to find a local minimum in the interval x1 < xopt < x2.

Note that the Brent and Golden methods do not guarantee success unless a valid bracket triple is provided. If a three-point bracket cannot be found, consider scipy.optimize.minimize. Also, all methods are intended only for local minimization. When the function of interest has more than one local minimum, consider global_optimization.

Custom minimizers

It may be useful to pass a custom minimization method, for example when using some library frontend to minimize_scalar. You can simply pass a callable as the method parameter.

The callable is called as method(fun, args, **kwargs, **options) where kwargs corresponds to any other parameters passed to minimize (such as bracket, tol, etc.), except the options dict, which has its contents also passed as method parameters pair by pair. The method shall return an OptimizeResult object.

The provided method callable must be able to accept (and possibly ignore) arbitrary parameters; the set of parameters accepted by minimize may expand in future versions and then these parameters will be passed to the method. You can find an example in the scipy.optimize tutorial.

New in version 0.11.0.

References

1: Press, W., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes in C. Cambridge University Press.
2: Forsythe, G.E., M. A. Malcolm, and C. B. Moler. “Computer Methods for Mathematical Computations.” Prentice-Hall Series in Automatic Computation 259 (1977).
3: Brent, Richard P. Algorithms for Minimization Without Derivatives. Courier Corporation, 2013.

Examples

Consider the problem of minimizing the following function.

>>> def f(x):
...     return (x - 2) * x * (x + 2)**2

Using the Brent method, we find the local minimum as:

>>> from scipy.optimize import minimize_scalar
>>> res = minimize_scalar(f)
>>> res.fun
-9.9149495908

The minimizer is:

>>> res.x
1.28077640403

Using the Bounded method, we find a local minimum with specified bounds as:

>>> res = minimize_scalar(f, bounds=(-3, -1), method='bounded')
>>> res.fun  # minimum
3.28365179850e-13
>>> res.x  # minimizer
-2.0000002026

class pygom.loss.partial[source]#

partial(func, *args, **keywords) - new function with partial application of the given arguments and keywords.

args#: tuple of arguments to future partial calls

func#: function object to use in future partial calls

keywords#: dictionary of keyword arguments to future partial calls

pygom.loss.profile(obj, alpha, theta=None, lb=None, ub=None, full_output=False)[source]#

Finds the profile confidence interval at the \(\alpha\) level under the \(\mathcal{X}^{2}\) assumption for the likelihood

Parameters

obj (ode object) – an object initialized from BaseLoss
alpha (numeric) – confidence level, \(0 < \alpha < 1\)
theta (array like, optional) – the MLE parameters. When None given, it tries to estimate the optimal using methods provided by obj
lb (array like, optional) – expected lower bound
ub (array like, optional) – expected upper bound
full_output (bool, optional) – if more output is desired

Returns

l (array like) – lower confidence interval
u (array like) – upper confidence interval

pygom.loss.residual_interpolant(ode, spline_list, t, theta, vec=True)[source]#

Returns the residuals between the first derivative of the interpolant and the function of the ode

Parameters

ode (DeterministicOde) – an ode object
spline_list (list) – list of scipy.interpolate.UnivariateSpline
t (array like) – time
theta (array list) – parameter value
vec (bool, optional) – if the matrix should be flattened to be a vector
aggregate (bool, optional) – sum the vector/matrix

Returns

r – the residuals

Return type

array list

pygom.loss.residual_sample(ode, fxApprox, xApprox, t, theta, vec=True)[source]#

Returns the residuals between the first derivative of the interpolant and the function of the ode using samples at time points t.

Parameters

ode (DeterministicOde) – an ode object
fxApprox (list) – list of approximated values for the first derivative
xApprox (list) – list of approximated values for the states
t (array like) – time
theta (array list) – parameter value
vec (bool, optional) – if the matrix should be flattened to be a vector

Returns

r – the residuals

Return type

array list

See also

cost_sample()

pygom.loss.simps(y, x=None, *, dx=1.0, axis=-1)#

Integrate y(x) using samples along the given axis and the composite Simpson’s rule. If x is None, spacing of dx is assumed.

Parameters

y (array_like) – Array to be integrated.
x (array_like, optional) – If given, the points at which y is sampled.
dx (float, optional) – Spacing of integration points along axis of x. Only used when x is None. Default is 1.
axis (int, optional) – Axis along which to integrate. Default is the last axis.

Returns

The estimated integral computed with the composite Simpson’s rule.

Return type

float

See also

quad: adaptive quadrature using QUADPACK
fixed_quad: fixed-order Gaussian quadrature
dblquad: double integrals
tplquad: triple integrals
romb: integrators for sampled data
cumulative_trapezoid: cumulative integration for sampled data
cumulative_simpson: cumulative integration using Simpson’s 1/3 rule

Notes

For an odd number of samples that are equally spaced the result is exact if the function is a polynomial of order 3 or less. If the samples are not equally spaced, then the result is exact only if the function is a polynomial of order 2 or less.

References

1: Cartwright, Kenneth V. Simpson’s Rule Cumulative Integration with MS Excel and Irregularly-spaced Data. Journal of Mathematical Sciences and Mathematics Education. 12 (2): 1-9

Examples

>>> from scipy import integrate
>>> import numpy as np
>>> x = np.arange(0, 10)
>>> y = np.arange(0, 10)

>>> integrate.simpson(y, x=x)
40.5

>>> y = np.power(x, 3)
>>> integrate.simpson(y, x=x)
1640.5
>>> integrate.quad(lambda x: x**3, 0, 9)[0]
1640.25

loss

Contents

loss#