Source code for signalz.generators.levy_noise
"""
.. versionadded:: 0.5
This function generates random variable with Levy alpha-stable distribution
(also reffered just as stable distribution).
The Levy distribution is defined by two parameters :math:`\\alpha`
and :math:`\\beta`. The Gaussian distribution is special case of
Levy distribution with :math:`\\alpha=2` and :math:`\\beta=0`.
Notes about Implementation
==============================
The implemented algorithm is known as Chambers-Mallows-Stuck method
:cite:`weron1996chambers`. Using two random variables (one with exponential
distribution and one with uniform distribution). The input random variables
are obtained with `numpy.random.exponentinal` and `numpy.random.uniform`.
Example Usage
==================
The following example produce 1000 samples of Levy noise located (mean value)
at -2 (`position`), with characteristic exponent index of 1.5 (`alpha`),
skeewness of 1 (`beta`) and diffusion of 1. (`sigma`).
.. code-block:: python
import signalz
x = signalz.levy_noise(1000, alpha=1.5, beta=0.5, sigma=1., position=-2)
References
============
.. bibliography:: levy_noise.bib
:style: plain
Function Documentation
======================================
"""
import numpy as np
from signalz.misc import check_type_or_raise
[docs]def levy_noise(n, alpha=2., beta=1., sigma=1., position=0.):
"""
This function produces Levy noise.
**Args:**
* `n` - length of the output data (int) - how many samples will be on output
**Kwargs:**
* `alpha` - characteristic exponent index (float) in range `0<alpha<2`
* `beta` - skeewness (float) in range `-1<beta<1`
* `sigma` - diffusion (float), in case of gaussian distribution it is
standard deviation
* `position` - position parameter (float)
**Returns:**
* vector of values representing the noise (1d array)
"""
# correct the inputs or throw error
alpha = float(alpha)
beta = float(beta)
check_type_or_raise(n, int, "n")
if not 0. <= alpha <= 2.:
raise ValueError("Alpha must be between 0 and 2")
if not -1. <= beta <= 1.:
raise ValueError("Beta must be between -1 and 1")
# generate random variables
v = np.random.uniform(-0.5 * np.pi, 0.5 * np.pi, n)
w = np.random.exponential(1, n)
# convert random variables to levy noise
if alpha == 1.:
arg1 = (0.5 * np.pi) + (beta * v)
arg2 = w * np.cos(v)
arg3 = 0.5 * np.pi * beta * sigma + np.log(sigma)
x = 2 / np.pi * (arg1 * np.tan(v) - (beta * np.log(arg2 / arg1)))
return (sigma * x) + arg3 + position
else:
arg1 = 0.5 * np.pi * alpha
b_ab = np.arctan(beta * np.tan(arg1)) / alpha
s_ab = (1 + (beta**2) * np.tan(arg1)**2)**(1/(2.*alpha))
arg2 = alpha * (v + b_ab)
n1 = np.sin(arg2)
d1 = np.cos(v)**(1/alpha)
n2 = np.cos(v - arg2)
x = s_ab * (n1/d1) * (n2/w)**((1-alpha)/alpha)
return (sigma * x) + position