"""Frank module."""
import sys
import numpy as np
import scipy.integrate as integrate
from scipy.optimize import least_squares
from copulas import EPSILON
from copulas.bivariate.base import Bivariate, CopulaTypes
from copulas.bivariate.utils import split_matrix
MIN_FLOAT_LOG = np.log(sys.float_info.min)
MAX_FLOAT_LOG = np.log(sys.float_info.max)
[docs]class Frank(Bivariate):
"""Class for Frank copula model."""
copula_type = CopulaTypes.FRANK
theta_interval = [-float('inf'), float('inf')]
invalid_thetas = [0]
[docs] def generator(self, t):
"""Return the generator function."""
a = (np.exp(-self.theta * t) - 1) / (np.exp(-self.theta) - 1)
return -np.log(a)
def _g(self, z):
r"""Assist in solving the Frank copula.
This functions encapsulates :math:`g(z) = e^{-\theta z} - 1` used on Frank copulas.
Argument:
z: np.ndarray
Returns:
np.ndarray
"""
return np.exp(-self.theta * z) - 1
[docs] def probability_density(self, X):
r"""Compute probability density function for given copula family.
The probability density(PDF) for the Frank family of copulas correspond to the formula:
.. math:: c(U,V) = \frac{\partial^2 C(u,v)}{\partial v \partial u} =
\frac{-\theta g(1)(1 + g(u + v))}{(g(u) g(v) + g(1)) ^ 2}
Where the g function is defined by:
.. math:: g(x) = e^{-\theta x} - 1
Args:
X: `np.ndarray`
Returns:
np.array: probability density
"""
self.check_fit()
U, V = split_matrix(X)
if self.theta == 0:
return U * V
else:
num = (-self.theta * self._g(1)) * (1 + self._g(U + V))
aux = self._g(U) * self._g(V) + self._g(1)
den = np.power(aux, 2)
return num / den
[docs] def cumulative_distribution(self, X):
r"""Compute the cumulative distribution function for the Frank copula.
The cumulative density(cdf), or distribution function for the Frank family of copulas
correspond to the formula:
.. math:: C(u,v) = −\frac{\ln({\frac{1 + g(u) g(v)}{g(1)}})}{\theta}
Args:
X: `np.ndarray`
Returns:
np.array: cumulative distribution
"""
self.check_fit()
U, V = split_matrix(X)
num = (np.exp(-self.theta * U) - 1) * (np.exp(-self.theta * V) - 1)
den = np.exp(-self.theta) - 1
return -1.0 / self.theta * np.log(1 + num / den)
[docs] def percent_point(self, y, V):
"""Compute the inverse of conditional cumulative distribution :math:`C(u|v)^{-1}`.
Args:
y: `np.ndarray` value of :math:`C(u|v)`.
v: `np.ndarray` given value of v.
"""
self.check_fit()
if self.theta == 0:
return V
else:
return super().percent_point(y, V)
[docs] def partial_derivative(self, X):
r"""Compute partial derivative of cumulative distribution.
The partial derivative of the copula(CDF) is the conditional CDF.
.. math:: F(v|u) = \frac{\partial}{\partial u}C(u,v) =
\frac{g(u)g(v) + g(v)}{g(u)g(v) + g(1)}
Args:
X (np.ndarray)
y (float)
Returns:
np.ndarray
"""
self.check_fit()
U, V = split_matrix(X)
if self.theta == 0:
return V
else:
num = self._g(U) * self._g(V) + self._g(U)
den = self._g(U) * self._g(V) + self._g(1)
return num / den
[docs] def compute_theta(self):
r"""Compute theta parameter using Kendall's tau.
On Frank copula, the relationship between tau and theta is defined by:
.. math:: \tau = 1 − \frac{4}{\theta} + \frac{4}{\theta^2}\int_0^\theta \!
\frac{t}{e^t -1} \mathrm{d}t.
In order to solve it, we can simplify it as
.. math:: 0 = 1 + \frac{4}{\theta}(D_1(\theta) - 1) - \tau
where the function D is the Debye function of first order, defined as:
.. math:: D_1(x) = \frac{1}{x}\int_0^x\frac{t}{e^t -1} \mathrm{d}t.
"""
result = least_squares(self._tau_to_theta, 1, bounds=(MIN_FLOAT_LOG, MAX_FLOAT_LOG))
return result.x[0]
def _tau_to_theta(self, alpha):
"""Relationship between tau and theta as a solvable equation."""
def debye(t):
return t / (np.exp(t) - 1)
debye_value = integrate.quad(debye, EPSILON, alpha)[0] / alpha
return 4 * (debye_value - 1) / alpha + 1 - self.tau