Skip to content

Commit

Permalink
scipy included
Browse files Browse the repository at this point in the history
  • Loading branch information
@amruthesht
amruthesht committed Jan 21, 2021
1 parent ce2b1cc commit 5c4ad76
Showing 1 changed file with 97 additions and 71 deletions.
168 changes: 97 additions & 71 deletions levy/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@
import numpy as np
from scipy.special import gamma
from scipy import optimize
from scipy.stats import levy_stable

__version__ = "1.1"

Expand Down Expand Up @@ -495,26 +496,80 @@ def levy(x, alpha, beta, mu=0.0, sigma=1.0, cdf=False):
res /= sigma
return float(res) if np.isscalar(x) else res

def modified_levy(x, alpha, beta, mu=0.0, sigma=1.0, a=0, b=np.inf):
if b == np.inf:
return levy(x, alpha, beta, mu=mu, sigma=sigma, cdf=False) / 2 / (1 - levy(a, alpha, beta, mu=mu, sigma=sigma, cdf=True))
else:
return levy(x, alpha, beta, mu=mu, sigma=sigma, cdf=False) / 2 / (levy(b, alpha, beta, mu=mu, sigma=sigma, cdf=True) - levy(a, alpha, beta, mu=mu, sigma=sigma, cdf=True))
def modified_levy(x, alpha, beta, mu=0.0, sigma=1.0, a=0, b=np.inf, cdf=False):
"""
Levy distribution with the tail replaced by the analytical (power law) approximation.
`alpha` in (0, 2] is the index of stability, or characteristic exponent.
`beta` in [-1, 1] is the skewness. `mu` in the reals and `sigma` > 0 are the
location and scale of the distribution (corresponding to `delta` and `gamma`
in Nolan's notation; note that sigma in levy corresponds to sqrt(2) sigma
in the Normal distribution).
'a' and 'b' are the cut off limits for the distribution.
Both a and b are > 0 and the distribution is built on the assumption.
It uses parametrization 0 (to get it from another parametrization, convert).
Example:
>>> x = np.array([1, 2, 3])
>>> modified_levy(x, 1.5, 0, a=0.1, b=3.1)
array([0.23897864, 0.09999616, 0.0372704])
def neglog_levy(x, alpha, beta, mu, sigma):
:param x: values where the function is evaluated
:type x: :class:`~numpy.ndarray`
:param alpha: alpha
:type alpha: float
:param beta: beta
:type beta: float
:param mu: mu
:type mu: float
:param sigma: sigma
:type sigma: float
:param a: a
:type a: float
:param b: b
:type b: float
:return: values of the pdf (or cdf if parameter 'cdf' is set to True) at 'x'
:rtype: :class:`~numpy.ndarray`
"""
Interpolate negative log densities of the Levy stable distribution
specified by `alpha` and `beta`. Small/negative densities are capped
at 1e-100 to preserve sanity.
if x == 0:
raise ValueError
else:
if cdf == False:
if b == np.inf:
return levy(x, alpha, beta, mu=mu, sigma=sigma, cdf=False) / 2 / (1 - levy(a, alpha, beta, mu=mu, sigma=sigma, cdf=True))
else:
return levy(x, alpha, beta, mu=mu, sigma=sigma, cdf=False) / 2 / (levy(b, alpha, beta, mu=mu, sigma=sigma, cdf=True) - levy(a, alpha, beta, mu=mu, sigma=sigma, cdf=True))
else:
if x < 0:
return (levy(x, alpha, beta, mu=mu, sigma=sigma, cdf=True) - levy(-b, alpha, beta, mu=mu, sigma=sigma, cdf=True)) / 2 / (levy(b, alpha, beta, mu=mu, sigma=sigma, cdf=True) - levy(a, alpha, beta, mu=mu, sigma=sigma, cdf=True))
else:
return (levy(x, alpha, beta, mu=mu, sigma=sigma, cdf=True) - levy(a, alpha, beta, mu=mu, sigma=sigma, cdf=True)) / 2 / (levy(b, alpha, beta, mu=mu, sigma=sigma, cdf=True) - levy(a, alpha, beta, mu=mu, sigma=sigma, cdf=True))

