123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913 |
- # -*- coding: utf-8 -*-
- """
- Definition of a hierarchy of classes for kernel functions to be used
- in convolution, e.g., for data smoothing (low pass filtering) or
- firing rate estimation.
- Symmetric kernels
- ~~~~~~~~~~~~~~~~~
- .. autosummary::
- :toctree: toctree/kernels/
- RectangularKernel
- TriangularKernel
- EpanechnikovLikeKernel
- GaussianKernel
- LaplacianKernel
- Asymmetric kernels
- ~~~~~~~~~~~~~~~~~~
- .. autosummary::
- :toctree: toctree/kernels/
- ExponentialKernel
- AlphaKernel
- Examples
- --------
- >>> import quantities as pq
- >>> kernel1 = GaussianKernel(sigma=100*pq.ms)
- >>> kernel2 = ExponentialKernel(sigma=8*pq.ms, invert=True)
- :copyright: Copyright 2016 by the Elephant team, see `doc/authors.rst`.
- :license: Modified BSD, see LICENSE.txt for details.
- """
- from __future__ import division, print_function, unicode_literals
- import math
- import numpy as np
- import quantities as pq
- import scipy.optimize
- import scipy.special
- import scipy.stats
- from elephant.utils import deprecated_alias
- __all__ = [
- 'RectangularKernel', 'TriangularKernel', 'EpanechnikovLikeKernel',
- 'GaussianKernel', 'LaplacianKernel', 'ExponentialKernel', 'AlphaKernel'
- ]
- class Kernel(object):
- r"""
- This is the base class for commonly used kernels.
- **General definition of a kernel:**
- A function :math:`K(x, y)` is called a kernel function if
- :math:`\int{K(x, y) g(x) g(y) \textrm{d}x \textrm{d}y} \ \geq 0 \quad
- \forall g \in L_2`
- **Currently implemented kernels are:**
- * rectangular
- * triangular
- * epanechnikovlike
- * gaussian
- * laplacian
- * exponential (asymmetric)
- * alpha function (asymmetric)
- In neuroscience, a popular application of kernels is in performing
- smoothing operations via convolution. In this case, the kernel has the
- properties of a probability density, i.e., it is positive and normalized
- to one. Popular choices are the rectangular or Gaussian kernels.
- Exponential and alpha kernels may also be used to represent the
- postsynaptic current/potentials in a linear (current-based) model.
- Parameters
- ----------
- sigma : pq.Quantity
- Standard deviation of the kernel.
- invert : bool, optional
- If True, asymmetric kernels (e.g., exponential or alpha kernels) are
- inverted along the time axis.
- Default: False.
- Raises
- ------
- TypeError
- If `sigma` is not `pq.Quantity`.
- If `sigma` is negative.
- If `invert` is not `bool`.
- """
- def __init__(self, sigma, invert=False):
- if not isinstance(sigma, pq.Quantity):
- raise TypeError("'sigma' must be a quantity")
- if sigma.magnitude < 0:
- raise ValueError("'sigma' cannot be negative")
- if not isinstance(invert, bool):
- raise ValueError("'invert' must be bool")
- self.sigma = sigma
- self.invert = invert
- def __repr__(self):
- return "{cls}(sigma={sigma}, invert={invert})".format(
- cls=self.__class__.__name__, sigma=self.sigma, invert=self.invert)
- @deprecated_alias(t='times')
- def __call__(self, times):
- """
- Evaluates the kernel at all points in the array `times`.
- Parameters
- ----------
- times : pq.Quantity
- A vector with time intervals on which the kernel is evaluated.
- Returns
- -------
- pq.Quantity
- Vector with the result of the kernel evaluations.
- Raises
- ------
- TypeError
- If `times` is not `pq.Quantity`.
- If the dimensionality of `times` and :attr:`sigma` are different.
- """
- self._check_time_input(times)
- return self._evaluate(times)
- def _evaluate(self, times):
- """
- Evaluates the kernel Probability Density Function, PDF.
- Parameters
- ----------
- times : pq.Quantity
- Vector with the interval on which the kernel is evaluated, not
- necessarily a time interval.
- Returns
- -------
- pq.Quantity
- Vector with the result of the kernel evaluation.
- """
- raise NotImplementedError(
- "The Kernel class should not be used directly, "
- "instead the subclasses for the single kernels.")
- def boundary_enclosing_area_fraction(self, fraction):
- """
- Calculates the boundary :math:`b` so that the integral from
- :math:`-b` to :math:`b` encloses a certain fraction of the
- integral over the complete kernel.
- By definition the returned value is hence non-negative, even if the
- whole probability mass of the kernel is concentrated over negative
- support for inverted kernels.
- Parameters
- ----------
- fraction : float
- Fraction of the whole area which has to be enclosed.
- Returns
- -------
- pq.Quantity
- Boundary of the kernel containing area `fraction` under the
- kernel density.
- Raises
- ------
- ValueError
- If `fraction` was chosen too close to one, such that in
- combination with integral approximation errors the calculation of
- a boundary was not possible.
- """
- raise NotImplementedError(
- "The Kernel class should not be used directly, "
- "instead the subclasses for the single kernels.")
- def _check_fraction(self, fraction):
- """
- Checks the input variable of the method
- :attr:`boundary_enclosing_area_fraction` for validity of type and
- value.
- Parameters
- ----------
- fraction : float or int
- Fraction of the area under the kernel function.
- Raises
- ------
- TypeError
- If `fraction` is neither a float nor an int.
- If `fraction` is not in the interval [0, 1).
- """
- if not isinstance(fraction, (float, int)):
- raise TypeError("`fraction` must be float or integer")
- if isinstance(self, (TriangularKernel, RectangularKernel)):
- valid = 0 <= fraction <= 1
- bracket = ']'
- else:
- valid = 0 <= fraction < 1
- bracket = ')'
- if not valid:
- raise ValueError("`fraction` must be in the interval "
- "[0, 1{}".format(bracket))
- def _check_time_input(self, t):
- if not isinstance(t, pq.Quantity):
- raise TypeError("The argument 't' of the kernel callable must be "
- "of type Quantity")
- if t.dimensionality.simplified != self.sigma.dimensionality.simplified:
- raise TypeError("The dimensionality of sigma and the input array "
- "to the callable kernel object must be the same. "
- "Otherwise a normalization to 1 of the kernel "
- "cannot be performed.")
- @deprecated_alias(t='time')
- def cdf(self, time):
- r"""
- Cumulative Distribution Function, CDF.
- Parameters
- ----------
- time : pq.Quantity
- The input time scalar.
- Returns
- -------
- float
- CDF at `time`.
- """
- raise NotImplementedError
- def icdf(self, fraction):
- r"""
- Inverse Cumulative Distribution Function, ICDF, also known as a
- quantile.
- Parameters
- ----------
- fraction : float
- The fraction of CDF to compute the quantile from.
- Returns
- -------
- pq.Quantity
- The time scalar `times` such that `CDF(t) = fraction`.
- """
- raise NotImplementedError
- @deprecated_alias(t='times')
- def median_index(self, times):
- r"""
- Estimates the index of the Median of the kernel.
- We define the Median index :math:`i` of a kernel as:
- .. math::
- t_i = \text{ICDF}\left( \frac{\text{CDF}(t_0) +
- \text{CDF}(t_{N-1})}{2} \right)
- where :math:`t_0` and :math:`t_{N-1}` are the first and last entries of
- the input array, CDF and ICDF stand for Cumulative Distribution
- Function and its Inverse, respectively.
- This function is not mandatory for symmetrical kernels but it is
- required when asymmetrical kernels have to be aligned at their median.
- Parameters
- ----------
- times : pq.Quantity
- Vector with the interval on which the kernel is evaluated.
- Returns
- -------
- int
- Index of the estimated value of the kernel median.
- Raises
- ------
- TypeError
- If the input array is not a time pq.Quantity array.
- ValueError
- If the input array is empty.
- If the input array is not sorted.
- See Also
- --------
- Kernel.cdf : cumulative distribution function
- Kernel.icdf : inverse cumulative distribution function
- """
- self._check_time_input(times)
- if len(times) == 0:
- raise ValueError("The input time array is empty.")
- if len(times) <= 2:
- # either left or right; choose left
- return 0
- is_sorted = (np.diff(times.magnitude) >= 0).all()
- if not is_sorted:
- raise ValueError("The input time array must be sorted (in "
- "ascending order).")
- cdf_mean = 0.5 * (self.cdf(times[0]) + self.cdf(times[-1]))
- if cdf_mean == 0.:
- # any index of the kernel non-support is valid; choose median
- return len(times) // 2
- icdf = self.icdf(fraction=cdf_mean)
- icdf = icdf.rescale(times.units).magnitude
- # icdf is guaranteed to be in (t_start, t_end) interval
- median_index = np.nonzero(times.magnitude >= icdf)[0][0]
- return median_index
- def is_symmetric(self):
- r"""
- True for symmetric kernels and False otherwise (asymmetric kernels).
- A kernel is symmetric if its PDF is symmetric w.r.t. time:
- .. math::
- \text{pdf}(-t) = \text{pdf}(t)
- Returns
- -------
- bool
- Whether the kernels is symmetric or not.
- """
- return isinstance(self, SymmetricKernel)
- @property
- def min_cutoff(self):
- """
- Half width of the kernel.
- Returns
- -------
- float
- The returned value varies according to the kernel type.
- """
- raise NotImplementedError
- class SymmetricKernel(Kernel):
- """
- Base class for symmetric kernels.
- """
- class RectangularKernel(SymmetricKernel):
- r"""
- Class for rectangular kernels.
- .. math::
- K(t) = \left\{\begin{array}{ll} \frac{1}{2 \tau}, & |t| < \tau \\
- 0, & |t| \geq \tau \end{array} \right.
- with :math:`\tau = \sqrt{3} \sigma` corresponding to the half width
- of the kernel.
- The parameter `invert` has no effect on symmetric kernels.
- Examples
- --------
- .. plot::
- :include-source:
- from elephant import kernels
- import quantities as pq
- import numpy as np
- import matplotlib.pyplot as plt
- time_array = np.linspace(-3, 3, num=100) * pq.s
- kernel = kernels.RectangularKernel(sigma=1*pq.s)
- kernel_time = kernel(time_array)
- plt.plot(time_array, kernel_time)
- plt.title("RectangularKernel with sigma=1s")
- plt.xlabel("time, s")
- plt.ylabel("kernel, 1/s")
- plt.show()
- """
- @property
- def min_cutoff(self):
- min_cutoff = np.sqrt(3.0)
- return min_cutoff
- def _evaluate(self, times):
- t_units = times.units
- t_abs = np.abs(times.magnitude)
- tau = math.sqrt(3) * self.sigma.rescale(t_units).magnitude
- kernel = (t_abs < tau) * 1 / (2 * tau)
- kernel = pq.Quantity(kernel, units=1 / t_units)
- return kernel
- @deprecated_alias(t='time')
- def cdf(self, time):
- self._check_time_input(time)
- tau = math.sqrt(3) * self.sigma.rescale(time.units).magnitude
- time = np.clip(time.magnitude, a_min=-tau, a_max=tau)
- cdf = (time + tau) / (2 * tau)
- return cdf
- def icdf(self, fraction):
- self._check_fraction(fraction)
- tau = math.sqrt(3) * self.sigma
- icdf = tau * (2 * fraction - 1)
- return icdf
- def boundary_enclosing_area_fraction(self, fraction):
- self._check_fraction(fraction)
- return np.sqrt(3.0) * self.sigma * fraction
- class TriangularKernel(SymmetricKernel):
- r"""
- Class for triangular kernels.
- .. math::
- K(t) = \left\{ \begin{array}{ll} \frac{1}{\tau} (1
- - \frac{|t|}{\tau}), & |t| < \tau \\
- 0, & |t| \geq \tau \end{array} \right.
- with :math:`\tau = \sqrt{6} \sigma` corresponding to the half width of
- the kernel.
- The parameter `invert` has no effect on symmetric kernels.
- Examples
- --------
- .. plot::
- :include-source:
- from elephant import kernels
- import quantities as pq
- import numpy as np
- import matplotlib.pyplot as plt
- time_array = np.linspace(-3, 3, num=1000) * pq.s
- kernel = kernels.TriangularKernel(sigma=1*pq.s)
- kernel_time = kernel(time_array)
- plt.plot(time_array, kernel_time)
- plt.title("TriangularKernel with sigma=1s")
- plt.xlabel("time, s")
- plt.ylabel("kernel, 1/s")
- plt.show()
- """
- @property
- def min_cutoff(self):
- min_cutoff = np.sqrt(6.0)
- return min_cutoff
- def _evaluate(self, times):
- tau = math.sqrt(6) * self.sigma.rescale(times.units).magnitude
- kernel = scipy.stats.triang.pdf(times.magnitude, c=0.5, loc=-tau,
- scale=2 * tau)
- kernel = pq.Quantity(kernel, units=1 / times.units)
- return kernel
- @deprecated_alias(t='time')
- def cdf(self, time):
- self._check_time_input(time)
- tau = math.sqrt(6) * self.sigma.rescale(time.units).magnitude
- cdf = scipy.stats.triang.cdf(time.magnitude, c=0.5, loc=-tau,
- scale=2 * tau)
- return cdf
- def icdf(self, fraction):
- self._check_fraction(fraction)
- tau = math.sqrt(6) * self.sigma.magnitude
- icdf = scipy.stats.triang.ppf(fraction, c=0.5, loc=-tau, scale=2 * tau)
- return icdf * self.sigma.units
- def boundary_enclosing_area_fraction(self, fraction):
- self._check_fraction(fraction)
- return np.sqrt(6.0) * self.sigma * (1 - np.sqrt(1 - fraction))
- class EpanechnikovLikeKernel(SymmetricKernel):
- r"""
- Class for Epanechnikov-like kernels.
- .. math::
- K(t) = \left\{\begin{array}{ll} (3 /(4 d)) (1 - (t / d)^2),
- & |t| < d \\
- 0, & |t| \geq d \end{array} \right.
- with :math:`d = \sqrt{5} \sigma` being the half width of the kernel.
- The Epanechnikov kernel under full consideration of its axioms has a half
- width of :math:`\sqrt{5}`. Ignoring one axiom also the respective kernel
- with half width = 1 can be called Epanechnikov kernel [1]_.
- However, arbitrary width of this type of kernel is here preferred to be
- called 'Epanechnikov-like' kernel.
- The parameter `invert` has no effect on symmetric kernels.
- References
- ----------
- .. [1] https://de.wikipedia.org/wiki/Epanechnikov-Kern
- Examples
- --------
- .. plot::
- :include-source:
- from elephant import kernels
- import quantities as pq
- import numpy as np
- import matplotlib.pyplot as plt
- time_array = np.linspace(-3, 3, num=100) * pq.s
- kernel = kernels.EpanechnikovLikeKernel(sigma=1*pq.s)
- kernel_time = kernel(time_array)
- plt.plot(time_array, kernel_time)
- plt.title("EpanechnikovLikeKernel with sigma=1s")
- plt.xlabel("time, s")
- plt.ylabel("kernel, 1/s")
- plt.show()
- """
- @property
- def min_cutoff(self):
- min_cutoff = np.sqrt(5.0)
- return min_cutoff
- def _evaluate(self, times):
- tau = math.sqrt(5) * self.sigma.rescale(times.units).magnitude
- t_div_tau = np.clip(times.magnitude / tau, a_min=-1, a_max=1)
- kernel = 3. / (4. * tau) * np.maximum(0., 1 - t_div_tau ** 2)
- kernel = pq.Quantity(kernel, units=1 / times.units)
- return kernel
- @deprecated_alias(t='time')
- def cdf(self, time):
- self._check_time_input(time)
- tau = math.sqrt(5) * self.sigma.rescale(time.units).magnitude
- t_div_tau = np.clip(time.magnitude / tau, a_min=-1, a_max=1)
- cdf = 3. / 4 * (t_div_tau - t_div_tau ** 3 / 3.) + 0.5
- return cdf
- def icdf(self, fraction):
- self._check_fraction(fraction)
- # CDF(t) = -1/4 t^3 + 3/4 t + 1/2
- coefs = [-1. / 4, 0, 3. / 4, 0.5 - fraction]
- roots = np.roots(coefs)
- icdf = next(root for root in roots if -1 <= root <= 1)
- tau = math.sqrt(5) * self.sigma
- return icdf * tau
- def boundary_enclosing_area_fraction(self, fraction):
- r"""
- Refer to :func:`Kernel.boundary_enclosing_area_fraction` for the
- documentation.
- Notes
- -----
- For Epanechnikov-like kernels, integration of its density within
- the boundaries 0 and :math:`b`, and then solving for :math:`b` leads
- to the problem of finding the roots of a polynomial of third order.
- The implemented formulas are based on the solution of this problem
- given in [1]_, where the following 3 solutions are given:
- * :math:`u_1 = 1`, solution on negative side;
- * :math:`u_2 = \frac{-1 + i\sqrt{3}}{2}`, solution for larger
- values than zero crossing of the density;
- * :math:`u_3 = \frac{-1 - i\sqrt{3}}{2}`, solution for smaller
- values than zero crossing of the density.
- The solution :math:`u_3` is the relevant one for the problem at hand,
- since it involves only positive area contributions.
- References
- ----------
- .. [1] https://en.wikipedia.org/wiki/Cubic_function
- """
- self._check_fraction(fraction)
- # Python's complex-operator cannot handle quantities, hence the
- # following construction on quantities is necessary:
- Delta_0 = complex(1.0 / (5.0 * self.sigma.magnitude ** 2), 0) / \
- self.sigma.units ** 2
- Delta_1 = complex(2.0 * np.sqrt(5.0) * fraction /
- (25.0 * self.sigma.magnitude ** 3), 0) / \
- self.sigma.units ** 3
- C = ((Delta_1 + (Delta_1 ** 2.0 - 4.0 * Delta_0 ** 3.0) ** (
- 1.0 / 2.0)) /
- 2.0) ** (1.0 / 3.0)
- u_3 = complex(-1.0 / 2.0, -np.sqrt(3.0) / 2.0)
- b = -5.0 * self.sigma ** 2 * (u_3 * C + Delta_0 / (u_3 * C))
- return b.real
- class GaussianKernel(SymmetricKernel):
- r"""
- Class for gaussian kernels.
- .. math::
- K(t) = (\frac{1}{\sigma \sqrt{2 \pi}}) \exp(-\frac{t^2}{2 \sigma^2})
- with :math:`\sigma` being the standard deviation.
- The parameter `invert` has no effect on symmetric kernels.
- Examples
- --------
- .. plot::
- :include-source:
- from elephant import kernels
- import quantities as pq
- import numpy as np
- import matplotlib.pyplot as plt
- time_array = np.linspace(-3, 3, num=100) * pq.s
- kernel = kernels.GaussianKernel(sigma=1*pq.s)
- kernel_time = kernel(time_array)
- plt.plot(time_array, kernel_time)
- plt.title("GaussianKernel with sigma=1s")
- plt.xlabel("time, s")
- plt.ylabel("kernel, 1/s")
- plt.show()
- """
- @property
- def min_cutoff(self):
- min_cutoff = 3.0
- return min_cutoff
- def _evaluate(self, times):
- sigma = self.sigma.rescale(times.units).magnitude
- kernel = scipy.stats.norm.pdf(times.magnitude, loc=0, scale=sigma)
- kernel = pq.Quantity(kernel, units=1 / times.units)
- return kernel
- @deprecated_alias(t='time')
- def cdf(self, time):
- self._check_time_input(time)
- sigma = self.sigma.rescale(time.units).magnitude
- cdf = scipy.stats.norm.cdf(time, loc=0, scale=sigma)
- return cdf
- def icdf(self, fraction):
- self._check_fraction(fraction)
- icdf = scipy.stats.norm.ppf(fraction, loc=0,
- scale=self.sigma.magnitude)
- return icdf * self.sigma.units
- def boundary_enclosing_area_fraction(self, fraction):
- self._check_fraction(fraction)
- return self.sigma * np.sqrt(2.0) * scipy.special.erfinv(fraction)
- class LaplacianKernel(SymmetricKernel):
- r"""
- Class for laplacian kernels.
- .. math::
- K(t) = \frac{1}{2 \tau} \exp\left(-\left|\frac{t}{\tau}\right|\right)
- with :math:`\tau = \sigma / \sqrt{2}`.
- The parameter `invert` has no effect on symmetric kernels.
- Examples
- --------
- .. plot::
- :include-source:
- from elephant import kernels
- import quantities as pq
- import numpy as np
- import matplotlib.pyplot as plt
- time_array = np.linspace(-3, 3, num=1000) * pq.s
- kernel = kernels.LaplacianKernel(sigma=1*pq.s)
- kernel_time = kernel(time_array)
- plt.plot(time_array, kernel_time)
- plt.title("LaplacianKernel with sigma=1s")
- plt.xlabel("time, s")
- plt.ylabel("kernel, 1/s")
- plt.show()
- """
- @property
- def min_cutoff(self):
- min_cutoff = 3.0
- return min_cutoff
- def _evaluate(self, times):
- tau = self.sigma.rescale(times.units).magnitude / math.sqrt(2)
- kernel = scipy.stats.laplace.pdf(times.magnitude, loc=0, scale=tau)
- kernel = pq.Quantity(kernel, units=1 / times.units)
- return kernel
- @deprecated_alias(t='time')
- def cdf(self, time):
- self._check_time_input(time)
- tau = self.sigma.rescale(time.units).magnitude / math.sqrt(2)
- cdf = scipy.stats.laplace.cdf(time.magnitude, loc=0, scale=tau)
- return cdf
- def icdf(self, fraction):
- self._check_fraction(fraction)
- tau = self.sigma.magnitude / math.sqrt(2)
- icdf = scipy.stats.laplace.ppf(fraction, loc=0, scale=tau)
- return icdf * self.sigma.units
- def boundary_enclosing_area_fraction(self, fraction):
- self._check_fraction(fraction)
- return -self.sigma * np.log(1.0 - fraction) / np.sqrt(2.0)
- # Potential further symmetric kernels from Wiki Kernels (statistics):
- # Quartic (biweight), Triweight, Tricube, Cosine, Logistics, Silverman
- class ExponentialKernel(Kernel):
- r"""
- Class for exponential kernels.
- .. math::
- K(t) = \left\{\begin{array}{ll} (1 / \tau) \exp{(-t / \tau)},
- & t \geq 0 \\
- 0, & t < 0 \end{array} \right.
- with :math:`\tau = \sigma`.
- Examples
- --------
- .. plot::
- :include-source:
- from elephant import kernels
- import quantities as pq
- import numpy as np
- import matplotlib.pyplot as plt
- time_array = np.linspace(-1, 4, num=100) * pq.s
- kernel = kernels.ExponentialKernel(sigma=1*pq.s)
- kernel_time = kernel(time_array)
- plt.plot(time_array, kernel_time)
- plt.title("ExponentialKernel with sigma=1s")
- plt.xlabel("time, s")
- plt.ylabel("kernel, 1/s")
- plt.show()
- """
- @property
- def min_cutoff(self):
- min_cutoff = 3.0
- return min_cutoff
- def _evaluate(self, times):
- tau = self.sigma.rescale(times.units).magnitude
- if self.invert:
- times = -times
- kernel = scipy.stats.expon.pdf(times.magnitude, loc=0, scale=tau)
- kernel = pq.Quantity(kernel, units=1 / times.units)
- return kernel
- @deprecated_alias(t='time')
- def cdf(self, time):
- self._check_time_input(time)
- tau = self.sigma.rescale(time.units).magnitude
- time = time.magnitude
- if self.invert:
- time = np.minimum(time, 0)
- return np.exp(time / tau)
- time = np.maximum(time, 0)
- return 1. - np.exp(-time / tau)
- def icdf(self, fraction):
- self._check_fraction(fraction)
- if self.invert:
- return self.sigma * np.log(fraction)
- return -self.sigma * np.log(1.0 - fraction)
- def boundary_enclosing_area_fraction(self, fraction):
- # the boundary b, which encloses a 'fraction' of CDF in [-b, b] range,
- # does not depend on the invert, if the kernel is cut at zero.
- # It's easier to compute 'b' for a kernel that has not been inverted.
- kernel = self.__class__(sigma=self.sigma, invert=False)
- return kernel.icdf(fraction)
- class AlphaKernel(Kernel):
- r"""
- Class for alpha kernels.
- .. math::
- K(t) = \left\{\begin{array}{ll} (1 / \tau^2)
- \ t\ \exp{(-t / \tau)}, & t > 0 \\
- 0, & t \leq 0 \end{array} \right.
- with :math:`\tau = \sigma / \sqrt{2}`.
- Examples
- --------
- .. plot::
- :include-source:
- from elephant import kernels
- import quantities as pq
- import numpy as np
- import matplotlib.pyplot as plt
- time_array = np.linspace(-1, 4, num=100) * pq.s
- kernel = kernels.AlphaKernel(sigma=1*pq.s)
- kernel_time = kernel(time_array)
- plt.plot(time_array, kernel_time)
- plt.title("AlphaKernel with sigma=1s")
- plt.xlabel("time, s")
- plt.ylabel("kernel, 1/s")
- plt.show()
- """
- @property
- def min_cutoff(self):
- min_cutoff = 3.0
- return min_cutoff
- def _evaluate(self, times):
- t_units = times.units
- tau = self.sigma.rescale(t_units).magnitude / math.sqrt(2)
- times = times.magnitude
- if self.invert:
- times = -times
- kernel = (times >= 0) * 1 / tau ** 2 * times * np.exp(-times / tau)
- kernel = pq.Quantity(kernel, units=1 / t_units)
- return kernel
- @deprecated_alias(t='time')
- def cdf(self, time):
- self._check_time_input(time)
- tau = self.sigma.rescale(time.units).magnitude / math.sqrt(2)
- cdf = self._cdf_stripped(time.magnitude, tau)
- return cdf
- def _cdf_stripped(self, t, tau):
- # CDF without time units
- if self.invert:
- t = np.minimum(t, 0)
- return np.exp(t / tau) * (tau - t) / tau
- t = np.maximum(t, 0)
- return 1 - np.exp(-t / tau) * (t + tau) / tau
- def icdf(self, fraction):
- self._check_fraction(fraction)
- tau = self.sigma.magnitude / math.sqrt(2)
- def cdf(x):
- # CDF fof the AlphaKernel, subtracted 'fraction'
- # evaluates the error of the root of cdf(x) = fraction
- return self._cdf_stripped(x, tau) - fraction
- # fraction is a good starting point for CDF approximation
- x0 = fraction if not self.invert else fraction - 1
- x_quantile = scipy.optimize.fsolve(cdf, x0=x0, xtol=1e-7)[0]
- x_quantile = pq.Quantity(x_quantile, units=self.sigma.units)
- return x_quantile
- def boundary_enclosing_area_fraction(self, fraction):
- # the boundary b, which encloses a 'fraction' of CDF in [-b, b] range,
- # does not depend on the invert, if the kernel is cut at zero.
- # It's easier to compute 'b' for a kernel that has not been inverted.
- kernel = self.__class__(sigma=self.sigma, invert=False)
- return kernel.icdf(fraction)
|