DCA - API

class dca.decline_curve_analysis.Arps(theta1: float, theta2: float, theta3: float)

Create an Arps instance.

Parameters:

theta1 (float) – Parameter number one. Interpreted as the logarithm of the integral from zero to infinity. This scales the generalized Pareto distribution.
theta2 (float) – Parameter number two.
theta3 (float) – Parameter number three.

Notes

We use the parametrization given in the paper ‘Bayesian Hierarchical Modeling’ by Se Yoon Lee et al by default. This parametrization is over unbounded variables \(\boldsymbol{\theta} = (\theta_1, \theta_2, \theta_3)\).

\[\mu(t; \boldsymbol{\theta}) = \exp \left( \theta_1 - \theta_2 - (1 + \exp(-\theta_3)) \log(1 + t \exp(\theta_3 - \theta_2)) \right)\]

It’s also possible to use the “original” parametrization with parameters \((q_1, h, D)\)

See also:: https://en.wikipedia.org/wiki/Generalized_Pareto_distribution

Examples

>>> model = Arps(1, -2, 1)
>>> time = np.arange(5)
>>> model(time)
array([20.08...,  0.31...,  0.12...,  0.072... ,  0.04...])

>>> model = Arps.from_original_parametrization(q_1=10, h=0.66, D=2)

__init__(theta1: float, theta2: float, theta3: float)

add(k)

Add a constant k in log space and return a new instance.

Examples

>>> model = Arps(10, 3, 1)
>>> model.add(2)
Arps(12.0, 3.0, 1.0)

eval(t): Evaluate \(\mu(t; \boldsymbol{\theta})\).

eval_grad_log(t)

Evaluate the gradient of the log-Arps curve.

Returns:: Array of shape (3, len(t)) with gradient information. To compute the gradient with respect to the parameters \(\boldsymbol{\theta}\), the caller must sum over \(t\).
Return type:: np.ndarray

Notes

This method evaluates the gradient of the log-Arps curve:

\[\nabla_{\boldsymbol{\theta}} \log \mu(t; \boldsymbol{\theta}) = \frac{1}{\mu(t; \boldsymbol{\theta})} \nabla_{\boldsymbol{\theta}} \mu(t; \boldsymbol{\theta})\]

Note that this method can be used to evaluate the gradient of the Arps curve, since:

\[\nabla_{\boldsymbol{\theta}} \mu(t; \boldsymbol{\theta}) = \mu(t; \boldsymbol{\theta}) \nabla_{\boldsymbol{\theta}} \log \mu(t; \boldsymbol{\theta})\]

Examples

>>> model = Arps(1, -2, 1)
>>> time = np.arange(5)
>>> gradient = model.eval_grad_log(time)
>>> gradient.shape
(3, 5)
>>> gradient.sum(axis=1)
array([5.        , 0.334..., 0.286...])

eval_integral(t): Evaluate the anti-derivative of \(\mu(t; \boldsymbol{\theta})\).

eval_integral_from_0_to_T(T): Evaluate \(\int_0^{T} \mu(t; \boldsymbol{\theta}) dt\).

eval_integral_from_0_to_inf(): Evaluate \(\int_0^{\infty} \mu(t; \boldsymbol{\theta}) dt\).

eval_integral_from_T1_to_T2(T1, T2): Evaluate \(\int_{T_1}^{T_2} \mu(t; \boldsymbol{\theta}) dt\).

eval_integral_from_T_to_inf(T): Evaluate \(\int_T^{\infty} \mu(t; \boldsymbol{\theta}) dt\).

eval_log(t): Evaluate \(\log \mu(t; \boldsymbol{\theta})\).

static eval_log_original_parametrization(t, q_1, h, D): Evaluate \(\log \mu(t; q_1, h, D)\).

static eval_original_parametrization(t, q_1, h, D): Evaluate \(\mu(t; q_1, h, D)\).

eval_squared_integral(t): Evaluate the anti-derivative of \(\mu(t; \boldsymbol{\theta})^2\).

classmethod from_original_parametrization(q_1, h, D): Initialize an instance using the original parametrization.

original_parametrization()