def scipy_levy(x, alpha, beta, mu=0.0, sigma=1.0, cdf=False):
if cdf == False:
return levy_stable.pdf(x, alpha, beta, loc=mu, scale=sigma)
else:
return levy_stable.cdf(x, alpha, beta, loc=mu, scale=sigma)

def scipy_modified_levy(x, alpha, beta, mu=0.0, sigma=1.0, a=0, b=np.inf, cdf=False):
"""
Levy distribution with the tail replaced by the analytical (power law) approximation.
`alpha` in (0, 2] is the index of stability, or characteristic exponent.
`beta` in [-1, 1] is the skewness. `mu` in the reals and `sigma` > 0 are the
location and scale of the distribution (corresponding to `delta` and `gamma`
in Nolan's notation; note that sigma in levy corresponds to sqrt(2) sigma
in the Normal distribution).
'a' and 'b' are the cut off limits for the distribution.
Both a and b are > 0 and the distribution is built on the assumption.
It uses parametrization 0 (to get it from another parametrization, convert).
Example:
>>> x = np.array([1,2,3])
>>> neglog_levy(x, 1.5, 0.0, 0.0, 1.0)
array([1.59929892, 2.47054131, 3.45747366])
>>> x = np.array([1, 2, 3])
>>> modified_levy(x, 1.5, 0, a=0.1, b=3.1)
array([0.23897864, 0.09999616, 0.0372704])
:param x: values where the function is evaluated
:type x: :class:`~numpy.ndarray`
Expand All @@ -526,13 +581,29 @@ def neglog_levy(x, alpha, beta, mu, sigma):
:type mu: float
:param sigma: sigma
:type sigma: float
:return: values of -log(pdf(x))
:param a: a
:type a: float
:param b: b
:type b: float
:return: values of the pdf (or cdf if parameter 'cdf' is set to True) at 'x'
:rtype: :class:`~numpy.ndarray`
"""
if x == 0:
raise ValueError
else:
if cdf == False:
if b == np.inf:
return levy_stable.pdf(x, alpha, beta, loc=mu, scale=sigma) / 2 / (1 - levy_stable.cdf(x, alpha, beta, loc=mu, scale=sigma))
else:
return levy_stable.pdf(x, alpha, beta, loc=mu, scale=sigma) / 2 / (levy_stable.cdf(b, alpha, beta, loc=mu, scale=sigma) - levy_stable.cdf(a, alpha, beta, loc=mu, scale=sigma))
else:
if x < 0:
return (levy_stable.cdf(x, alpha, beta, loc=mu, scale=sigma) - levy_stable.cdf(-b, alpha, beta, loc=mu, scale=sigma)) / 2 / (levy_stable.cdf(b, alpha, beta, loc=mu, scale=sigma) - levy_stable.cdf(a, alpha, beta, loc=mu, scale=sigma))
else:
return (levy_stable.cdf(x, alpha, beta, loc=mu, scale=sigma) - levy_stable.cdf(a, alpha, beta, loc=mu, scale=sigma)) / 2 / (levy_stable.cdf(b, alpha, beta, loc=mu, scale=sigma) - levy_stable.cdf(a, alpha, beta, loc=mu, scale=sigma))

return -np.log(np.maximum(1e-100, levy(x, alpha, beta, mu, sigma)))

def neglog_modified_levy(x, alpha, beta, mu, sigma, a, b):
def neglog_levy(x, alpha, beta, mu, sigma, a=0, b=np.inf, dist_type=1):
"""
Interpolate negative log densities of the Levy stable distribution
specified by `alpha` and `beta`. Small/negative densities are capped
Expand All @@ -558,61 +629,16 @@ def neglog_modified_levy(x, alpha, beta, mu, sigma, a, b):
:return: values of -log(pdf(x))
:rtype: :class:`~numpy.ndarray`
"""

