qp Demo
Note: This is a legacy notebook created for qp version < 1. Some of this notebook may no longer be accurate for qp version >= 1. It is kept availablle for reference.
Alex Malz, Phil Marshall, Eric Charles
In this notebook we use the qp module to approximate some simple, standard, 1D PDFs using sets of quantiles, samples, and histograms, and assess their relative accuracy.
We also show how such analyses can be extended to use “composite” PDFs made up of mixtures of standard distributions.
import numpy as np
import os
import scipy.stats as sps
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
%matplotlib inline
Requirements
To run qp, you will need to first install the module by following the instructions here.
import qp
Background: the scipy.stats module
The scipy.stats module is the standard for manipulating distribtions so is a natural place to start for implementing 1D PDF parameterizations.
It allows you do define a wide variety of distibutions and uses numpy array broadcasting for efficiency.
Gaussian (Normal) example
Here are some examples of things you can do with the scipy.stats module, using a Gaussian or Normal distribution.
loc and scale are the means and standard deviations of the underlying Gaussians.
Note the distinction between passing arguments to norm and passing arguments to pdf to access multiple distributions and their PDF values at multiple points.
# evaluate a single distribution's PDF at one value
print("PDF at one point for one distribution:",
sps.norm(loc=0, scale=1).pdf(0.5))
# evaluate a single distribution's PDF at multiple value
print("PDF at three points for one distribution:",
sps.norm(loc=0, scale=1).pdf([0.5, 1., 1.5]))
# evalute three distributions' PDFs at one shared value
print("PDF at one point for three distributions:",
sps.norm(loc=[0., 1., 2.], scale=1).pdf(0.5))
# evalute three distributions' PDFs each at one different value
print("PDF at one different point for three distributions:",
sps.norm(loc=[0., 1., 2.], scale=1).pdf([0.5, 1., 1.5]))
# evalute three distributions' PDFs each at four different values
# (note the change in shape of the argument)
print("PDF at four different points for three distributions:\n",
sps.norm(loc=[0., 1., 2.], scale=1).pdf([[0.5],[1.],[1.5],[2]]))
# evalute three distributions' PDFs at each of four different values
# (note the change in shape of the argument)
print("PDF at four different points for three distributions: broadcast reversed\n",
sps.norm(loc=[[0.], [1.], [2.]], scale=1).pdf([0.5,1.,1.5,2]))
PDF at one point for one distribution: 0.35206532676429947
PDF at three points for one distribution: [0.35206533 0.24197072 0.1295176 ]
PDF at one point for three distributions: [0.35206533 0.35206533 0.1295176 ]
PDF at one different point for three distributions: [0.35206533 0.39894228 0.35206533]
PDF at four different points for three distributions:
[[0.35206533 0.35206533 0.1295176 ]
[0.24197072 0.39894228 0.24197072]
[0.1295176 0.35206533 0.35206533]
[0.05399097 0.24197072 0.39894228]]
PDF at four different points for three distributions: broadcast reversed
[[0.35206533 0.24197072 0.1295176 0.05399097]
[0.35206533 0.39894228 0.35206533 0.24197072]
[0.1295176 0.24197072 0.35206533 0.39894228]]
The scipy.stats classes
In the scipy.stats module, all of the distributions are sub-classes of scipy.stats.rv_continuous.
You make an object of a particular sub-type, and then ‘freeze’ it by passing it shape parameters.
print("This is the generic normal distribution class: ",
sps._continuous_distns.norm_gen)
ng = sps._continuous_distns.norm_gen()
print("This is an instance of the generic normal distribution class",
ng)
norm_sp = ng(loc=0, scale=1)
print("This is a frozen normal distribution, with specific paramters",
norm_sp, norm_sp.kwds)
print("The frozen object know what generic distribution it comes from",
norm_sp.dist)
This is the generic normal distribution class: <class 'scipy.stats._continuous_distns.norm_gen'>
This is an instance of the generic normal distribution class <scipy.stats._continuous_distns.norm_gen object at 0x7f54b6557d00>
This is a frozen normal distribution, with specific paramters <scipy.stats._distn_infrastructure.rv_continuous_frozen object at 0x7f54b6556320> {'loc': 0, 'scale': 1}
The frozen object know what generic distribution it comes from <scipy.stats._continuous_distns.norm_gen object at 0x7f54b6557730>
Properties of distributions
scipy.stats lets you evaluate multiple properties of distributions. These include:
pdf: Probability Density Function
cdf: Cumulative Distribution Function
ppf: Percent Point Function (Inverse of CDF)
sf: Survival Function (1-CDF)
isf: Inverse Survival Function (Inverse of SF)
rvs: Random Variates (i.e., sampled values)
stats: Return mean, variance, optionally: (Fisher’s) skew, or (Fisher’s) kurtosis
moment: non-central moments of the distribution
print("PDF = ", norm_sp.pdf(0.5))
print("CDF = ", norm_sp.cdf(0.5))
print("PPF = ", norm_sp.ppf(0.6))
print("SF = ", norm_sp.sf(0.6))
print("ISF = ", norm_sp.isf(0.5))
print("RVS = ", norm_sp.rvs())
print("stats = ", norm_sp.stats())
print("M2 = ", norm_sp.moment(2))
PDF = 0.35206532676429947
CDF = 0.6914624612740131
PPF = 0.2533471031357997
SF = 0.2742531177500736
ISF = 0.0
RVS = -0.9926012323717978
stats = (np.float64(0.0), np.float64(1.0))
M2 = 1.0
qp parameterizations and visualization functionality
The next part of this notebook shows how we can extend the functionality of scipy.stats to implement distributions that are based on parameterizations of 1D PDFs, like histograms, interpolations, splines, or mixture models.
Parameterizations from scipy.stats
qp automatically generates classes for all of the scipy.stats.rv_continuous distributions, providing feed-through access to all scipy.stats.rv_continuous objects but adds on additional attributes and methods specific to parameterization conversions.
qp.stats.keys()
odict_keys(['alpha', 'anglit', 'arcsine', 'argus', 'beta', 'betaprime', 'bradford', 'burr', 'burr12', 'cauchy', 'chi', 'chi2', 'cosine', 'crystalball', 'dgamma', 'dpareto_lognorm', 'dweibull', 'erlang', 'expon', 'exponnorm', 'exponpow', 'exponweib', 'f', 'fatiguelife', 'fisk', 'foldcauchy', 'foldnorm', 'gamma', 'gausshyper', 'genexpon', 'genextreme', 'gengamma', 'genhalflogistic', 'genhyperbolic', 'geninvgauss', 'genlogistic', 'gennorm', 'genpareto', 'gibrat', 'gompertz', 'gumbel_l', 'gumbel_r', 'halfcauchy', 'halfgennorm', 'halflogistic', 'halfnorm', 'hypsecant', 'invgamma', 'invgauss', 'invweibull', 'irwinhall', 'jf_skew_t', 'johnsonsb', 'johnsonsu', 'kappa3', 'kappa4', 'ksone', 'kstwo', 'kstwobign', 'landau', 'laplace', 'laplace_asymmetric', 'levy', 'levy_l', 'levy_stable', 'loggamma', 'logistic', 'loglaplace', 'lognorm', 'loguniform', 'lomax', 'maxwell', 'mielke', 'moyal', 'nakagami', 'ncf', 'nct', 'ncx2', 'norm', 'norminvgauss', 'pareto', 'pearson3', 'powerlaw', 'powerlognorm', 'powernorm', 'rayleigh', 'rdist', 'recipinvgauss', 'reciprocal', 'rel_breitwigner', 'rice', 'semicircular', 'skewcauchy', 'skewnorm', 'studentized_range', 't', 'trapezoid', 'trapz', 'triang', 'truncexpon', 'truncnorm', 'truncpareto', 'truncweibull_min', 'tukeylambda', 'uniform', 'vonmises', 'vonmises_line', 'wald', 'weibull_max', 'weibull_min', 'wrapcauchy', 'spline', 'hist', 'interp', 'interp_irregular', 'quant', 'mixmod', 'sparse', 'packed_interp'])
help(qp.stats.lognorm_gen)
Help on class lognorm in module qp.core.factory:
class lognorm(qp.parameterizations.base.Pdf_gen_wrap, scipy.stats._continuous_distns.lognorm_gen)
| lognorm(*args, **kwargs)
|
| A lognormal continuous random variable.
|
| %(before_notes)s
|
| Notes
| -----
| The probability density function for `lognorm` is:
|
| .. math::
|
| f(x, s) = \frac{1}{s x \sqrt{2\pi}}
| \exp\left(-\frac{\log^2(x)}{2s^2}\right)
|
| for :math:`x > 0`, :math:`s > 0`.
|
| `lognorm` takes ``s`` as a shape parameter for :math:`s`.
|
| %(after_notes)s
|
| Suppose a normally distributed random variable ``X`` has mean ``mu`` and
| standard deviation ``sigma``. Then ``Y = exp(X)`` is lognormally
| distributed with ``s = sigma`` and ``scale = exp(mu)``.
|
| %(example)s
|
| The logarithm of a log-normally distributed random variable is
| normally distributed:
|
| >>> import numpy as np
| >>> import matplotlib.pyplot as plt
| >>> from scipy import stats
| >>> fig, ax = plt.subplots(1, 1)
| >>> mu, sigma = 2, 0.5
| >>> X = stats.norm(loc=mu, scale=sigma)
| >>> Y = stats.lognorm(s=sigma, scale=np.exp(mu))
| >>> x = np.linspace(*X.interval(0.999))
| >>> y = Y.rvs(size=10000)
| >>> ax.plot(x, X.pdf(x), label='X (pdf)')
| >>> ax.hist(np.log(y), density=True, bins=x, label='log(Y) (histogram)')
| >>> ax.legend()
| >>> plt.show()
|
| Method resolution order:
| lognorm
| qp.parameterizations.base.Pdf_gen_wrap
| qp.parameterizations.base.Pdf_gen
| scipy.stats._continuous_distns.lognorm_gen
| scipy.stats._distn_infrastructure.rv_continuous
| scipy.stats._distn_infrastructure.rv_generic
| builtins.object
|
| Methods defined here:
|
| create_ensemble(data: 'Mapping', ancil: 'Optional[Mapping]' = None) -> 'Ensemble'
| Creates an Ensemble of distribution(s) in the given parameterization.
|
| Input data format:
| data = {'arg1': values, 'arg2': values ...} where 'arg1', 'arg2'... are the arguments for the parameterization.
| The length of the values should be the number of distributions being created in the Ensemble, with a minimum value of 1.
|
|
| Parameters
| ----------
| data : Mapping
| The dictionary of data for the distributions.
| ancil : Optional[Mapping], optional
| A dictionary of metadata for the distributions, where any arrays have the same length as the number of distributions, by default None
|
| Returns
| -------
| Ensemble
| An Ensemble object containing all of the given distributions.
|
| Examples
| --------
|
| To create an Ensemble with two Gaussian distributions and their associated ids:
|
| >>> import qp
| >>> data = {'loc': np.array([[0.45],[0.55]]) , 'scale': np.array([[0.2],[0.15]])}
| >>> ancil = {'ids': [20,25]}
| >>> ens = qp.stats.norm.create_ensemble(data,ancil)
| >>> ens.metadata
| {'pdf_name': array([b'norm'], dtype='|S4'), 'pdf_version': array([0])}
|
| freeze = _my_freeze(self, *args, **kwds)
|
| ----------------------------------------------------------------------
| Data and other attributes defined here:
|
| name = 'lognorm'
|
| version = 0
|
| ----------------------------------------------------------------------
| Methods inherited from qp.parameterizations.base.Pdf_gen_wrap:
|
| __init__(self, *args, **kwargs)
| C'tor
|
| ----------------------------------------------------------------------
| Class methods inherited from qp.parameterizations.base.Pdf_gen_wrap:
|
| add_mappings() from builtins.type
| Add this classes mappings to the conversion dictionary
|
| create(**kwds) from builtins.type
| Create and return a `scipy.stats.rv_frozen` object using the
| keyword arguments provided
|
| create_gen(**kwds) from builtins.type
| Create and return a `scipy.stats.rv_continuous` object using the
| keyword arguments provided
|
| get_allocation_kwds(npdf, **kwargs) from builtins.type
| Return kwds necessary to create 'empty' hdf5 file with npdf entries
| for iterative writeout
|
| ----------------------------------------------------------------------
| Class methods inherited from qp.parameterizations.base.Pdf_gen:
|
| add_method_dicts() from builtins.type
| Add empty method dicts
|
| creation_method(method=None) from builtins.type
| Return the method used to create a PDF of this type
|
| extraction_method(method=None) from builtins.type
| Return the method used to extract data to create a PDF of this type
|
| plot(pdf, **kwargs) from builtins.type
| Plot the pdf as a curve
|
| plot_native(pdf, **kwargs) from builtins.type
| Plot the PDF in a way that is particular to this type of distribution
|
| This defaults to plotting it as a curve, but this can be overwritten
|
| print_method_maps(stream=<ipykernel.iostream.OutStream object at 0x7f54f45ac340>) from builtins.type
| Print the maps showing the methods
|
| reader_method(version=None) from builtins.type
| Return the method used to convert data read from a file PDF of this type
|
| ----------------------------------------------------------------------
| Readonly properties inherited from qp.parameterizations.base.Pdf_gen:
|
| metadata
| Return the metadata for this set of PDFs
|
| objdata
| Return the object data for this set of PDFs
|
| ----------------------------------------------------------------------
| Data descriptors inherited from qp.parameterizations.base.Pdf_gen:
|
| __dict__
| dictionary for instance variables (if defined)
|
| __weakref__
| list of weak references to the object (if defined)
|
| ----------------------------------------------------------------------
| Methods inherited from scipy.stats._continuous_distns.lognorm_gen:
|
| fit(self, data, *args, **kwds)
| Return estimates of shape (if applicable), location, and scale
| parameters from data. The default estimation method is Maximum
| Likelihood Estimation (MLE), but Method of Moments (MM)
| is also available.
|
| Starting estimates for the fit are given by input arguments;
| for any arguments not provided with starting estimates,
| ``self._fitstart(data)`` is called to generate such.
|
| One can hold some parameters fixed to specific values by passing in
| keyword arguments ``f0``, ``f1``, ..., ``fn`` (for shape parameters)
| and ``floc`` and ``fscale`` (for location and scale parameters,
| respectively).
|
| Parameters
| ----------
| data : array_like or `CensoredData` instance
| Data to use in estimating the distribution parameters.
| arg1, arg2, arg3,... : floats, optional
| Starting value(s) for any shape-characterizing arguments (those not
| provided will be determined by a call to ``_fitstart(data)``).
| No default value.
| **kwds : floats, optional
| - `loc`: initial guess of the distribution's location parameter.
| - `scale`: initial guess of the distribution's scale parameter.
|
| Special keyword arguments are recognized as holding certain
| parameters fixed:
|
| - f0...fn : hold respective shape parameters fixed.
| Alternatively, shape parameters to fix can be specified by name.
| For example, if ``self.shapes == "a, b"``, ``fa`` and ``fix_a``
| are equivalent to ``f0``, and ``fb`` and ``fix_b`` are
| equivalent to ``f1``.
|
| - floc : hold location parameter fixed to specified value.
|
| - fscale : hold scale parameter fixed to specified value.
|
| - optimizer : The optimizer to use. The optimizer must take
| ``func`` and starting position as the first two arguments,
| plus ``args`` (for extra arguments to pass to the
| function to be optimized) and ``disp``.
| The ``fit`` method calls the optimizer with ``disp=0`` to suppress output.
| The optimizer must return the estimated parameters.
|
| - method : The method to use. The default is "MLE" (Maximum
| Likelihood Estimate); "MM" (Method of Moments)
| is also available.
|
| Raises
| ------
| TypeError, ValueError
| If an input is invalid
| `~scipy.stats.FitError`
| If fitting fails or the fit produced would be invalid
|
| Returns
| -------
| parameter_tuple : tuple of floats
| Estimates for any shape parameters (if applicable), followed by
| those for location and scale. For most random variables, shape
| statistics will be returned, but there are exceptions (e.g.
| ``norm``).
|
| Notes
| -----
| With ``method="MLE"`` (default), the fit is computed by minimizing
| the negative log-likelihood function. A large, finite penalty
| (rather than infinite negative log-likelihood) is applied for
| observations beyond the support of the distribution.
|
| With ``method="MM"``, the fit is computed by minimizing the L2 norm
| of the relative errors between the first *k* raw (about zero) data
| moments and the corresponding distribution moments, where *k* is the
| number of non-fixed parameters.
| More precisely, the objective function is::
|
| (((data_moments - dist_moments)
| / np.maximum(np.abs(data_moments), 1e-8))**2).sum()
|
| where the constant ``1e-8`` avoids division by zero in case of
| vanishing data moments. Typically, this error norm can be reduced to
| zero.
| Note that the standard method of moments can produce parameters for
| which some data are outside the support of the fitted distribution;
| this implementation does nothing to prevent this.
|
| For either method,
| the returned answer is not guaranteed to be globally optimal; it
| may only be locally optimal, or the optimization may fail altogether.
| If the data contain any of ``np.nan``, ``np.inf``, or ``-np.inf``,
| the `fit` method will raise a ``RuntimeError``.
|
| When passing a ``CensoredData`` instance to ``data``, the log-likelihood
| function is defined as:
|
| .. math::
|
| l(\pmb{\theta}; k) & = \sum
| \log(f(k_u; \pmb{\theta}))
| + \sum
| \log(F(k_l; \pmb{\theta})) \\
| & + \sum
| \log(1 - F(k_r; \pmb{\theta})) \\
| & + \sum
| \log(F(k_{\text{high}, i}; \pmb{\theta})
| - F(k_{\text{low}, i}; \pmb{\theta}))
|
| where :math:`f` and :math:`F` are the pdf and cdf, respectively, of the
| function being fitted, :math:`\pmb{\theta}` is the parameter vector,
| :math:`u` are the indices of uncensored observations,
| :math:`l` are the indices of left-censored observations,
| :math:`r` are the indices of right-censored observations,
| subscripts "low"/"high" denote endpoints of interval-censored observations, and
| :math:`i` are the indices of interval-censored observations.
|
| When `method='MLE'` and
| the location parameter is fixed by using the `floc` argument,
| this function uses explicit formulas for the maximum likelihood
| estimation of the log-normal shape and scale parameters, so the
| `optimizer`, `loc` and `scale` keyword arguments are ignored.
| If the location is free, a likelihood maximum is found by
| setting its partial derivative wrt to location to 0, and
| solving by substituting the analytical expressions of shape
| and scale (or provided parameters).
| See, e.g., equation 3.1 in
| A. Clifford Cohen & Betty Jones Whitten (1980)
| Estimation in the Three-Parameter Lognormal Distribution,
| Journal of the American Statistical Association, 75:370, 399-404
| https://doi.org/10.2307/2287466
|
|
| Examples
| --------
|
| Generate some data to fit: draw random variates from the `beta`
| distribution
|
| >>> import numpy as np
| >>> from scipy.stats import beta
| >>> a, b = 1., 2.
| >>> rng = np.random.default_rng(172786373191770012695001057628748821561)
| >>> x = beta.rvs(a, b, size=1000, random_state=rng)
|
| Now we can fit all four parameters (``a``, ``b``, ``loc`` and
| ``scale``):
|
| >>> a1, b1, loc1, scale1 = beta.fit(x)
| >>> a1, b1, loc1, scale1
| (1.0198945204435628, 1.9484708982737828, 4.372241314917588e-05, 0.9979078845964814)
|
| The fit can be done also using a custom optimizer:
|
| >>> from scipy.optimize import minimize
| >>> def custom_optimizer(func, x0, args=(), disp=0):
| ... res = minimize(func, x0, args, method="slsqp", options={"disp": disp})
| ... if res.success:
| ... return res.x
| ... raise RuntimeError('optimization routine failed')
| >>> a1, b1, loc1, scale1 = beta.fit(x, method="MLE", optimizer=custom_optimizer)
| >>> a1, b1, loc1, scale1
| (1.0198821087258905, 1.948484145914738, 4.3705304486881485e-05, 0.9979104663953395)
|
| We can also use some prior knowledge about the dataset: let's keep
| ``loc`` and ``scale`` fixed:
|
| >>> a1, b1, loc1, scale1 = beta.fit(x, floc=0, fscale=1)
| >>> loc1, scale1
| (0, 1)
|
| We can also keep shape parameters fixed by using ``f``-keywords. To
| keep the zero-th shape parameter ``a`` equal 1, use ``f0=1`` or,
| equivalently, ``fa=1``:
|
| >>> a1, b1, loc1, scale1 = beta.fit(x, fa=1, floc=0, fscale=1)
| >>> a1
| 1
|
| Not all distributions return estimates for the shape parameters.
| ``norm`` for example just returns estimates for location and scale:
|
| >>> from scipy.stats import norm
| >>> x = norm.rvs(a, b, size=1000, random_state=123)
| >>> loc1, scale1 = norm.fit(x)
| >>> loc1, scale1
| (0.92087172783841631, 2.0015750750324668)
|
| ----------------------------------------------------------------------
| Methods inherited from scipy.stats._distn_infrastructure.rv_continuous:
|
| __getstate__(self)
|
| cdf(self, x, *args, **kwds)
| Cumulative distribution function of the given RV.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| cdf : ndarray
| Cumulative distribution function evaluated at `x`
|
| expect(self, func=None, args=(), loc=0, scale=1, lb=None, ub=None, conditional=False, **kwds)
| Calculate expected value of a function with respect to the
| distribution by numerical integration.
|
| The expected value of a function ``f(x)`` with respect to a
| distribution ``dist`` is defined as::
|
| ub
| E[f(x)] = Integral(f(x) * dist.pdf(x)),
| lb
|
| where ``ub`` and ``lb`` are arguments and ``x`` has the ``dist.pdf(x)``
| distribution. If the bounds ``lb`` and ``ub`` correspond to the
| support of the distribution, e.g. ``[-inf, inf]`` in the default
| case, then the integral is the unrestricted expectation of ``f(x)``.
| Also, the function ``f(x)`` may be defined such that ``f(x)`` is ``0``
| outside a finite interval in which case the expectation is
| calculated within the finite range ``[lb, ub]``.
|
| Parameters
| ----------
| func : callable, optional
| Function for which integral is calculated. Takes only one argument.
| The default is the identity mapping f(x) = x.
| args : tuple, optional
| Shape parameters of the distribution.
| loc : float, optional
| Location parameter (default=0).
| scale : float, optional
| Scale parameter (default=1).
| lb, ub : scalar, optional
| Lower and upper bound for integration. Default is set to the
| support of the distribution.
| conditional : bool, optional
| If True, the integral is corrected by the conditional probability
| of the integration interval. The return value is the expectation
| of the function, conditional on being in the given interval.
| Default is False.
|
| Additional keyword arguments are passed to the integration routine.
|
| Returns
| -------
| expect : float
| The calculated expected value.
|
| Notes
| -----
| The integration behavior of this function is inherited from
| `scipy.integrate.quad`. Neither this function nor
| `scipy.integrate.quad` can verify whether the integral exists or is
| finite. For example ``cauchy(0).mean()`` returns ``np.nan`` and
| ``cauchy(0).expect()`` returns ``0.0``.
|
| Likewise, the accuracy of results is not verified by the function.
| `scipy.integrate.quad` is typically reliable for integrals that are
| numerically favorable, but it is not guaranteed to converge
| to a correct value for all possible intervals and integrands. This
| function is provided for convenience; for critical applications,
| check results against other integration methods.
|
| The function is not vectorized.
|
| Examples
| --------
|
| To understand the effect of the bounds of integration consider
|
| >>> from scipy.stats import expon
| >>> expon(1).expect(lambda x: 1, lb=0.0, ub=2.0)
| 0.6321205588285578
|
| This is close to
|
| >>> expon(1).cdf(2.0) - expon(1).cdf(0.0)
| 0.6321205588285577
|
| If ``conditional=True``
|
| >>> expon(1).expect(lambda x: 1, lb=0.0, ub=2.0, conditional=True)
| 1.0000000000000002
|
| The slight deviation from 1 is due to numerical integration.
|
| The integrand can be treated as a complex-valued function
| by passing ``complex_func=True`` to `scipy.integrate.quad` .
|
| >>> import numpy as np
| >>> from scipy.stats import vonmises
| >>> res = vonmises(loc=2, kappa=1).expect(lambda x: np.exp(1j*x),
| ... complex_func=True)
| >>> res
| (-0.18576377217422957+0.40590124735052263j)
|
| >>> np.angle(res) # location of the (circular) distribution
| 2.0
|
| fit_loc_scale(self, data, *args)
| Estimate loc and scale parameters from data using 1st and 2nd moments.
|
| Parameters
| ----------
| data : array_like
| Data to fit.
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
|
| Returns
| -------
| Lhat : float
| Estimated location parameter for the data.
| Shat : float
| Estimated scale parameter for the data.
|
| isf(self, q, *args, **kwds)
| Inverse survival function (inverse of `sf`) at q of the given RV.
|
| Parameters
| ----------
| q : array_like
| upper tail probability
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| x : ndarray or scalar
| Quantile corresponding to the upper tail probability q.
|
| logcdf(self, x, *args, **kwds)
| Log of the cumulative distribution function at x of the given RV.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| logcdf : array_like
| Log of the cumulative distribution function evaluated at x
|
| logpdf(self, x, *args, **kwds)
| Log of the probability density function at x of the given RV.
|
| This uses a more numerically accurate calculation if available.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| logpdf : array_like
| Log of the probability density function evaluated at x
|
| logsf(self, x, *args, **kwds)
| Log of the survival function of the given RV.
|
| Returns the log of the "survival function," defined as (1 - `cdf`),
| evaluated at `x`.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| logsf : ndarray
| Log of the survival function evaluated at `x`.
|
| pdf(self, x, *args, **kwds)
| Probability density function at x of the given RV.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| pdf : ndarray
| Probability density function evaluated at x
|
| ppf(self, q, *args, **kwds)
| Percent point function (inverse of `cdf`) at q of the given RV.
|
| Parameters
| ----------
| q : array_like
| lower tail probability
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| x : array_like
| quantile corresponding to the lower tail probability q.
|
| sf(self, x, *args, **kwds)
| Survival function (1 - `cdf`) at x of the given RV.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| sf : array_like
| Survival function evaluated at x
|
| ----------------------------------------------------------------------
| Methods inherited from scipy.stats._distn_infrastructure.rv_generic:
|
| __call__(self, *args, **kwds)
| Freeze the distribution for the given arguments.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution. Should include all
| the non-optional arguments, may include ``loc`` and ``scale``.
|
| Returns
| -------
| rv_frozen : rv_frozen instance
| The frozen distribution.
|
| __setstate__(self, state)
|
| entropy(self, *args, **kwds)
| Differential entropy of the RV.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| Location parameter (default=0).
| scale : array_like, optional (continuous distributions only).
| Scale parameter (default=1).
|
| Notes
| -----
| Entropy is defined base `e`:
|
| >>> import numpy as np
| >>> from scipy.stats._distn_infrastructure import rv_discrete
| >>> drv = rv_discrete(values=((0, 1), (0.5, 0.5)))
| >>> np.allclose(drv.entropy(), np.log(2.0))
| True
|
| interval(self, confidence, *args, **kwds)
| Confidence interval with equal areas around the median.
|
| Parameters
| ----------
| confidence : array_like of float
| Probability that an rv will be drawn from the returned range.
| Each value should be in the range [0, 1].
| arg1, arg2, ... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| location parameter, Default is 0.
| scale : array_like, optional
| scale parameter, Default is 1.
|
| Returns
| -------
| a, b : ndarray of float
| end-points of range that contain ``100 * alpha %`` of the rv's
| possible values.
|
| Notes
| -----
| This is implemented as ``ppf([p_tail, 1-p_tail])``, where
| ``ppf`` is the inverse cumulative distribution function and
| ``p_tail = (1-confidence)/2``. Suppose ``[c, d]`` is the support of a
| discrete distribution; then ``ppf([0, 1]) == (c-1, d)``. Therefore,
| when ``confidence=1`` and the distribution is discrete, the left end
| of the interval will be beyond the support of the distribution.
| For discrete distributions, the interval will limit the probability
| in each tail to be less than or equal to ``p_tail`` (usually
| strictly less).
|
| mean(self, *args, **kwds)
| Mean of the distribution.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| mean : float
| the mean of the distribution
|
| median(self, *args, **kwds)
| Median of the distribution.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| Location parameter, Default is 0.
| scale : array_like, optional
| Scale parameter, Default is 1.
|
| Returns
| -------
| median : float
| The median of the distribution.
|
| See Also
| --------
| rv_discrete.ppf
| Inverse of the CDF
|
| moment(self, order, *args, **kwds)
| non-central moment of distribution of specified order.
|
| Parameters
| ----------
| order : int, order >= 1
| Order of moment.
| arg1, arg2, arg3,... : float
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| nnlf(self, theta, x)
| Negative loglikelihood function.
| Notes
| -----
| This is ``-sum(log pdf(x, theta), axis=0)`` where `theta` are the
| parameters (including loc and scale).
|
| rvs(self, *args, **kwds)
| Random variates of given type.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| Location parameter (default=0).
| scale : array_like, optional
| Scale parameter (default=1).
| size : int or tuple of ints, optional
| Defining number of random variates (default is 1).
| random_state : {None, int, `numpy.random.Generator`,
| `numpy.random.RandomState`}, optional
|
| If `random_state` is None (or `np.random`), the
| `numpy.random.RandomState` singleton is used.
| If `random_state` is an int, a new ``RandomState`` instance is
| used, seeded with `random_state`.
| If `random_state` is already a ``Generator`` or ``RandomState``
| instance, that instance is used.
|
| Returns
| -------
| rvs : ndarray or scalar
| Random variates of given `size`.
|
| stats(self, *args, **kwds)
| Some statistics of the given RV.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional (continuous RVs only)
| scale parameter (default=1)
| moments : str, optional
| composed of letters ['mvsk'] defining which moments to compute:
| 'm' = mean,
| 'v' = variance,
| 's' = (Fisher's) skew,
| 'k' = (Fisher's) kurtosis.
| (default is 'mv')
|
| Returns
| -------
| stats : sequence
| of requested moments.
|
| std(self, *args, **kwds)
| Standard deviation of the distribution.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| std : float
| standard deviation of the distribution
|
| support(self, *args, **kwargs)
| Support of the distribution.
|
| Parameters
| ----------
| arg1, arg2, ... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| location parameter, Default is 0.
| scale : array_like, optional
| scale parameter, Default is 1.
|
| Returns
| -------
| a, b : array_like
| end-points of the distribution's support.
|
| var(self, *args, **kwds)
| Variance of the distribution.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| var : float
| the variance of the distribution
|
| ----------------------------------------------------------------------
| Data descriptors inherited from scipy.stats._distn_infrastructure.rv_generic:
|
| random_state
| Get or set the generator object for generating random variates.
|
| If `random_state` is None (or `np.random`), the
| `numpy.random.RandomState` singleton is used.
| If `random_state` is an int, a new ``RandomState`` instance is used,
| seeded with `random_state`.
| If `random_state` is already a ``Generator`` or ``RandomState``
| instance, that instance is used.
help(qp.stats.lognorm)
Help on class lognorm in module qp.core.factory:
class lognorm(qp.parameterizations.base.Pdf_gen_wrap, scipy.stats._continuous_distns.lognorm_gen)
| lognorm(*args, **kwargs)
|
| A lognormal continuous random variable.
|
| %(before_notes)s
|
| Notes
| -----
| The probability density function for `lognorm` is:
|
| .. math::
|
| f(x, s) = \frac{1}{s x \sqrt{2\pi}}
| \exp\left(-\frac{\log^2(x)}{2s^2}\right)
|
| for :math:`x > 0`, :math:`s > 0`.
|
| `lognorm` takes ``s`` as a shape parameter for :math:`s`.
|
| %(after_notes)s
|
| Suppose a normally distributed random variable ``X`` has mean ``mu`` and
| standard deviation ``sigma``. Then ``Y = exp(X)`` is lognormally
| distributed with ``s = sigma`` and ``scale = exp(mu)``.
|
| %(example)s
|
| The logarithm of a log-normally distributed random variable is
| normally distributed:
|
| >>> import numpy as np
| >>> import matplotlib.pyplot as plt
| >>> from scipy import stats
| >>> fig, ax = plt.subplots(1, 1)
| >>> mu, sigma = 2, 0.5
| >>> X = stats.norm(loc=mu, scale=sigma)
| >>> Y = stats.lognorm(s=sigma, scale=np.exp(mu))
| >>> x = np.linspace(*X.interval(0.999))
| >>> y = Y.rvs(size=10000)
| >>> ax.plot(x, X.pdf(x), label='X (pdf)')
| >>> ax.hist(np.log(y), density=True, bins=x, label='log(Y) (histogram)')
| >>> ax.legend()
| >>> plt.show()
|
| Method resolution order:
| lognorm
| qp.parameterizations.base.Pdf_gen_wrap
| qp.parameterizations.base.Pdf_gen
| scipy.stats._continuous_distns.lognorm_gen
| scipy.stats._distn_infrastructure.rv_continuous
| scipy.stats._distn_infrastructure.rv_generic
| builtins.object
|
| Methods defined here:
|
| create_ensemble(data: 'Mapping', ancil: 'Optional[Mapping]' = None) -> 'Ensemble'
| Creates an Ensemble of distribution(s) in the given parameterization.
|
| Input data format:
| data = {'arg1': values, 'arg2': values ...} where 'arg1', 'arg2'... are the arguments for the parameterization.
| The length of the values should be the number of distributions being created in the Ensemble, with a minimum value of 1.
|
|
| Parameters
| ----------
| data : Mapping
| The dictionary of data for the distributions.
| ancil : Optional[Mapping], optional
| A dictionary of metadata for the distributions, where any arrays have the same length as the number of distributions, by default None
|
| Returns
| -------
| Ensemble
| An Ensemble object containing all of the given distributions.
|
| Examples
| --------
|
| To create an Ensemble with two Gaussian distributions and their associated ids:
|
| >>> import qp
| >>> data = {'loc': np.array([[0.45],[0.55]]) , 'scale': np.array([[0.2],[0.15]])}
| >>> ancil = {'ids': [20,25]}
| >>> ens = qp.stats.norm.create_ensemble(data,ancil)
| >>> ens.metadata
| {'pdf_name': array([b'norm'], dtype='|S4'), 'pdf_version': array([0])}
|
| freeze = _my_freeze(self, *args, **kwds)
|
| ----------------------------------------------------------------------
| Data and other attributes defined here:
|
| name = 'lognorm'
|
| version = 0
|
| ----------------------------------------------------------------------
| Methods inherited from qp.parameterizations.base.Pdf_gen_wrap:
|
| __init__(self, *args, **kwargs)
| C'tor
|
| ----------------------------------------------------------------------
| Class methods inherited from qp.parameterizations.base.Pdf_gen_wrap:
|
| add_mappings() from builtins.type
| Add this classes mappings to the conversion dictionary
|
| create(**kwds) from builtins.type
| Create and return a `scipy.stats.rv_frozen` object using the
| keyword arguments provided
|
| create_gen(**kwds) from builtins.type
| Create and return a `scipy.stats.rv_continuous` object using the
| keyword arguments provided
|
| get_allocation_kwds(npdf, **kwargs) from builtins.type
| Return kwds necessary to create 'empty' hdf5 file with npdf entries
| for iterative writeout
|
| ----------------------------------------------------------------------
| Class methods inherited from qp.parameterizations.base.Pdf_gen:
|
| add_method_dicts() from builtins.type
| Add empty method dicts
|
| creation_method(method=None) from builtins.type
| Return the method used to create a PDF of this type
|
| extraction_method(method=None) from builtins.type
| Return the method used to extract data to create a PDF of this type
|
| plot(pdf, **kwargs) from builtins.type
| Plot the pdf as a curve
|
| plot_native(pdf, **kwargs) from builtins.type
| Plot the PDF in a way that is particular to this type of distribution
|
| This defaults to plotting it as a curve, but this can be overwritten
|
| print_method_maps(stream=<ipykernel.iostream.OutStream object at 0x7f54f45ac340>) from builtins.type
| Print the maps showing the methods
|
| reader_method(version=None) from builtins.type
| Return the method used to convert data read from a file PDF of this type
|
| ----------------------------------------------------------------------
| Readonly properties inherited from qp.parameterizations.base.Pdf_gen:
|
| metadata
| Return the metadata for this set of PDFs
|
| objdata
| Return the object data for this set of PDFs
|
| ----------------------------------------------------------------------
| Data descriptors inherited from qp.parameterizations.base.Pdf_gen:
|
| __dict__
| dictionary for instance variables (if defined)
|
| __weakref__
| list of weak references to the object (if defined)
|
| ----------------------------------------------------------------------
| Methods inherited from scipy.stats._continuous_distns.lognorm_gen:
|
| fit(self, data, *args, **kwds)
| Return estimates of shape (if applicable), location, and scale
| parameters from data. The default estimation method is Maximum
| Likelihood Estimation (MLE), but Method of Moments (MM)
| is also available.
|
| Starting estimates for the fit are given by input arguments;
| for any arguments not provided with starting estimates,
| ``self._fitstart(data)`` is called to generate such.
|
| One can hold some parameters fixed to specific values by passing in
| keyword arguments ``f0``, ``f1``, ..., ``fn`` (for shape parameters)
| and ``floc`` and ``fscale`` (for location and scale parameters,
| respectively).
|
| Parameters
| ----------
| data : array_like or `CensoredData` instance
| Data to use in estimating the distribution parameters.
| arg1, arg2, arg3,... : floats, optional
| Starting value(s) for any shape-characterizing arguments (those not
| provided will be determined by a call to ``_fitstart(data)``).
| No default value.
| **kwds : floats, optional
| - `loc`: initial guess of the distribution's location parameter.
| - `scale`: initial guess of the distribution's scale parameter.
|
| Special keyword arguments are recognized as holding certain
| parameters fixed:
|
| - f0...fn : hold respective shape parameters fixed.
| Alternatively, shape parameters to fix can be specified by name.
| For example, if ``self.shapes == "a, b"``, ``fa`` and ``fix_a``
| are equivalent to ``f0``, and ``fb`` and ``fix_b`` are
| equivalent to ``f1``.
|
| - floc : hold location parameter fixed to specified value.
|
| - fscale : hold scale parameter fixed to specified value.
|
| - optimizer : The optimizer to use. The optimizer must take
| ``func`` and starting position as the first two arguments,
| plus ``args`` (for extra arguments to pass to the
| function to be optimized) and ``disp``.
| The ``fit`` method calls the optimizer with ``disp=0`` to suppress output.
| The optimizer must return the estimated parameters.
|
| - method : The method to use. The default is "MLE" (Maximum
| Likelihood Estimate); "MM" (Method of Moments)
| is also available.
|
| Raises
| ------
| TypeError, ValueError
| If an input is invalid
| `~scipy.stats.FitError`
| If fitting fails or the fit produced would be invalid
|
| Returns
| -------
| parameter_tuple : tuple of floats
| Estimates for any shape parameters (if applicable), followed by
| those for location and scale. For most random variables, shape
| statistics will be returned, but there are exceptions (e.g.
| ``norm``).
|
| Notes
| -----
| With ``method="MLE"`` (default), the fit is computed by minimizing
| the negative log-likelihood function. A large, finite penalty
| (rather than infinite negative log-likelihood) is applied for
| observations beyond the support of the distribution.
|
| With ``method="MM"``, the fit is computed by minimizing the L2 norm
| of the relative errors between the first *k* raw (about zero) data
| moments and the corresponding distribution moments, where *k* is the
| number of non-fixed parameters.
| More precisely, the objective function is::
|
| (((data_moments - dist_moments)
| / np.maximum(np.abs(data_moments), 1e-8))**2).sum()
|
| where the constant ``1e-8`` avoids division by zero in case of
| vanishing data moments. Typically, this error norm can be reduced to
| zero.
| Note that the standard method of moments can produce parameters for
| which some data are outside the support of the fitted distribution;
| this implementation does nothing to prevent this.
|
| For either method,
| the returned answer is not guaranteed to be globally optimal; it
| may only be locally optimal, or the optimization may fail altogether.
| If the data contain any of ``np.nan``, ``np.inf``, or ``-np.inf``,
| the `fit` method will raise a ``RuntimeError``.
|
| When passing a ``CensoredData`` instance to ``data``, the log-likelihood
| function is defined as:
|
| .. math::
|
| l(\pmb{\theta}; k) & = \sum
| \log(f(k_u; \pmb{\theta}))
| + \sum
| \log(F(k_l; \pmb{\theta})) \\
| & + \sum
| \log(1 - F(k_r; \pmb{\theta})) \\
| & + \sum
| \log(F(k_{\text{high}, i}; \pmb{\theta})
| - F(k_{\text{low}, i}; \pmb{\theta}))
|
| where :math:`f` and :math:`F` are the pdf and cdf, respectively, of the
| function being fitted, :math:`\pmb{\theta}` is the parameter vector,
| :math:`u` are the indices of uncensored observations,
| :math:`l` are the indices of left-censored observations,
| :math:`r` are the indices of right-censored observations,
| subscripts "low"/"high" denote endpoints of interval-censored observations, and
| :math:`i` are the indices of interval-censored observations.
|
| When `method='MLE'` and
| the location parameter is fixed by using the `floc` argument,
| this function uses explicit formulas for the maximum likelihood
| estimation of the log-normal shape and scale parameters, so the
| `optimizer`, `loc` and `scale` keyword arguments are ignored.
| If the location is free, a likelihood maximum is found by
| setting its partial derivative wrt to location to 0, and
| solving by substituting the analytical expressions of shape
| and scale (or provided parameters).
| See, e.g., equation 3.1 in
| A. Clifford Cohen & Betty Jones Whitten (1980)
| Estimation in the Three-Parameter Lognormal Distribution,
| Journal of the American Statistical Association, 75:370, 399-404
| https://doi.org/10.2307/2287466
|
|
| Examples
| --------
|
| Generate some data to fit: draw random variates from the `beta`
| distribution
|
| >>> import numpy as np
| >>> from scipy.stats import beta
| >>> a, b = 1., 2.
| >>> rng = np.random.default_rng(172786373191770012695001057628748821561)
| >>> x = beta.rvs(a, b, size=1000, random_state=rng)
|
| Now we can fit all four parameters (``a``, ``b``, ``loc`` and
| ``scale``):
|
| >>> a1, b1, loc1, scale1 = beta.fit(x)
| >>> a1, b1, loc1, scale1
| (1.0198945204435628, 1.9484708982737828, 4.372241314917588e-05, 0.9979078845964814)
|
| The fit can be done also using a custom optimizer:
|
| >>> from scipy.optimize import minimize
| >>> def custom_optimizer(func, x0, args=(), disp=0):
| ... res = minimize(func, x0, args, method="slsqp", options={"disp": disp})
| ... if res.success:
| ... return res.x
| ... raise RuntimeError('optimization routine failed')
| >>> a1, b1, loc1, scale1 = beta.fit(x, method="MLE", optimizer=custom_optimizer)
| >>> a1, b1, loc1, scale1
| (1.0198821087258905, 1.948484145914738, 4.3705304486881485e-05, 0.9979104663953395)
|
| We can also use some prior knowledge about the dataset: let's keep
| ``loc`` and ``scale`` fixed:
|
| >>> a1, b1, loc1, scale1 = beta.fit(x, floc=0, fscale=1)
| >>> loc1, scale1
| (0, 1)
|
| We can also keep shape parameters fixed by using ``f``-keywords. To
| keep the zero-th shape parameter ``a`` equal 1, use ``f0=1`` or,
| equivalently, ``fa=1``:
|
| >>> a1, b1, loc1, scale1 = beta.fit(x, fa=1, floc=0, fscale=1)
| >>> a1
| 1
|
| Not all distributions return estimates for the shape parameters.
| ``norm`` for example just returns estimates for location and scale:
|
| >>> from scipy.stats import norm
| >>> x = norm.rvs(a, b, size=1000, random_state=123)
| >>> loc1, scale1 = norm.fit(x)
| >>> loc1, scale1
| (0.92087172783841631, 2.0015750750324668)
|
| ----------------------------------------------------------------------
| Methods inherited from scipy.stats._distn_infrastructure.rv_continuous:
|
| __getstate__(self)
|
| cdf(self, x, *args, **kwds)
| Cumulative distribution function of the given RV.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| cdf : ndarray
| Cumulative distribution function evaluated at `x`
|
| expect(self, func=None, args=(), loc=0, scale=1, lb=None, ub=None, conditional=False, **kwds)
| Calculate expected value of a function with respect to the
| distribution by numerical integration.
|
| The expected value of a function ``f(x)`` with respect to a
| distribution ``dist`` is defined as::
|
| ub
| E[f(x)] = Integral(f(x) * dist.pdf(x)),
| lb
|
| where ``ub`` and ``lb`` are arguments and ``x`` has the ``dist.pdf(x)``
| distribution. If the bounds ``lb`` and ``ub`` correspond to the
| support of the distribution, e.g. ``[-inf, inf]`` in the default
| case, then the integral is the unrestricted expectation of ``f(x)``.
| Also, the function ``f(x)`` may be defined such that ``f(x)`` is ``0``
| outside a finite interval in which case the expectation is
| calculated within the finite range ``[lb, ub]``.
|
| Parameters
| ----------
| func : callable, optional
| Function for which integral is calculated. Takes only one argument.
| The default is the identity mapping f(x) = x.
| args : tuple, optional
| Shape parameters of the distribution.
| loc : float, optional
| Location parameter (default=0).
| scale : float, optional
| Scale parameter (default=1).
| lb, ub : scalar, optional
| Lower and upper bound for integration. Default is set to the
| support of the distribution.
| conditional : bool, optional
| If True, the integral is corrected by the conditional probability
| of the integration interval. The return value is the expectation
| of the function, conditional on being in the given interval.
| Default is False.
|
| Additional keyword arguments are passed to the integration routine.
|
| Returns
| -------
| expect : float
| The calculated expected value.
|
| Notes
| -----
| The integration behavior of this function is inherited from
| `scipy.integrate.quad`. Neither this function nor
| `scipy.integrate.quad` can verify whether the integral exists or is
| finite. For example ``cauchy(0).mean()`` returns ``np.nan`` and
| ``cauchy(0).expect()`` returns ``0.0``.
|
| Likewise, the accuracy of results is not verified by the function.
| `scipy.integrate.quad` is typically reliable for integrals that are
| numerically favorable, but it is not guaranteed to converge
| to a correct value for all possible intervals and integrands. This
| function is provided for convenience; for critical applications,
| check results against other integration methods.
|
| The function is not vectorized.
|
| Examples
| --------
|
| To understand the effect of the bounds of integration consider
|
| >>> from scipy.stats import expon
| >>> expon(1).expect(lambda x: 1, lb=0.0, ub=2.0)
| 0.6321205588285578
|
| This is close to
|
| >>> expon(1).cdf(2.0) - expon(1).cdf(0.0)
| 0.6321205588285577
|
| If ``conditional=True``
|
| >>> expon(1).expect(lambda x: 1, lb=0.0, ub=2.0, conditional=True)
| 1.0000000000000002
|
| The slight deviation from 1 is due to numerical integration.
|
| The integrand can be treated as a complex-valued function
| by passing ``complex_func=True`` to `scipy.integrate.quad` .
|
| >>> import numpy as np
| >>> from scipy.stats import vonmises
| >>> res = vonmises(loc=2, kappa=1).expect(lambda x: np.exp(1j*x),
| ... complex_func=True)
| >>> res
| (-0.18576377217422957+0.40590124735052263j)
|
| >>> np.angle(res) # location of the (circular) distribution
| 2.0
|
| fit_loc_scale(self, data, *args)
| Estimate loc and scale parameters from data using 1st and 2nd moments.
|
| Parameters
| ----------
| data : array_like
| Data to fit.
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
|
| Returns
| -------
| Lhat : float
| Estimated location parameter for the data.
| Shat : float
| Estimated scale parameter for the data.
|
| isf(self, q, *args, **kwds)
| Inverse survival function (inverse of `sf`) at q of the given RV.
|
| Parameters
| ----------
| q : array_like
| upper tail probability
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| x : ndarray or scalar
| Quantile corresponding to the upper tail probability q.
|
| logcdf(self, x, *args, **kwds)
| Log of the cumulative distribution function at x of the given RV.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| logcdf : array_like
| Log of the cumulative distribution function evaluated at x
|
| logpdf(self, x, *args, **kwds)
| Log of the probability density function at x of the given RV.
|
| This uses a more numerically accurate calculation if available.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| logpdf : array_like
| Log of the probability density function evaluated at x
|
| logsf(self, x, *args, **kwds)
| Log of the survival function of the given RV.
|
| Returns the log of the "survival function," defined as (1 - `cdf`),
| evaluated at `x`.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| logsf : ndarray
| Log of the survival function evaluated at `x`.
|
| pdf(self, x, *args, **kwds)
| Probability density function at x of the given RV.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| pdf : ndarray
| Probability density function evaluated at x
|
| ppf(self, q, *args, **kwds)
| Percent point function (inverse of `cdf`) at q of the given RV.
|
| Parameters
| ----------
| q : array_like
| lower tail probability
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| x : array_like
| quantile corresponding to the lower tail probability q.
|
| sf(self, x, *args, **kwds)
| Survival function (1 - `cdf`) at x of the given RV.
|
| Parameters
| ----------
| x : array_like
| quantiles
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| sf : array_like
| Survival function evaluated at x
|
| ----------------------------------------------------------------------
| Methods inherited from scipy.stats._distn_infrastructure.rv_generic:
|
| __call__(self, *args, **kwds)
| Freeze the distribution for the given arguments.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution. Should include all
| the non-optional arguments, may include ``loc`` and ``scale``.
|
| Returns
| -------
| rv_frozen : rv_frozen instance
| The frozen distribution.
|
| __setstate__(self, state)
|
| entropy(self, *args, **kwds)
| Differential entropy of the RV.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| Location parameter (default=0).
| scale : array_like, optional (continuous distributions only).
| Scale parameter (default=1).
|
| Notes
| -----
| Entropy is defined base `e`:
|
| >>> import numpy as np
| >>> from scipy.stats._distn_infrastructure import rv_discrete
| >>> drv = rv_discrete(values=((0, 1), (0.5, 0.5)))
| >>> np.allclose(drv.entropy(), np.log(2.0))
| True
|
| interval(self, confidence, *args, **kwds)
| Confidence interval with equal areas around the median.
|
| Parameters
| ----------
| confidence : array_like of float
| Probability that an rv will be drawn from the returned range.
| Each value should be in the range [0, 1].
| arg1, arg2, ... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| location parameter, Default is 0.
| scale : array_like, optional
| scale parameter, Default is 1.
|
| Returns
| -------
| a, b : ndarray of float
| end-points of range that contain ``100 * alpha %`` of the rv's
| possible values.
|
| Notes
| -----
| This is implemented as ``ppf([p_tail, 1-p_tail])``, where
| ``ppf`` is the inverse cumulative distribution function and
| ``p_tail = (1-confidence)/2``. Suppose ``[c, d]`` is the support of a
| discrete distribution; then ``ppf([0, 1]) == (c-1, d)``. Therefore,
| when ``confidence=1`` and the distribution is discrete, the left end
| of the interval will be beyond the support of the distribution.
| For discrete distributions, the interval will limit the probability
| in each tail to be less than or equal to ``p_tail`` (usually
| strictly less).
|
| mean(self, *args, **kwds)
| Mean of the distribution.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| mean : float
| the mean of the distribution
|
| median(self, *args, **kwds)
| Median of the distribution.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| Location parameter, Default is 0.
| scale : array_like, optional
| Scale parameter, Default is 1.
|
| Returns
| -------
| median : float
| The median of the distribution.
|
| See Also
| --------
| rv_discrete.ppf
| Inverse of the CDF
|
| moment(self, order, *args, **kwds)
| non-central moment of distribution of specified order.
|
| Parameters
| ----------
| order : int, order >= 1
| Order of moment.
| arg1, arg2, arg3,... : float
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| nnlf(self, theta, x)
| Negative loglikelihood function.
| Notes
| -----
| This is ``-sum(log pdf(x, theta), axis=0)`` where `theta` are the
| parameters (including loc and scale).
|
| rvs(self, *args, **kwds)
| Random variates of given type.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| Location parameter (default=0).
| scale : array_like, optional
| Scale parameter (default=1).
| size : int or tuple of ints, optional
| Defining number of random variates (default is 1).
| random_state : {None, int, `numpy.random.Generator`,
| `numpy.random.RandomState`}, optional
|
| If `random_state` is None (or `np.random`), the
| `numpy.random.RandomState` singleton is used.
| If `random_state` is an int, a new ``RandomState`` instance is
| used, seeded with `random_state`.
| If `random_state` is already a ``Generator`` or ``RandomState``
| instance, that instance is used.
|
| Returns
| -------
| rvs : ndarray or scalar
| Random variates of given `size`.
|
| stats(self, *args, **kwds)
| Some statistics of the given RV.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional (continuous RVs only)
| scale parameter (default=1)
| moments : str, optional
| composed of letters ['mvsk'] defining which moments to compute:
| 'm' = mean,
| 'v' = variance,
| 's' = (Fisher's) skew,
| 'k' = (Fisher's) kurtosis.
| (default is 'mv')
|
| Returns
| -------
| stats : sequence
| of requested moments.
|
| std(self, *args, **kwds)
| Standard deviation of the distribution.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| std : float
| standard deviation of the distribution
|
| support(self, *args, **kwargs)
| Support of the distribution.
|
| Parameters
| ----------
| arg1, arg2, ... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information).
| loc : array_like, optional
| location parameter, Default is 0.
| scale : array_like, optional
| scale parameter, Default is 1.
|
| Returns
| -------
| a, b : array_like
| end-points of the distribution's support.
|
| var(self, *args, **kwds)
| Variance of the distribution.
|
| Parameters
| ----------
| arg1, arg2, arg3,... : array_like
| The shape parameter(s) for the distribution (see docstring of the
| instance object for more information)
| loc : array_like, optional
| location parameter (default=0)
| scale : array_like, optional
| scale parameter (default=1)
|
| Returns
| -------
| var : float
| the variance of the distribution
|
| ----------------------------------------------------------------------
| Data descriptors inherited from scipy.stats._distn_infrastructure.rv_generic:
|
| random_state
| Get or set the generator object for generating random variates.
|
| If `random_state` is None (or `np.random`), the
| `numpy.random.RandomState` singleton is used.
| If `random_state` is an int, a new ``RandomState`` instance is used,
| seeded with `random_state`.
| If `random_state` is already a ``Generator`` or ``RandomState``
| instance, that instance is used.
Native plotting
If you have a single distribution you can plot it, the qp.plotting.plot_native function will find a nice way to represent the data used to construct the distribution.
loc1 = np.array([[0]])
scale1 = np.array([[1]])
norm_dist1 = qp.stats.norm(loc=loc1, scale=scale1)
fig, axes = qp.plotting.plot_native(norm_dist1, xlim=(-5., 5.))
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[9], line 3
1 loc1 = np.array([[0]])
2 scale1 = np.array([[1]])
----> 3 norm_dist1 = qp.stats.norm(loc=loc1, scale=scale1)
4 fig, axes = qp.plotting.plot_native(norm_dist1, xlim=(-5., 5.))
File ~/checkouts/readthedocs.org/user_builds/qp/envs/latest/lib/python3.10/site-packages/qp/parameterizations/base.py:409, in Pdf_gen_wrap.__init__(self, *args, **kwargs)
407 # pylint: disable=no-member,protected-access
408 super().__init__(*args, **kwargs)
--> 409 self._other_init(*args, **kwargs)
TypeError: rv_continuous.__init__() got an unexpected keyword argument 'loc'
# fig, axes = qp.stats.norm.plot_native(norm_dist1, xlim=(-5., 5.))
qp histogram (piecewise constant) parameterization
This represents a set of distributions made by interpolating a set of histograms with shared binning. To construct this you need to give the bin edges (shape=(N)) and the bin values (shape=(npdf, N-1)).
Note that the native visual representation is different from the Normal distribution.
# Convert to a histogram by computing the bin values by taking the intergral of the CDF
xvals = np.linspace(-5, 5, 11)
cdf = norm_dist1.cdf(xvals)
bin_vals = cdf[:,1:] - cdf[:,0:-1]
# Construct histogram PDF using the bin edges and the bin values
hist_dist = qp.hist(bins=xvals, pdfs=bin_vals)
yvals = hist_dist.pdf(xvals)
# Construct a single PDF for plotting
hist_dist1 = qp.hist(bins=xvals, pdfs=np.atleast_2d(bin_vals[0]))
fig, axes = qp.plotting.plot_native(hist_dist1, xlim=(-5., 5.))
leg = fig.legend()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 3
1 # Convert to a histogram by computing the bin values by taking the intergral of the CDF
2 xvals = np.linspace(-5, 5, 11)
----> 3 cdf = norm_dist1.cdf(xvals)
4 bin_vals = cdf[:,1:] - cdf[:,0:-1]
5 # Construct histogram PDF using the bin edges and the bin values
NameError: name 'norm_dist1' is not defined
What if you want to evaluate a vector of input values, where each input value is different for each PDF? In that case you need the shape of the vector of input value to match the implicit shape of the PDFs, which in this case is (2,1)
xvals_x = np.array([[-1.], [1.]])
yvals_x = hist_dist.pdf(xvals_x)
print ("For an input vector of shape %s the output shape is %s" % (xvals_x.shape, yvals_x.shape))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 2
1 xvals_x = np.array([[-1.], [1.]])
----> 2 yvals_x = hist_dist.pdf(xvals_x)
3 print ("For an input vector of shape %s the output shape is %s" % (xvals_x.shape, yvals_x.shape))
NameError: name 'hist_dist' is not defined
qp quantile parameterization
This represents a set of distributions made by interpolating the locations at which various distributions reach a given set of quantiles. To construct this you need to give the quantiles edges (shape=(N)) and the location values (shape=(npdf, N)).
Note that the native visual representation is different.
# Define the quantile values to compute the locations for
quants = np.linspace(0.01, 0.99, 7)
# Compute the corresponding locations
locs = norm_dist1.ppf(quants)
# Construct the distribution using the quantile value and locations
quant_dist = qp.quant(quants=quants, locs=locs)
quant_vals = quant_dist.pdf(xvals)
print("The input and output shapes are:", xvals.shape, quant_vals.shape)
# Construct a single PDF for plotting
quant_dist1 = qp.quant(quants=np.atleast_1d(quants), locs=np.atleast_2d(locs[0]))
fig, axes = qp.plotting.plot_native(quant_dist1, xlim=(-5., 5.), label="quantiles")
leg = fig.legend()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 4
2 quants = np.linspace(0.01, 0.99, 7)
3 # Compute the corresponding locations
----> 4 locs = norm_dist1.ppf(quants)
5 # Construct the distribution using the quantile value and locations
6 quant_dist = qp.quant(quants=quants, locs=locs)
NameError: name 'norm_dist1' is not defined
print(quants)
print(quant_dist.dist.quants)
[0.01 0.17333333 0.33666667 0.5 0.66333333 0.82666667
0.99 ]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 2
1 print(quants)
----> 2 print(quant_dist.dist.quants)
NameError: name 'quant_dist' is not defined
qp interpolated parameterization
This represents a set of distributions made by interpolating a set of x and y values. To construct this you need to give the x and y values (both of shape=(npdf, N))
Note that the native visual representation is pretty similar to the original one for the Gaussian.
# Define the x-grid locations
xvals = np.linspace(-5, 5, 11)
# Compute the corresponding y values
yvals = norm_dist1.pdf(xvals)
# Construct the PDFs using the x grid and y values
interp_dist = qp.interp(xvals=xvals, yvals=yvals)
interp_vals = interp_dist.pdf(xvals)
print("The input and output shapes are:", xvals.shape, interp_vals.shape)
# Construct a single PDF for plotting
interp_dist1 = qp.interp(xvals=xvals, yvals=np.atleast_2d(yvals[0]))
fig, axes = qp.plotting.plot_native(interp_dist1, xlim=(-5., 5.), label="interpolated")
leg = fig.legend()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 4
2 xvals = np.linspace(-5, 5, 11)
3 # Compute the corresponding y values
----> 4 yvals = norm_dist1.pdf(xvals)
5 # Construct the PDFs using the x grid and y values
6 interp_dist = qp.interp(xvals=xvals, yvals=yvals)
NameError: name 'norm_dist1' is not defined
qp spline parameterization constructed from kernel density estimate (samples) parameterization
This represents a set of distributions made by producing a kernel density estimate from a set of samples.
To construct this you need to give the samples edges (shape=(npdf, Nsamples)).
Note again that the the native visual represenation is different.
# Take 100 random samples from each of 2 PDFs
samples = norm_dist1.rvs(size=(2, 1000))
# Define points at which to evaluate the kernal density estimate (KDE)
xvals_kde = np.linspace(-5., 5., 51)
# Use a utility function to construct the KDE, sample it, and they construct a spline
kde_dist = qp.spline_from_samples(xvals=xvals_kde, samples=samples)
kde_vals = kde_dist.pdf(xvals_kde)
print("The input and output shapes are:", xvals.shape, kde_vals.shape)
# Construct a single PDF for plotting
kde_dist1 = qp.spline_from_samples(xvals=xvals_kde, samples=np.atleast_2d(samples[0]))
fig, axes = qp.plotting.plot_native(kde_dist1, xlim=(-5., 5.), label="kde")
leg = fig.legend()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 2
1 # Take 100 random samples from each of 2 PDFs
----> 2 samples = norm_dist1.rvs(size=(2, 1000))
3 # Define points at which to evaluate the kernal density estimate (KDE)
4 xvals_kde = np.linspace(-5., 5., 51)
NameError: name 'norm_dist1' is not defined
qp spline parameterization
This represents a set of distributions made building a set of splines. Though the parameterization is defined by the spline knots, you can construct this from x and y values (both of shape=(npdf, N)).
Note that the native visual representation is pretty similar to the original one for the Gaussian.
Note also that the spline knots are stored.
# To make a spline you need the spline knots, you can get those from the xval, yval values
splx, sply, spln = qp.spline_gen.build_normed_splines(np.expand_dims(xvals,0), yvals)
spline_dist_orig = qp.spline(splx=splx, sply=sply, spln=spln)
# Or we can do these two steps together using one function
spline_dist = qp.spline_from_xy(xvals=np.expand_dims(xvals,0), yvals=yvals)
spline_vals = spline_dist.pdf(xvals)
print("The input and output shapes are:", xvals.shape, spline_vals.shape)
print("Spline knots", spline_dist.dist.splx, spline_dist.dist.sply, spline_dist.dist.spln)
# Construct a single PDF for plotting
spline_dist1 = qp.spline_from_xy(xvals=np.atleast_2d(xvals), yvals=np.atleast_2d(yvals))
print(spline_dist1.dist.splx.shape)
fig, axes = qp.plotting.plot_native(spline_dist1, xlim=(-5., 5.), label="spline")
leg = fig.legend()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 2
1 # To make a spline you need the spline knots, you can get those from the xval, yval values
----> 2 splx, sply, spln = qp.spline_gen.build_normed_splines(np.expand_dims(xvals,0), yvals)
3 spline_dist_orig = qp.spline(splx=splx, sply=sply, spln=spln)
4 # Or we can do these two steps together using one function
NameError: name 'yvals' is not defined
Overplotting
You can visually compare the represenations by plotting them all on the same figure.
fig, axes = qp.plotting.plot_native(norm_dist1, xlim=(-5., 5.), label="norm")
qp.plotting.plot_native(hist_dist1, axes=axes)
qp.plotting.plot_native(quant_dist1, axes=axes)
qp.plotting.plot_native(interp_dist1, axes=axes, label="interp")
# qp.plotting.plot_native(kde_dist1, axes=axes)
# qp.plotting.plot_native(spline_dist1, axes=axes, label="spline")
leg = fig.legend()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 fig, axes = qp.plotting.plot_native(norm_dist1, xlim=(-5., 5.), label="norm")
2 qp.plotting.plot_native(hist_dist1, axes=axes)
3 qp.plotting.plot_native(quant_dist1, axes=axes)
NameError: name 'norm_dist1' is not defined
The qp.Ensemble Class
This is the basic element of qp - an object representing a set of probability density functions. This class is stored in the module ensemble.py.
To create a qp.Ensemble you need to specify the class used to represent the PDFs, and provide that data for the specific set of PDFs.
Ensembles of distributions
qp no longer distinguishes between distributions and ensembles thereof – a single distribution is just a special case of an ensemble with only one member, which takes advantage of computational efficiencies in scipy.
The shape of the array returned by a call to the pdf function of a distribution depends on the shape of the parameters and evaluate points.
For distributions that take multiple input arrays, qp uses te convention that the rows are the individual distributions and the columns are the values of the parameters defining the distributions under a known parameterization.
# This is a trivial extension, with the number of pdfs as a member of the `scipy.stats.norm_gen` distribution.
loc = np.array([[0],[1]])
scale = np.array([[1],[1]])
norm_dist = qp.stats.norm(loc=loc, scale=scale)
xvals = np.linspace(-5, 5, 51)
yvals = norm_dist.pdf(xvals)
print("This object represents %i pdfs" % norm_dist.npdf)
print("The input and output shapes are:", xvals.shape, yvals.shape)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[19], line 4
2 loc = np.array([[0],[1]])
3 scale = np.array([[1],[1]])
----> 4 norm_dist = qp.stats.norm(loc=loc, scale=scale)
5 xvals = np.linspace(-5, 5, 51)
6 yvals = norm_dist.pdf(xvals)
File ~/checkouts/readthedocs.org/user_builds/qp/envs/latest/lib/python3.10/site-packages/qp/parameterizations/base.py:409, in Pdf_gen_wrap.__init__(self, *args, **kwargs)
407 # pylint: disable=no-member,protected-access
408 super().__init__(*args, **kwargs)
--> 409 self._other_init(*args, **kwargs)
TypeError: rv_continuous.__init__() got an unexpected keyword argument 'loc'
print ("For an input vector of shape %s the output shape is %s" % (xvals.shape, yvals.shape))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 1
----> 1 print ("For an input vector of shape %s the output shape is %s" % (xvals.shape, yvals.shape))
NameError: name 'yvals' is not defined
# In this case we return an array were the rows are the evaluation points and the columns the different PDFs
#vector_pdf = qp.stats.norm(loc=[0., 1., 2], scale=1.)
#vector_pdf.pdf([[0.], [0.5]])
# This is the same, except we use `numpy.expand_dims` to shape the input array of evaluation points
# vector_pdf = qp.stats.norm(loc=[0., 1., 2], scale=1.)
# vector_pdf.pdf(np.expand_dims(np.array([0., 0.5]), -1))
# In this case we return an array were the rows are pdfs and the columns the evaluation points
#vector_pdf = qp.stats.norm(loc=[[0.], [1.], [2]], scale=1.)
#vector_pdf.pdf([0., 0.5])
# This is the same, except we use `numpy.expand_dims` to shape the input array of pdf parameters
# vector_pdf = qp.stats.norm(loc=np.expand_dims([0., 1., 2], -1), scale=1.)
# vector_pdf.pdf([0., 0.5])
Here we will create 100 Gaussians with means distributed between -1 and 1, and widths distributed between 0.9 and 1.1.
locs = 2* (np.random.uniform(size=(100,1))-0.5)
scales = 1 + 0.2*(np.random.uniform(size=(100,1))-0.5)
ens_n = qp.Ensemble(qp.stats.norm, data=dict(loc=locs, scale=scales))
Using the ensemble
All of the methods of the distributions (pdf, cdf etc.) work the same way for an ensemble as for underlying classes.
To isolate a single distribution in the ensemble, use the square brackets operator [].
vals_n = ens_n.pdf(xvals)
print("The shapes are: ", xvals.shape, vals_n.shape)
fig, axes = qp.plotting.plot_native(ens_n[15], xlim=(-5.,5.))
The shapes are: (11,) (100, 11)
Converting the ensemble
The qp.Ensemble.convert_to function lets you convert ensembles to other representations. To do this you have to provide the original ensemble, the class you want to convert to, and any some keyword arguments to specify details about how to convert to the new class, here are some examples.
bins = np.linspace(-5, 5, 11)
quants = np.linspace(0.01, 0.99, 7)
print("Making hist")
ens_h = ens_n.convert_to(qp.hist_gen, bins=bins)
print("Making interp")
ens_i = ens_n.convert_to(qp.interp_gen, xvals=bins)
print("Making spline")
ens_s = ens_n.convert_to(qp.spline_gen, xvals=bins, method="xy")
#print("Making spline from samples")
#ens_s = ens_n.convert_to(qp.spline_gen, xvals=bins, samples=1000, method="samples")
print("Making quants")
ens_q = ens_n.convert_to(qp.quant_gen, quants=quants)
print("Making mixmod")
ens_m = ens_n.convert_to(qp.mixmod_gen, samples=1000, ncomps=3)
#print("Making flexcode")
#ens_f = ens_n.convert_to(qp.flex_gen, grid=bins, basis_system='cosine')
Making hist
Making interp
Making spline
Making quants
Making mixmod
The qp.convert function also works the more or less the same way, but with slightly different syntax, where you can use the name of the class instead of the class object.
print("Making hist")
ens_h2 = qp.convert(ens_n, "hist", bins=bins)
print("Making interp")
ens_i2 = qp.convert(ens_n, "interp", xvals=bins)
print("Making spline")
ens_s2 = qp.convert(ens_n, "spline", xvals=bins, method="xy")
print("Making quants")
ens_q2 = qp.convert(ens_n, "quant", quants=quants)
print("Making mixmod")
ens_m2 = qp.convert(ens_n, "mixmod", samples=1000, ncomps=3)
Making hist
Making interp
Making spline
Making quants
Making mixmod
Comparing Parametrizations
qp supports quantitative comparisons between different distributions, across parametrizations.
Qualitative Comparisons: Plotting
Let’s visualize the PDF object in order to original and the other representaions. The solid, black line shows the true PDF evaluated between the bounds. The green rugplot shows the locations of the 1000 samples we took. The vertical, dotted, blue lines show the percentiles we asked for, and the hotizontal, dotted, red lines show the 10 equally spaced bins we asked for. Note that the quantiles refer to the probability distribution between the bounds, because we are not able to integrate numerically over an infinite range. Interpolations of each parametrization are given as dashed lines in their corresponding colors. Note that the interpolations of the quantile and histogram parametrizations are so close to each other that the difference is almost imperceptible!
fig, axes = qp.plotting.plot_native(ens_n[15], xlim=(-5.,5.))
qp.plotting.plot_native(ens_h[15], axes=axes)
qp.plotting.plot_native(ens_q[15], axes=axes, label='quantile')
qp.plotting.plot_native(ens_i[15], axes=axes, label='interp')
# qp.plotting.plot_native(ens_s[15], axes=axes, label='spline')
qp.plotting.plot_native(ens_m[15], axes=axes, label='mixmod')
#qp.qp_plot_native(ens_f[15], axes=axes, label='flex')
leg = fig.legend()
We can also interpolate the function onto an evenly spaced grid point and cache those values with the gridded function.
grid = np.linspace(-3., 3., 100)
gridded = ens_n.pdf(grid)
cached_gridded = ens_n.gridded(grid)[1]
check = gridded - cached_gridded
print(check.min(), check.max())
0.0 0.0
Quantitative Comparisons
symm_lims = np.array([-1., 1.])
all_lims = [symm_lims, 2.*symm_lims, 3.*symm_lims]
Next, let’s compare the different parametrizations to the truth using the Kullback-Leibler Divergence (KLD). The KLD is a measure of how close two probability distributions are to one another – a smaller value indicates closer agreement. It is measured in units of bits of information, the information lost in going from the second distribution to the first distribution. The KLD calculator here takes in a shared grid upon which to evaluate the true distribution and the interpolated approximation of that distribution and returns the KLD of the approximation relative to the truth, which is not in general the same as the KLD of the truth relative to the approximation. Below, we’ll calculate the KLD of the approximation relative to the truth over different ranges, showing that it increases as it includes areas where the true distribution and interpolated distributions diverge.
# for a single pair of pdfs. (the 15th in each ensemble)
klds = qp.metrics.calculate_kld(ens_n, ens_s, limits=symm_lims)
print(klds)
[0.00366959 0.00371182 0.00315443 0.00044191 0.00343655 0.00453162
0.00166763 0.00350199 0.00060963 0.00325694 0.00379914 0.00413461
0.00146783 0.00375452 0.00340785 0.00038341 0.00089308 0.0024219
0.00440481 0.00269532 0.00149761 0.00435205 0.00282545 0.00335401
0.00251219 0.000994 0.00108705 0.00223901 0.00389532 0.00268311
0.00247407 0.00025512 0.00403371 0.0008904 0.00446524 0.00273664
0.00177683 0.00376072 0.00366289 0.00328361 0.00344289 0.00265748
0.00453964 0.0010088 0.00355582 0.00237827 0.00420921 0.00378332
0.00265234 0.00292326 0.00114919 0.00156324 0.00286138 0.00256105
0.0008692 0.00310536 0.00420783 0.00339083 0.00249314 0.00176598
0.00336663 0.00120526 0.00037089 0.00363469 0.00272068 0.00033138
0.00037614 0.00377794 0.00442546 0.00106841 0.00394342 0.00179919
0.00440221 0.0037757 0.00393809 0.00344012 0.00401084 0.00333697
0.00230136 0.00306801 0.00272743 0.00015566 0.00411218 0.00117511
0.0033651 0.00127422 0.0035337 0.00227786 0.00104108 0.00398201
0.00043188 0.00042997 0.00202351 0.00360442 0.00157278 0.0014487
0.00031734 0.00222945 0.00019697 0.00338896]
# Loop over all the other ensemble types
ensembles = [ens_n, ens_h, ens_i, ens_s, ens_q, ens_m]
for ensemble in ensembles[1:]:
D = []
for lims in all_lims:
klds = qp.metrics.calculate_kld(ens_n, ensemble, limits=lims)
D.append("%.2e +- %.2e" % (klds.mean(), klds.std()))
print(ensemble.gen_class.name + ' approximation: KLD over 1, 2, 3, sigma ranges = ' + str(D))
hist approximation: KLD over 1, 2, 3, sigma ranges = ['1.17e-02 +- 3.34e-03', '2.77e-02 +- 4.93e-03', '3.71e-02 +- 4.76e-03']
interp approximation: KLD over 1, 2, 3, sigma ranges = ['3.09e-02 +- 1.21e-02', '2.33e-02 +- 2.62e-03', '1.02e-02 +- 1.95e-03']
broken KLD: (array([-1.09200913e-03, -1.15274828e-03, -1.17316806e-04, 2.98469288e-03,
-6.30907754e-04, -3.27401498e-03, 2.71612685e-03, -6.97162313e-04,
5.42843806e-03, -2.37848264e-04, -1.38010311e-03, -2.12890982e-03,
3.42548914e-03, -1.19370195e-03, -5.13373147e-04, 3.36432998e-03,
4.43792025e-03, 2.30845034e-03, -2.63806241e-03, 6.47747562e-04,
2.66695759e-03, -2.61763706e-03, 4.61862291e-04, -4.40600617e-04,
1.02123880e-03, 3.48204398e-03, 2.11391009e-03, 2.13382392e-03,
-1.55021904e-03, 7.65561900e-04, 2.36325008e-03, 3.23039085e-03,
-1.93376926e-03, 5.16386655e-03, -3.04683135e-03, 1.08420747e-03,
3.14575135e-03, -1.29257708e-03, -1.07431344e-03, -3.25646342e-04,
-3.71052367e-04, 6.91371397e-04, -3.41174332e-03, 2.45948559e-03,
-8.52707067e-04, 1.00001821e-03, -2.42977350e-03, -6.00932678e-04,
9.13221711e-04, 2.42634833e-04, 3.47581453e-03, 2.44748472e-03,
9.50377518e-04, 1.59897229e-03, 4.94628865e-03, -8.78471671e-05,
-2.33422538e-03, -1.55724663e-04, 1.28270236e-03, 4.16764574e-03,
-5.02728766e-04, 2.27560787e-03, 3.82784630e-03, -6.80967129e-04,
2.00069620e-03, 4.19457645e-03, 2.56250321e-03, -1.32236357e-03,
-3.03565704e-03, 2.61841077e-03, -1.49759973e-03, 2.86792639e-03,
-2.89862484e-03, -1.31529205e-03, -1.66032306e-03, -6.24571040e-04,
-1.89155568e-03, -3.02536267e-04, 2.69388532e-03, -2.64840493e-05,
4.25718168e-04, 3.73140001e-03, -2.08044997e-03, 3.78867807e-03,
-2.55138101e-04, 3.36589252e-03, -1.66046473e-04, 3.11267336e-03,
2.96681526e-03, -1.70890071e-03, 2.97669618e-03, 2.35673789e-03,
3.33547667e-03, -8.46405026e-04, 1.64430508e-03, 4.81518695e-03,
4.26533461e-03, 1.45638790e-03, 4.17213692e-03, -5.24988662e-04]), array([[0.06374948, 0.06494121, 0.06614891, ..., 0.05519975, 0.05414232,
0.05310009],
[0.07680358, 0.07819532, 0.0796044 , ..., 0.04032144, 0.03947154,
0.03863573],
[0.03819755, 0.03897505, 0.03976483, ..., 0.10268134, 0.10112654,
0.0995864 ],
...,
[0.01611605, 0.01651205, 0.0169162 , ..., 0.16182971, 0.15977997,
0.15774147],
[0.00434846, 0.00448184, 0.00461884, ..., 0.2460105 , 0.24356705,
0.24112345],
[0.08503405, 0.0864493 , 0.08788006, ..., 0.0451456 , 0.04426107,
0.04338989]], shape=(100, 400)), array([[0.06371 , 0.06499321, 0.0662958 , ..., 0.05534102, 0.05419514,
0.05306721],
[0.07675814, 0.07825416, 0.07977194, ..., 0.04045474, 0.03952586,
0.03861286],
[0.03815968, 0.03899093, 0.03983578, ..., 0.1027111 , 0.10109 ,
0.09948766],
...,
[0.01609088, 0.01649852, 0.01691371, ..., 0.16159592, 0.15953752,
0.15749503],
[0.00433764, 0.00444268, 0.00454967, ..., 0.24513562, 0.24282825,
0.24052353],
[0.0849625 , 0.08645462, 0.08796623, ..., 0.04524142, 0.04428982,
0.04335338]], shape=(100, 400)), np.float64(0.010025062656641603))
spline approximation: KLD over 1, 2, 3, sigma ranges = ['2.58e-03 +- 1.30e-03', '2.22e-16 +- 0.00e+00', '1.08e-03 +- 2.90e-04']
broken KLD: (array([ 1.00049724e-02, 8.85752542e-03, 8.38931798e-03, -8.53260793e-03,
1.02965083e-02, 6.94717859e-03, -4.67332162e-03, 8.18564662e-03,
-1.90419414e-02, 7.53436472e-03, 9.86289169e-03, 7.59348682e-03,
-8.42713836e-03, 7.63290701e-03, 7.95773348e-03, -1.05560432e-02,
-1.39414721e-02, -5.06776033e-03, 5.38609022e-03, 4.64508107e-03,
-4.24981236e-03, 6.29298839e-03, 5.38943822e-03, 8.51278243e-03,
2.45881612e-03, -8.85354250e-03, -2.26594562e-03, -2.83423889e-03,
8.21855430e-03, 3.50597756e-03, -5.84745151e-03, -1.06979046e-02,
8.92474338e-03, -1.75901749e-02, 6.92323545e-03, 7.98254533e-04,
-7.51315744e-03, 9.91802243e-03, 9.84777452e-03, 8.60570082e-03,
4.85031887e-03, 4.52259369e-03, 7.78266113e-03, -3.83035560e-03,
9.86249647e-03, 3.34490741e-03, 9.03658145e-03, 1.98323787e-03,
2.45181109e-03, 7.17423827e-03, -8.60798624e-03, -3.18803608e-03,
8.44195167e-04, -1.10020924e-03, -1.65211158e-02, 9.89314984e-03,
7.61781726e-03, 3.60043654e-03, 7.96638094e-04, -1.39189314e-02,
1.01088592e-02, -2.66036058e-03, -1.26776699e-02, 4.36974694e-03,
-5.10185883e-03, -1.44631252e-02, -7.02387404e-03, 9.53407449e-03,
8.02017301e-03, -4.39946735e-03, 5.82736834e-03, -5.84863371e-03,
7.40164738e-03, 9.46565796e-03, 8.32378146e-03, 9.57842547e-03,
9.37653455e-03, 6.35007704e-03, -7.05488985e-03, 9.41928577e-03,
7.67932840e-03, -1.34282037e-02, 7.77251564e-03, -1.03003944e-02,
5.06400565e-03, -7.94954285e-03, 1.96502702e-03, -9.88787604e-03,
-6.13550586e-03, 7.27552215e-03, -8.55125718e-03, -5.80947032e-03,
-9.87775247e-03, 7.11249886e-03, -7.63544291e-05, -1.65263120e-02,
-1.48295059e-02, 6.26127710e-04, -1.50251519e-02, 9.27540784e-03]), array([[0.25283534, 0.25512887, 0.25741852, ..., 0.23536772, 0.23306003,
0.23075285],
[0.28302329, 0.28532901, 0.28762492, ..., 0.20541361, 0.20307159,
0.20073632],
[0.18373431, 0.18581891, 0.18791031, ..., 0.30050403, 0.29859386,
0.29666924],
...,
[0.11465316, 0.11638709, 0.11813616, ..., 0.36121834, 0.35997311,
0.35869855],
[0.05377097, 0.05486613, 0.0559779 , ..., 0.40035185, 0.40039892,
0.40040525],
[0.28126508, 0.28335575, 0.28543573, ..., 0.20527824, 0.20309454,
0.2009156 ]], shape=(100, 200)), array([[0.25685301, 0.25882082, 0.26078864, ..., 0.24193245, 0.23996463,
0.23799681],
[0.28374402, 0.28578452, 0.28782501, ..., 0.21694512, 0.21490463,
0.21286413],
[0.19608222, 0.19792241, 0.19976261, ..., 0.29801471, 0.2966305 ,
0.29493759],
...,
[0.1291595 , 0.13108354, 0.13300759, ..., 0.35210965, 0.35066236,
0.34921508],
[0.11834003, 0.11834003, 0.11834003, ..., 0.3888567 , 0.3888567 ,
0.3888567 ],
[0.28096232, 0.28284854, 0.28473477, ..., 0.2153248 , 0.21343858,
0.21155235]], shape=(100, 200)), np.float64(0.010050251256281407))
broken KLD: (array([-0.02568675, -0.02940021, -0.02428436, 0.184063 , -0.02132552,
-0.02180899, 0.16816214, -0.02840228, 0.23189591, -0.02781392,
-0.02821229, -0.03912496, 0.19924294, -0.03304185, -0.02799765,
0.19842221, 0.21842512, 0.19338238, 0.05977594, -0.02978387,
0.16530405, 0.01223661, -0.02931695, -0.02579599, 0.02459352,
0.19736262, 0.1104659 , 0.16182258, -0.03377389, -0.00208328,
0.19369538, 0.19336772, -0.03460878, 0.23156688, -0.02052733,
0.1057792 , 0.2009164 , -0.02741525, -0.02601504, -0.02478098,
0.01708606, -0.02968984, -0.04765311, 0.140494 , -0.0243492 ,
-0.02966008, -0.03787057, 0.13849071, 0.03616483, -0.02576151,
0.1967015 , 0.14632276, 0.11645049, 0.1452956 , 0.22501094,
-0.01866641, -0.04043473, 0.06814691, 0.08604625, 0.22218306,
-0.0209717 , 0.12854375, 0.2101029 , 0.05604655, 0.19872623,
0.21307988, 0.16692342, -0.02870002, -0.04443366, 0.15042952,
0.00822971, 0.18202141, -0.0447638 , -0.02883851, -0.0342364 ,
-0.02359573, -0.03319279, -0.03178193, 0.20287446, -0.02001076,
-0.02242605, 0.20577207, -0.03838398, 0.20672298, 0.0040183 ,
0.19677392, 0.12581776, 0.21299806, 0.17310835, -0.03717583,
0.1841523 , 0.14819814, 0.21093313, -0.03242522, 0.04237985,
0.23591497, 0.21541091, 0.0705692 , 0.21673729, -0.0238712 ]), array([[0.06374948, 0.06494121, 0.06614891, ..., 0.05519975, 0.05414232,
0.05310009],
[0.07680358, 0.07819532, 0.0796044 , ..., 0.04032144, 0.03947154,
0.03863573],
[0.03819755, 0.03897505, 0.03976483, ..., 0.10268134, 0.10112654,
0.0995864 ],
...,
[0.01611605, 0.01651205, 0.0169162 , ..., 0.16182971, 0.15977997,
0.15774147],
[0.00434846, 0.00448184, 0.00461884, ..., 0.2460105 , 0.24356705,
0.24112345],
[0.08503405, 0.0864493 , 0.08788006, ..., 0.0451456 , 0.04426107,
0.04338989]], shape=(100, 400)), array([[0.11485988, 0.11485988, 0.11485988, ..., 0.11485988, 0.11485988,
0.11485988],
[0.11696171, 0.11696171, 0.11696171, ..., 0.11696171, 0.11696171,
0.11696171],
[0.11107291, 0.11107291, 0.11107291, ..., 0.11550925, 0.11367367,
0.11183809],
...,
[0. , 0. , 0. , ..., 0.17688247, 0.17496325,
0.17304403],
[0. , 0. , 0. , ..., 0.25225699, 0.25017336,
0.24808973],
[0.11245347, 0.11245347, 0.11245347, ..., 0.11245347, 0.11245347,
0.11245347]], shape=(100, 400)), np.float64(0.010025062656641603))
quant approximation: KLD over 1, 2, 3, sigma ranges = ['2.22e-16 +- 0.00e+00', '2.22e-16 +- 0.00e+00', '3.13e-01 +- 8.85e-02']
broken KLD: (array([ 9.09849966e-03, -3.56037110e-06, 3.34297142e-02, 1.58370273e-02,
-2.31006056e-03, 1.12735012e-02, 4.71153635e-03, -1.13801385e-02,
3.64006509e-04, 8.08939324e-03, 2.16150110e-02, 8.93082896e-04,
-2.07472757e-02, -2.11998452e-03, 1.99662614e-02, -6.70140751e-03,
-2.92576740e-02, 1.55077411e-02, 2.94254616e-02, 6.35229411e-04,
-5.17016329e-03, 2.39464090e-03, 2.47209732e-02, -8.66179242e-03,
-8.23011110e-03, 3.09904987e-02, 1.13391828e-02, -1.68925085e-02,
2.19714681e-02, -8.83654746e-03, 4.42393960e-03, -5.10274953e-04,
3.03288475e-02, -1.64746499e-02, 1.01429953e-02, 6.26747672e-03,
1.40699311e-02, 6.39444216e-03, -4.14515039e-04, 1.36627232e-02,
6.68033360e-03, 2.28333148e-03, 1.23851680e-02, 7.81064294e-03,
1.54309708e-02, 4.41882174e-03, -1.67703239e-03, 3.40741368e-03,
-2.93980737e-02, 3.36015800e-02, -9.55547769e-03, -1.75522662e-02,
1.74917990e-02, 1.41062219e-02, 3.34077937e-03, 1.40847919e-02,
2.77540656e-02, 6.27832475e-03, -5.28241771e-03, -5.10657448e-03,
2.80431030e-02, -5.99397741e-03, 2.43165971e-02, 5.90406977e-03,
-1.12202588e-02, 1.68984256e-03, 3.15484460e-03, -7.36633986e-04,
4.11597520e-03, 3.25110866e-02, 3.78361550e-03, 1.05639022e-02,
3.04597440e-03, 1.70970751e-02, 2.25679894e-02, 2.35614471e-02,
1.94361650e-02, 7.07682204e-03, 1.28326517e-02, 4.59997344e-03,
-3.30786244e-03, -7.44517642e-04, 2.48332482e-02, -5.37214535e-03,
1.78900042e-02, 6.82406865e-03, 1.28147734e-03, -5.36975375e-04,
-1.76745119e-02, 1.50924241e-02, -3.75368260e-03, -1.87781007e-02,
-7.32778186e-03, 2.28451285e-02, -5.85915714e-03, 9.85372898e-03,
1.12536887e-02, 5.02626744e-03, 4.39949793e-03, 3.52963089e-02]), array([[0.25283534, 0.25512887, 0.25741852, ..., 0.23536772, 0.23306003,
0.23075285],
[0.28302329, 0.28532901, 0.28762492, ..., 0.20541361, 0.20307159,
0.20073632],
[0.18373431, 0.18581891, 0.18791031, ..., 0.30050403, 0.29859386,
0.29666924],
...,
[0.11465316, 0.11638709, 0.11813616, ..., 0.36121834, 0.35997311,
0.35869855],
[0.05377097, 0.05486613, 0.0559779 , ..., 0.40035185, 0.40039892,
0.40040525],
[0.28126508, 0.28335575, 0.28543573, ..., 0.20527824, 0.20309454,
0.2009156 ]], shape=(100, 200)), array([[0.22127026, 0.22246365, 0.22367884, ..., 0.23263353, 0.23079908,
0.2289759 ],
[0.24784355, 0.25062226, 0.25342573, ..., 0.22705923, 0.22495713,
0.22283999],
[0.21068171, 0.21180669, 0.21290816, ..., 0.28791655, 0.28539498,
0.28287145],
...,
[0.11968078, 0.12125437, 0.12283204, ..., 0.36780158, 0.36624963,
0.36467242],
[0.06332938, 0.06483999, 0.06637119, ..., 0.39695503, 0.39738716,
0.39774792],
[0.25988681, 0.26166959, 0.26347176, ..., 0.18690649, 0.18581692,
0.18473782]], shape=(100, 200)), np.float64(0.010050251256281407))
broken KLD: (array([-0.00580723, 0.0078711 , 0.01920279, 0.01708116, -0.00995442,
0.00295074, 0.00497018, 0.01101056, 0.00157328, 0.00282242,
-0.00032306, 0.01198066, 0.00168589, 0.00078715, 0.00884906,
-0.00984602, 0.00255404, 0.02464648, 0.00968176, 0.00085468,
0.00817605, 0.00585882, 0.01628805, 0.00059489, 0.00300147,
0.02125731, 0.00335574, 0.01885313, 0.01979052, 0.01621723,
0.01161756, 0.01808923, 0.00416054, 0.0133484 , 0.00576453,
0.00276682, 0.00383595, 0.00512808, -0.00456278, 0.01083553,
-0.00305952, 0.01371852, 0.00167282, 0.01448163, 0.02682956,
-0.00470869, -0.00801741, 0.00162821, 0.00601322, 0.01307176,
0.02038086, 0.01898675, 0.01834009, 0.00693048, 0.01611156,
0.02309598, 0.00506604, -0.00066992, 0.01104295, 0.00799444,
0.00763062, 0.00510977, 0.02308457, 0.00890217, -0.01098457,
0.01701861, 0.00534485, -0.00473083, -0.00291977, 0.02212077,
0.00673783, 0.01714977, -0.00139679, 0.01399519, 0.0083615 ,
-0.00517029, 0.00681127, 0.01445057, 0.00699739, 0.00509852,
0.00344916, 0.01403833, -0.00454519, 0.0139689 , 0.00646104,
0.00964889, 0.00421189, 0.00157965, 0.01351069, -0.00266975,
0.01030515, 0.00123504, 0.00456983, 0.01341526, 0.00279861,
0.01132574, 0.01072234, -0.00352136, 0.0220254 , 0.0146481 ]), array([[0.06374948, 0.06494121, 0.06614891, ..., 0.05519975, 0.05414232,
0.05310009],
[0.07680358, 0.07819532, 0.0796044 , ..., 0.04032144, 0.03947154,
0.03863573],
[0.03819755, 0.03897505, 0.03976483, ..., 0.10268134, 0.10112654,
0.0995864 ],
...,
[0.01611605, 0.01651205, 0.0169162 , ..., 0.16182971, 0.15977997,
0.15774147],
[0.00434846, 0.00448184, 0.00461884, ..., 0.2460105 , 0.24356705,
0.24112345],
[0.08503405, 0.0864493 , 0.08788006, ..., 0.0451456 , 0.04426107,
0.04338989]], shape=(100, 400)), array([[0.07693287, 0.07857245, 0.08022603, ..., 0.05356344, 0.05218226,
0.0508229 ],
[0.09215045, 0.0932665 , 0.09438189, ..., 0.03332354, 0.03231764,
0.03133468],
[0.04508609, 0.04631214, 0.04756006, ..., 0.11431891, 0.11317188,
0.11202622],
...,
[0.01356333, 0.01399718, 0.01444196, ..., 0.17524689, 0.1731799 ,
0.17111301],
[0.00186477, 0.00195448, 0.00204802, ..., 0.22704465, 0.2253091 ,
0.22359265],
[0.10443489, 0.10609392, 0.1077582 , ..., 0.06569636, 0.06446379,
0.06324157]], shape=(100, 400)), np.float64(0.010025062656641603))
broken KLD: (array([ 4.43061454e-03, 6.96350971e-03, 9.92006854e-03, 3.10095545e-03,
1.21033680e-03, 4.32684830e-03, -2.89015410e-03, 8.25521507e-03,
1.20974674e-03, 2.33587844e-05, 3.61276471e-03, 8.37024189e-03,
1.02709176e-02, 6.66465646e-04, 2.99514801e-03, -8.03957422e-03,
9.15797712e-03, 8.18492224e-03, 5.36886118e-03, -5.15175154e-04,
-1.36260478e-03, 3.30998025e-03, 9.58761924e-03, 6.01822493e-03,
1.46513546e-03, 1.76658207e-03, 9.57459665e-03, 8.52287767e-03,
8.13658707e-03, 7.36271074e-03, 3.00520038e-03, 9.30940966e-03,
5.37604799e-03, 1.00211878e-02, 6.26280356e-03, -5.37568765e-04,
-1.61219978e-04, 2.41727693e-03, 2.93057808e-03, 1.81228846e-03,
8.62848126e-04, 3.50436949e-03, 1.18293921e-03, 1.14896735e-02,
1.12397304e-02, -3.94707607e-03, 6.42719668e-03, 3.57894940e-03,
8.05444816e-03, 4.79554832e-03, 8.70551637e-03, 8.90032034e-03,
5.59686931e-03, 5.60821718e-03, 3.90130446e-03, 7.28182771e-03,
4.37061265e-03, 2.39817889e-03, 3.71279587e-03, 4.13201513e-03,
4.96149622e-03, 5.39532428e-03, 1.26187250e-02, 4.45130270e-03,
1.50686712e-03, 1.19677058e-02, 1.43516232e-02, 2.74605468e-03,
3.62116779e-03, 1.50296361e-02, 2.80298015e-03, 7.72883634e-03,
5.89566886e-03, 4.40107981e-03, 7.52807917e-03, 3.69483830e-03,
3.55066485e-03, 7.12430020e-03, 1.28216758e-03, 2.36052029e-03,
1.78507982e-03, 1.27270713e-02, 4.97560015e-03, 5.31324102e-03,
3.33575946e-03, 1.41558686e-02, 3.91205844e-03, 1.37361016e-03,
8.05073409e-03, 1.02505199e-03, 1.04929367e-03, -3.58763301e-03,
6.41418074e-03, 3.30654060e-03, -3.44965275e-04, 5.69020507e-03,
5.80329539e-03, -3.42236524e-03, 1.20753908e-02, 4.78772803e-03]), array([[0.00622268, 0.00639946, 0.00658062, ..., 0.00501176, 0.00486932,
0.00473049],
[0.00779079, 0.00801042, 0.00823542, ..., 0.0029586 , 0.00286788,
0.00277967],
[0.00326943, 0.00336571, 0.00346451, ..., 0.01444508, 0.01410063,
0.01376317],
...,
[0.0008957 , 0.00092626, 0.00095777, ..., 0.02866605, 0.02804142,
0.02742785],
[0.00012842, 0.0001337 , 0.00013918, ..., 0.05520262, 0.05410607,
0.05302593],
[0.01035191, 0.01062039, 0.01089483, ..., 0.00399798, 0.00388416,
0.00377324]], shape=(100, 600)), array([[2.46228437e-03, 2.58259895e-03, 2.70806597e-03, ...,
1.00823661e-03, 9.55347718e-04, 9.04979559e-04],
[1.06403191e-02, 1.09847588e-02, 1.13379376e-02, ...,
4.84974038e-04, 4.59308710e-04, 4.34898006e-04],
[9.17429096e-04, 9.65429232e-04, 1.01569424e-03, ...,
2.08417832e-02, 2.03171659e-02, 1.98022231e-02],
...,
[2.08240300e-04, 2.19361503e-04, 2.31028936e-04, ...,
2.23586517e-02, 2.16843370e-02, 2.10259769e-02],
[5.24516562e-06, 5.62821438e-06, 6.03781300e-06, ...,
7.58496650e-02, 7.44552327e-02, 7.30718833e-02],
[7.56509765e-03, 7.84945513e-03, 8.14273939e-03, ...,
3.60373202e-03, 3.46412686e-03, 3.32923974e-03]], shape=(100, 600)), np.float64(0.01001669449081803))
mixmod approximation: KLD over 1, 2, 3, sigma ranges = ['2.22e-16 +- 0.00e+00', '2.22e-16 +- 0.00e+00', '2.22e-16 +- 0.00e+00']
The progression of KLD values should follow that of the root mean square error (RMSE), another measure of how close two functions are to one another. The RMSE also increases as it includes areas where the true distribution and interpolated distribution diverge. Unlike the KLD, the RMSE is symmetric, meaning the distance measured is not that of one distribution from the other but of the symmetric distance between them.
for ensemble in ensembles[1:]:
D = []
for lims in all_lims:
rmses = qp.metrics.calculate_rmse(ens_n, ensemble, limits=lims)
D.append("%.2e +- %.2e" % (rmses.mean(), rmses.std()))
print(ensemble.gen_class.name + ' approximation: RMSE over 1, 2, 3, sigma ranges = ' + str(D))
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[35], line 4
2 D = []
3 for lims in all_lims:
----> 4 rmses = qp.metrics.calculate_rmse(ens_n, ensemble, limits=lims)
5 D.append("%.2e +- %.2e" % (rmses.mean(), rmses.std()))
6 print(ensemble.gen_class.name + ' approximation: RMSE over 1, 2, 3, sigma ranges = ' + str(D))
File ~/checkouts/readthedocs.org/user_builds/qp/envs/latest/lib/python3.10/site-packages/qp/metrics/metrics.py:160, in calculate_rmse(p, q, limits, dx)
134 """
135 Calculates the Root Mean Square Error between two qp.Ensemble objects.
136
(...)
157 TO DO: change dx to N
158 """
159 if p.shape != q.shape:
--> 160 raise ValueError(
161 "Cannot calculate RMSE between two ensembles with different shapes"
162 )
164 # Make a grid from the limits and resolution
165 grid = _calculate_grid_parameters(limits, dx)
ValueError: Cannot calculate RMSE between two ensembles with different shapes
Both the KLD and RMSE metrics suggest that the quantile approximation is better in the high density region, but samples work better when the tails are included. We might expect the answer to the question of which approximation to use to depend on the application, and whether the tails need to be captured or not.
Storing and retreiving ensembles
You can store and retrieve ensembles from disk using the qp.Ensemble.write_to and qp.read methods.
These work in two steps, first they convert the Ensemble data to astropy.table objects, and then they write the tables. This means you can store the data in any format support by astropy.
tabs = ens_n.build_tables()
print(tabs.keys())
print()
print("Meta Data")
print(tabs['meta'])
print()
print("Object Data")
print(tabs['data'])
dict_keys(['meta', 'data'])
Meta Data
{'pdf_name': array([b'norm'], dtype='|S4'), 'pdf_version': array([0])}
Object Data
{'loc': array([[-0.04815224],
[-0.17455599],
[ 0.26994988],
[-0.9438302 ],
[-0.05271168],
[-0.15623651],
[ 0.76972841],
[ 0.24415093],
[ 0.95407152],
[-0.30514607],
[-0.02730188],
[ 0.19467777],
[-0.82224738],
[ 0.24489467],
[-0.26728835],
[ 0.96010884],
[-0.91048158],
[ 0.70021527],
[ 0.26876681],
[ 0.46147515],
[ 0.78417283],
[-0.23203026],
[ 0.42699502],
[ 0.24138619],
[ 0.54127299],
[ 0.87715334],
[ 0.80604695],
[ 0.69035976],
[ 0.19231423],
[-0.49782248],
[-0.70359358],
[-0.97797273],
[ 0.11163251],
[ 0.92192972],
[ 0.17434699],
[ 0.57254554],
[ 0.78336644],
[ 0.03074603],
[ 0.08198767],
[ 0.24332715],
[-0.3964033 ],
[-0.46817377],
[-0.07263986],
[ 0.83744027],
[-0.10296287],
[-0.5252825 ],
[-0.02098667],
[-0.4672512 ],
[ 0.53110688],
[ 0.35139861],
[ 0.85812597],
[ 0.76388526],
[-0.56271395],
[-0.63418735],
[-0.92109247],
[-0.16395277],
[-0.18085865],
[ 0.44234279],
[-0.58981155],
[ 0.82576381],
[ 0.10336635],
[ 0.79837245],
[-0.96850503],
[ 0.39644255],
[ 0.67543962],
[ 0.9779488 ],
[-0.944671 ],
[ 0.10000184],
[ 0.0948455 ],
[ 0.83619061],
[ 0.30743757],
[ 0.76713118],
[ 0.15608875],
[ 0.10878276],
[-0.17899016],
[-0.1526427 ],
[-0.05521282],
[ 0.34979033],
[ 0.7299227 ],
[-0.20975738],
[ 0.34423819],
[-0.9991165 ],
[ 0.18728735],
[-0.8645259 ],
[-0.39549337],
[ 0.8401505 ],
[ 0.48451351],
[ 0.75513166],
[ 0.85435285],
[-0.23480032],
[ 0.9452581 ],
[-0.92929214],
[-0.77718546],
[ 0.28828602],
[-0.69227265],
[ 0.86623961],
[-0.9805002 ],
[ 0.61461543],
[ 0.99653732],
[-0.18491633]]), 'scale': array([[1.02652848],
[1.00808147],
[1.06152735],
[1.04771426],
[1.05664805],
[0.90089844],
[0.99609112],
[1.02188863],
[0.92650824],
[1.03834458],
[1.01034733],
[0.9495649 ],
[0.96864358],
[0.98890826],
[1.02866674],
[1.02912959],
[0.95453265],
[0.94030169],
[0.9016703 ],
[1.05678472],
[1.01381402],
[0.91479379],
[1.05372917],
[1.04191843],
[1.04175948],
[0.99731179],
[1.08103332],
[0.97612083],
[0.98120037],
[1.04060429],
[0.93035021],
[1.03933086],
[0.97414753],
[0.9221706 ],
[0.90824023],
[0.99002036],
[0.95844986],
[1.01527575],
[1.02538558],
[1.05075973],
[0.98326831],
[1.05906413],
[0.90627333],
[1.05926846],
[1.03755034],
[1.07110147],
[0.95589884],
[0.91577658],
[1.02695227],
[1.06886591],
[0.98842929],
[1.02408925],
[0.97885433],
[0.97301947],
[0.93280155],
[1.08917596],
[0.94214005],
[0.9713751 ],
[1.01425569],
[0.905162 ],
[1.06217644],
[1.06216786],
[1.00669878],
[0.9592217 ],
[0.92133703],
[0.99110986],
[1.07375423],
[1.00887777],
[0.92192841],
[1.04538593],
[0.94932131],
[0.97370558],
[0.91923894],
[1.00832953],
[0.97787262],
[1.04681549],
[0.98143504],
[1.0137707 ],
[0.93074056],
[1.08637729],
[1.0983295 ],
[1.01730973],
[0.95373408],
[0.96971057],
[0.99363926],
[0.98641192],
[0.93597119],
[0.91013057],
[1.02455649],
[0.96179219],
[1.04843366],
[1.08551655],
[0.92364239],
[0.99754726],
[1.09266752],
[0.90063662],
[0.98843194],
[1.03813924],
[0.99634025],
[1.0484953 ]])}
Here is a loopback test showing that we get the same results before and after a write/read cycle.
suffix_list = ['_n', '_h', '_i', '_s', '_q', '_m']
filetypes = ['fits', 'hf5']
for ens, suffix in zip(ensembles, suffix_list):
for ft in filetypes:
outfile = "test%s.%s" % (suffix, ft)
metafile = "test%s_meta.%s" % (suffix, ft)
pdf_1 = ens.pdf(bins)
ens.write_to(outfile)
ens_r = qp.read(outfile)
pdf_2 = ens_r.pdf(bins)
check = pdf_1 - pdf_2
print(suffix, ft, check.min(), check.max())
os.unlink(outfile)
try:
os.unlink(metafile)
except Exception:
pass
_n fits 0.0 0.0
_n hf5 0.0 0.0
_h fits -5.551115123125783e-17 1.1102230246251565e-16
_h hf5 -5.551115123125783e-17 1.1102230246251565e-16
_i fits -1.1102230246251565e-16 5.551115123125783e-17
_i hf5 -1.1102230246251565e-16 5.551115123125783e-17
_s fits 0.0 0.0
_s hf5 0.0 0.0
_q fits 0.0 0.0
_q hf5 0.0 0.0
_m fits -1.1102230246251565e-16 0.0
_m hf5 -1.1102230246251565e-16 0.0
Finally, we can compare the moments of each approximation and compare those to the moments of the true distribution.
which_moments = range(3)
all_moments = []
for ens in ensembles:
moments = []
for n in which_moments:
moms = qp.metrics.calculate_moment(ens, n, limits=(-3, 3))
moments.append("%.2e +- %.2e" % (moms.mean(), moms.std()))
all_moments.append(moments)
print('moments: '+str(which_moments))
for ens, mom in zip(ensembles, all_moments):
print(ens.gen_class.name+': '+str(mom))
moments: range(0, 3)
norm: ['9.91e-01 +- 7.13e-03', '1.05e-01 +- 5.51e-01', '1.25e+00 +- 2.71e-01']
hist: ['9.92e-01 +- 7.33e-03', '1.04e-01 +- 5.50e-01', '1.41e+00 +- 2.61e-01']
interp: ['9.87e-01 +- 9.39e-03', '1.02e-01 +- 5.35e-01', '1.36e+00 +- 2.40e-01']
spline: ['9.90e-01 +- 7.28e-03', '1.05e-01 +- 5.52e-01', '1.25e+00 +- 2.75e-01']
quant: ['1.07e+00 +- 1.60e-02', '1.15e-01 +- 5.87e-01', '1.57e+00 +- 2.54e-01']
mixmod: ['9.92e-01 +- 8.51e-03', '1.08e-01 +- 5.56e-01', '1.28e+00 +- 2.73e-01']