Return parameters \((q_1, h, D)\) in the original parametrization of the Arps curve, given by:

\[\mu(t; q_1, h, D) = \frac{q_1}{(1 + hDt)^{1/h}}\]

Examples

>>> Arps(1, -2, 1).original_parametrization()
(20.08..., 0.73..., 27.47...)

shift_and_scale(shift=0, scale=1)

Scale and shift the curve along the horizontal axis.

Equivalent to fitting on a transformed grid:: \(t' = t \times \text{scale} + \text{shift}\).

Examples

>>> curve = Arps(10, 5, 1)
>>> x_grid = np.array([0, 1, 2, 3, 5])
>>> curve.eval_log(x_grid)
array([5.        , 4.97517309, 4.95078876, 4.92683151, 4.8801402 ])
>>> curve.shift_and_scale(scale=2, shift=3).eval_log((x_grid * 2 + 3))
array([5.        , 4.97517309, 4.95078876, 4.92683151, 4.8801402 ])

property thetas: Return the thetas as a tuple.

class dca.decline_curve_analysis.Exponential(theta1: float, theta2: float)

Create an Exponential instance.

Parameters:

theta1 (float) – Parameter number one. Interpreted as the logarithm of the integral from zero to infinity. This scales the exponential distribution.
theta2 (float) – Parameter number two. This is the rate parameter. Compared to Wikipedia’s article on the exponential distribution, lambda = exp(theta2).

Notes

The Exponential class holds an exponential curve, \(\mu(t; \boldsymbol{\theta})\) given by,

\[\log \left( \mu(t; \theta_1, \theta_2) \right) = \theta_1 - \theta_2 - t \exp(-\theta_2),\]

which assures that the curve is positive and decreasing.

See also:: https://en.wikipedia.org/wiki/Exponential_distribution

Examples

>>> model = Exponential(1, -1)
>>> time = np.arange(5)
>>> model(time)
array([7.38905610e+00, 4.87589299e-01, 3.21750601e-02, 2.12316902e-03,
       1.40103753e-04])

Verify that Exponential(theta1, theta2) == Arps(theta1, theta2, -∞):

>>> import numpy as np
>>> exponential = Exponential(1, -1.3)
>>> arps = Arps(1, -1.3, -20)  # Small value of theta3
>>> np.allclose(exponential(time), arps(time))
True

This example shows the difference between midpoint evaluation and integrating the exponential decline curve within each time period. Since the decline curve is convex, using the midpoint rule underestimates, but very slightly:

>>> exponential = Exponential.from_original_parametrization(C=1, k=0.1)
>>> time_on = np.array([1, 0.5, 1, 0.4, 1])
>>> midpoints = np.cumsum(time_on) - time_on / 2
>>> exponential.eval(midpoints) * time_on
array([0.95122942, 0.44124845, 0.81873075, 0.3053518 , 0.71177032])
>>> T2 = np.cumsum(time_on)
>>> T1 = T2 - time_on
>>> exponential.eval_integral_from_T1_to_T2(T1, T2)
array([0.95162582, 0.44129442, 0.81907193, 0.30537215, 0.71206693])

__init__(theta1: float, theta2: float)

add(k)

Add a constant k in log space and return a new instance.

Examples

>>> model = Exponential(10, 3)
>>> model.add(2)
Exponential(12.0, 3.0)

eval(t): Evaluate \(\mu(t; \boldsymbol{\theta})\).

eval_grad_log(t)

Evaluate the gradient of the log-Exponential curve.

Returns:: Array of shape (2, len(t)) with gradient information. To compute the gradient with respect to the parameters \(\boldsymbol{\theta}\), the caller must sum over \(t\).
Return type:: np.array

Notes

This method evaluates the gradient of the log-Exponential curve:

\[\nabla_{\boldsymbol{\theta}} \log \mu(t; \boldsymbol{\theta}) = \frac{1}{\mu(t; \boldsymbol{\theta})} \nabla_{\boldsymbol{\theta}} \mu(t; \boldsymbol{\theta})\]

Note that this method can be used to evaluate the gradient of the Exponential curve, since:

\[\nabla_{\boldsymbol{\theta}} \mu(t; \boldsymbol{\theta}) = \mu(t; \boldsymbol{\theta}) \nabla_{\boldsymbol{\theta}} \log \mu(t; \boldsymbol{\theta})\]

Examples

>>> model = Exponential(5, 0.66)
>>> time = np.arange(5)
>>> gradient = model.eval_grad_log(time)
>>> gradient.shape
(2, 5)
>>> gradient.sum(axis=1)
array([5.        , 0.16851334])

eval_integral(t): Evaluate the anti-derivative of \(\mu(t; \boldsymbol{\theta})\).

eval_integral_from_0_to_T(T): Evaluate \(\int_0^{T} \mu(t; \boldsymbol{\theta}) dt\).

eval_integral_from_0_to_inf(): Evaluate \(\int_0^{\infty} \mu(t; \boldsymbol{\theta}) dt\).

eval_integral_from_T1_to_T2(T1, T2): Evaluate \(\int_{T_1}^{T_2} \mu(t; \boldsymbol{\theta}) dt\).

eval_integral_from_T_to_inf(T): Evaluate \(\int_T^{\infty} \mu(t; \boldsymbol{\theta}) dt\).

eval_log(t): Evaluate \(\log \mu(t; \boldsymbol{\theta})\).

eval_squared_integral(t): Evaluate the anti-derivative of \(\mu(t; \boldsymbol{\theta})^2\).

classmethod from_original_parametrization(C, k)

Initialize an instance using the original parametrization.

Examples

>>> C, k = 10, 0.1
>>> theta1, theta2 = Exponential.from_original_parametrization(C, k).thetas
>>> Exponential(theta1, theta2).original_parametrization()
(10.0000000..., 0.1000000...)

original_parametrization()

Return parameters \((C, k)\) in the original parametrization of the Exponential curve, given by:

\[\mu(t; C, k) = C \exp(-k t)\]

Examples

>>> time = np.linspace(0, 10)
>>> C, k = 10, 0.1
>>> curve = Exponential.from_original_parametrization(C, k)
>>> np.allclose(C * np.exp(-k * time), curve(time))
True

shift_and_scale(shift=0, scale=1)

Scale and shift the curve along the horizontal axis.

Equivalent to fitting on a transformed grid:: \(t' = t \times \text{scale} + \text{shift}\).

Examples

>>> curve = Exponential(10, 1)
>>> x_grid = np.array([0, 1, 2, 3, 5])
>>> curve.eval_log(x_grid)
array([9.        , 8.63212056, 8.26424112, 7.89636168, 7.16060279])
>>> curve.shift_and_scale(scale=2, shift=3).eval_log((x_grid * 2 + 3))
array([9.        , 8.63212056, 8.26424112, 7.89636168, 7.16060279])

property thetas: Return the thetas as a tuple.

class dca.decline_curve_analysis.CurveLoss(*, curve_func, mu=None, precision_matrix=None, p=1.4, half_life=None, curve_func_grad=None)

Create a curve loss instance.

Parameters:

curve_func (callable) – A callable with signature (x, theta1, theta2, theta3), returning y. The thetas will be optimized by curve fitting.
mu (np.ndarray, optional) – A vector of prior means for curve fitting parameters. The default is None (interpreted as no prior).
precision_matrix (np.ndarray, optional) – A matrix of prior precisions for curve fitting parameters. (The inverse of the covariance matrix). The default is None (interpreted as no prior).
p (float, optional) – The p in the p-norm used in the loss function. The default is 1.4.
half_life (float, optional) – Half life used to create weights for exponential decay, weighting the most recent data points more. For instance, if half_life=100, then a data point today is worth twice as much as a data point from 100 days ago, and four times as much as a data point 200 days ago. The default is None (interpreted as infinity, or no decay).
curve_func_grad (callable) – A callable with signature (x, theta1, theta2, theta3), which returns a numpy array of shape (3, len(x)) representing the gradient at each point in x.

Notes

Loss function used for curve fitting.

We use a class to set up common hyperparmeters.

Curves are fit using a multivariate Gaussian prior. An instance of this class is meant to be used with scipy.optimize.minimize.

The loss function is:

\[\mathcal{L}(\boldsymbol{\theta}) = \sum_i w_i (\gamma) \lvert \ln y_i - \ln f(t_i; \boldsymbol{\theta}) \rvert ^p + \alpha \left( \boldsymbol{\theta} - \boldsymbol{\mu}_{\boldsymbol{\theta}} \right)^T \boldsymbol{\Sigma}_{\boldsymbol{\theta}}^{-1} \left( \boldsymbol{\theta} - \boldsymbol{\mu}_{\boldsymbol{\theta}} \right)\]

where,

\(w_i(\gamma)\) is a weighting function that scales up the contribution of the most recent data,

the choice of whether we use a logarithm or not is determined by cross-validation,

\(1 \leq p \leq 2\) is a hyperparameter determined by cross-validation,

a multivariate normal prior is put on \(\boldsymbol{\theta}\) (the second term with the quadratic function)

the strength of the prior is given by \(\alpha\) and determined by cross-validation.

Regularization of the prior is obatined by scaling the precision_matrix with a scalar value.

Examples

>>> import numpy as np
>>> from scipy.optimize import minimize
>>> mu = np.array([14.615444, 1.879282, 0.004616])
>>> precision_matrix = np.array([[ 0.80, -0.32, -0.11],
...                              [-0.32,  0.67, -0.19],
...                              [-0.11, -0.19,  0.26]])
>>> def log_Arps(t, theta1, theta2, theta3):
...     return Arps(theta1, theta2, theta3).eval_log(t)
>>> loss_function = CurveLoss(mu=mu, precision_matrix=precision_matrix,
...                           curve_func=log_Arps)
>>> x = np.arange(5)
>>> y = np.array([10, 8, 6, 4, 1])
>>> optimization_result = minimize(loss_function, x0=np.array([1, 1, 1]),
...                                method="BFGS", args=(x, y,),)
>>> np.allclose(optimization_result.x,
...             np.array([14.38081254, 4.30188387, 7.62215325]))
True

Return type:: None.

__call__(thetas, x, y, weights=None)

Get the loss of the curve.

Parameters:

thetas (np.ndarray) – Parameters to be optimized.
x (np.ndarray) – The x-values (typically time).
y (np.ndarray) – The y-values.
weights (np.ndarray) – Weights used for each data point. Will be multiplied with existing weights given by the half life parameter, if they exist.

Returns:

Negative log posterior, proportional to log(likelihood) + log(prior). Minimizing this is the same as maximizing the likelihood.

Return type:

float

__init__(*, curve_func, mu=None, precision_matrix=None, p=1.4, half_life=None, curve_func_grad=None)

grad(thetas, x, y, weights=None)

Gradient of the loss function.

Returns an array of shape (num_parameters,). For instance, if the DCA curve is the Arps curve, then a length-3 array is returned, which represents the objective function differentiated with respect to theta_1, theta_2 and theta_3.

class dca.models.AR1Model(curve_cls)

Class for an autoregressive regression model, consisting of:

A non linear-regression and
errors that follow an AR(1) model

We need two or three things to make a statistical model like this work:

Simulation
Likelihood
Gradients (not needed, but makes optimizer hyperparams faster/robust)

In probabilistic programming languages like Stan, you only need to specify the forward simulation model. Stan will infer the likelihood and the gradients automatically. Stan also has the option of running MCMC instead of optimizing. Using Stan instead of writing out models “by hand” could be a good option, and samples from the posterior via MCMC could be useful. We chose not to use Stan because of difficulties with C++ compiles across platforms, setting up environments on Windows, etc. Maybe in the future!

An important thing to know about this model is that it underspecified without a prior on the autocorrelation coefficient phi. A declining trend can either be interpreted as:

an actual declining trend, with some autocorrlated noise on top
an arbitrary trend, with huge autocorrelation

In case (2), the whole trend is interpreted as one big autocorrelated noise event. This happens as phi goes to one (typically phi>0.99), and it’s not what we want. It’s important to regularize phi with an appropriate prior.

References

Examples

>>> from dca import Arps
>>> model = AR1Model(curve_cls=Arps)
>>> model = model.update(theta=(5, 3, 1), sigma=1, phi=0.5, p=2)
>>> tau = np.ones(10)
>>> t = np.linspace(0, 10, num=10)
>>> simulations = model.simulate(tau=tau, t=t, seed=42, simulations=5)
>>> simulations.shape
(5, 10)

A good and a bad guess:

>>> float(model.negative_ll(theta=(5, 3, 1), sigma=1, phi=0.5, p=2, tau=tau,
...                   t=t, y=simulations.mean(axis=0)).sum())
7.35...
>>> float(model.negative_ll(theta=(0, 3, 1), sigma=1, phi=0.5, p=2, tau=tau,
...                   t=t, y=simulations.mean(axis=0)).sum())
89.6...

__init__(curve_cls)

static get_weights(*, t, tau, half_life=inf)

Return weights for exponential decay going back in time. If half_life is set to infinity, then all weights become unity.

Examples

>>> tau = np.ones(4)  # Assuming tau is constant for simplicity
>>> t = np.cumsum(tau) - tau / 2
>>> AR1Model.get_weights(t=t, tau=tau)
array([1., 1., 1., 1.])
>>> AR1Model.get_weights(t=t, tau=tau, half_life=2)
array([0.552..., 0.781..., 1.104... , 1.562...])

gradient(theta, sigma, phi, p, *, tau, t, y, half_life=inf, prev_eta=None)

Returns the gradient of the negative-log-likelihood with respect to:

theta
sigma
phi
p

Returns an array of shape (len(theta) + 3, len(t)).

negative_ll(theta, sigma, phi, p, *, tau, t, y, half_life=inf, prev_eta=None)

Returns the elementwise negative log-likelihood of the model.

Likelihood of \(\log(y_0)\) if prev_eta=None:

\[\begin{split}P(\log y_0) &\sim \mathcal{N}(\text{loc}, \text{scale}) \\ \quad \text{loc} &= \log(f(t_i; \theta)) \\ \quad \text{scale} &= \frac{\text{std}(\text{GN}(\text{std}=\sigma, \beta=p))}{\sqrt{1 - \phi^2}}\end{split}\]

Likelihood of \(\log(y_i)\) given \(\log(y_{i-1})\):

\[\begin{split}P(\log y_i \mid \log y_{i-1}) &\sim \text{GN}(\text{loc}, \text{scale}, \beta=p) \\ \quad \text{loc} &= \log(f(t_i; \theta)) + \phi \left[ \log(y_{i-1}) - \log(f(t_{i-1}; \theta)) \right] \\ \quad \text{scale} &= \sqrt{\tau_i} \sigma\end{split}\]

Returns a vector with entries \(-w_i \log P(\log y_i \mid \log y_{i-1})\).

optimize(*, tau, t, y, prior_strength=0.0, half_life=inf, neg_logpdf_prior=None): Optimize over parameters (theta, sigma, phi, p).

pointwise_std(*, p, sigma, phi)

Compute the pointwise standard deviation.

\[\frac{\sqrt{\Gamma(3 / \beta)}}{\sqrt{\Gamma(1 / \beta)}} \frac{\sigma}{\sqrt{1 - \phi^2}}\]

See:

simulate(*, tau, t, seed=None, prev_eta=None, simulations=1)

Simulate production rates \(y\) from the model:

\[\begin{split}\log(y_i) &= \log(f(t_i)) + \eta_i \\ \quad \eta_i &= \phi \cdot \eta_{i-1} + \sqrt{\tau_i} \cdot \epsilon_i \\ \quad \epsilon_i &\sim \text{GN}(\mu=0, \alpha=\sigma, \beta)\end{split}\]

where:

\(y_i = \frac{p_i}{\tau_i} = \text{production per time}\)
\(\tau_i = \text{time\_on}\)
GN is the Generalized Normal with p-norm given by \(\beta\)

To get simulated production, multiply by \(\text{time\_on}\).

The returned numpy array has shape (simulations, len(x)).

update(**kwargs): Update parameters (theta, sigma, phi, p) to values on the original, bounded scale. Setting them to None means they will be inferred by optimization.