return -np.log(np.maximum(1e-100, modified_levy(x, alpha, beta, mu, sigma, a, b)))


def fit_levy(x, par='0', **kwargs):
"""
Estimate parameters of Levy stable distribution given data x, using
Maximum Likelihood estimation.
By default, searches all possible Levy stable distributions. However
you may restrict the search by specifying the values of one or more
parameters. Notice that the parameters to be fixed can be chosen in
all the available parametrizations {0, 1, A, B, M}.
Examples:
>>> np.random.seed(0)
>>> x = random(1.5, 0, 0, 1, shape=(200,))
>>> fit_levy(x) # -- Fit a stable distribution to x
(par=0, alpha=1.52, beta=-0.08, mu=0.05, sigma=0.99, 402.37150603509247)
>>> fit_levy(x, beta=0.0) # -- Fit a symmetric stable distribution to x
(par=0, alpha=1.53, beta=0.00, mu=0.03, sigma=0.99, 402.43833088693725)
>>> fit_levy(x, beta=0.0, mu=0.0) # -- Fit a symmetric distribution centered on zero to x
(par=0, alpha=1.53, beta=0.00, mu=0.00, sigma=0.99, 402.4736618823546)
>>> fit_levy(x, alpha=1.0, beta=0.0) # -- Fit a Cauchy distribution to x
(par=0, alpha=1.00, beta=0.00, mu=0.10, sigma=0.90, 416.54249079255976)
:param x: values to be fitted
:type x: :class:`~numpy.ndarray`
:param par: parametrization
:type par: str
:return: a tuple with a `Parameters` object and the negative log likelihood of the data.
:rtype: tuple
"""

values = {par_name: None if par_name not in kwargs else kwargs[par_name] for i, par_name in
enumerate(par_names[par])}

parameters = Parameters(par=par, **values)
temp = Parameters(par=par, **values)

def neglog_density(param):
temp.x = param
alpha, beta, mu, sigma = temp.get('0')
return np.sum(neglog_levy(x, alpha, beta, mu, sigma))

bounds = tuple(par_bounds[i] for i in parameters.variables)
res = optimize.minimize(neglog_density, parameters.x, method='L-BFGS-B', bounds=bounds)
parameters.x = res.x

return parameters, neglog_density(parameters.x)

def fit_modified_levy(x, a=0, b=np.inf, par='0', **kwargs):
if dist_type == 1:
return -np.log(np.maximum(1e-100, levy(x, alpha, beta, mu, sigma)))
elif dist_type == 2:
return -np.log(np.maximum(1e-100, modified_levy(x, alpha, beta, mu, sigma, a, b)))
elif dist_type == 3:
return -np.log(np.maximum(1e-100, scipy_levy(x, alpha, beta, mu, sigma)))
elif dist_type == 4:
return -np.log(np.maximum(1e-100, scipy_modified_levy(x, alpha, beta, mu, sigma)))

def fit_levy(x, a=0, b=np.inf, par='0', dist_type=1, **kwargs):
"""
Estimate parameters of Levy stable distribution given data x, using
Maximum Likelihood estimation.
Expand Down Expand Up @@ -654,7 +680,7 @@ def fit_modified_levy(x, a=0, b=np.inf, par='0', **kwargs):
def neglog_density(param):
temp.x = param
alpha, beta, mu, sigma = temp.get('0')
return np.sum(neglog_modified_levy(x, alpha, beta, mu, sigma, a, b))
return np.sum(neglog_levy(x, alpha, beta, mu, sigma, a, b, dist_type))

bounds = tuple(par_bounds[i] for i in parameters.variables)
res = optimize.minimize(neglog_density, parameters.x, method='L-BFGS-B', bounds=bounds)
Expand Down

0 comments on commit 5c4ad76

Please sign in to comment.