diff --git a/levy/__init__.py b/levy/__init__.py index 5217159..9fed9ca 100644 --- a/levy/__init__.py +++ b/levy/__init__.py @@ -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" @@ -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` @@ -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 @@ -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. @@ -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)