from statsmodels.compat.pandas import deprecate_kwarg
import numpy as np
from statsmodels.tools.rng_qrng import check_random_state
def _make_index(prob, size, rng=None):
"""
Returns a boolean index for given probabilities.
Parameters
----------
prob : array_like
Probability of sampling from each distribution in dist.
size : int
The length of the returned sample.
rng : int, array_like of int, numpy.random.Generator, or numpy.random.RandomState, optional
If `rng` is None, a new ``Generator`` is created using fresh
entropy from the operating system. If `rng` is an int or array
of ints, a new ``Generator`` is created, seeded with `rng`. If
`rng` is already a ``Generator`` or ``RandomState`` instance,
that instance is used.
Notes
-----
prob = [.75,.25] means that there is a 75% chance of the first column
being True and a 25% chance of the second column being True. The
columns are mutually exclusive.
"""
rng = check_random_state(rng)
rv = rng.uniform(size=(size, 1))
cumprob = np.cumsum(prob)
return np.logical_and(np.r_[0, cumprob[:-1]] <= rv, rv < cumprob)
[docs]
def mixture_rvs(prob, size, dist, kwargs=None, rng=None):
"""
Sample from a mixture of distributions.
Parameters
----------
prob : array_like
Probability of sampling from each distribution in dist
size : int
The length of the returned sample.
dist : array_like
An iterable of distributions objects from scipy.stats.
kwargs : tuple of dicts, optional
A tuple of dicts. Each dict in kwargs can have keys loc, scale, and
args to be passed to the respective distribution in dist. If not
provided, the distribution defaults are used.
rng : int, array_like of int, numpy.random.Generator, or numpy.random.RandomState, optional
If `rng` is None, the legacy global (singleton) ``RandomState``
provided by ``numpy.random`` is used; this behavior is
deprecated and will change to creating a new ``Generator``
using fresh entropy from the operating system in a future
release. If `rng` is an int or array of ints, a new
``RandomState`` instance is created, seeded with `rng`. If `rng`
is already a ``Generator`` or ``RandomState`` instance, that
instance is used.
Returns
-------
ndarray
Sample from the mixture distribution, with length `size`.
Examples
--------
Say we want 5000 random variables from mixture of normals with two
distributions norm(-1,.5) and norm(1,.5) and we want to sample from the
first with probability .75 and the second with probability .25.
>>> from scipy import stats
>>> prob = [.75,.25]
>>> Y = mixture_rvs(prob, 5000, dist=[stats.norm, stats.norm],
... kwargs = (dict(loc=-1,scale=.5),dict(loc=1,scale=.5)))
"""
if len(prob) != len(dist):
raise ValueError("You must provide as many probabilities as distributions")
if not np.allclose(np.sum(prob), 1):
raise ValueError("prob does not sum to 1")
if kwargs is None:
kwargs = ({},) * len(prob)
idx = _make_index(prob, size, rng)
sample = np.empty(size)
for i in range(len(prob)):
sample_idx = idx[..., i]
sample_size = sample_idx.sum()
loc = kwargs[i].get("loc", 0)
scale = kwargs[i].get("scale", 1)
args = kwargs[i].get("args", ())
sample[sample_idx] = dist[i].rvs(
*args, **dict(loc=loc, scale=scale, size=sample_size, random_state=rng)
)
return sample
[docs]
class MixtureDistribution:
"""univariate mixture distribution
for simple case for now (unbound support)
does not yet inherit from scipy.stats.distributions
adding pdf to mixture_rvs, some restrictions on broadcasting
Currently it does not hold any state, all arguments included in each method.
"""
# def __init__(self, prob, size, dist, kwargs=None):
[docs]
@deprecate_kwarg("random_state", "rng")
def rvs(self, prob, size, dist, kwargs=None, rng=None):
"""
Sample from a mixture of distributions.
Parameters
----------
prob : array_like
Probability of sampling from each distribution in dist
size : int
The length of the returned sample.
dist : array_like
An iterable of distributions objects from scipy.stats.
kwargs : tuple of dicts, optional
A tuple of dicts. Each dict in kwargs can have keys loc, scale, and
args to be passed to the respective distribution in dist. If not
provided, the distribution defaults are used.
rng : int, array_like of int, numpy.random.Generator, or numpy.random.RandomState, optional
If `rng` is None, the legacy global (singleton) ``RandomState``
provided by ``numpy.random`` is used; this behavior is
deprecated and will change to creating a new ``Generator``
using fresh entropy from the operating system in a future
release. If `rng` is an int or array of ints, a new
``RandomState`` instance is created, seeded with `rng`. If
`rng` is already a ``Generator`` or ``RandomState`` instance,
that instance is used.
rng : int, array_like of int, numpy.random.Generator, or numpy.random.RandomState, optional
.. deprecated:: 0.15
random_state has been deprecated. In-line with SPEC-007, use
rng for passing a random number generator or seed.
Returns
-------
ndarray
Sample from the mixture distribution.
"""
return mixture_rvs(prob, size, dist, kwargs=kwargs, rng=rng)
[docs]
def pdf(self, x, prob, dist, kwargs=None):
"""
pdf a mixture of distributions.
Parameters
----------
x : array_like
Array containing locations where the PDF should be evaluated
prob : array_like
Probability of sampling from each distribution in dist
dist : array_like
An iterable of distributions objects from scipy.stats.
kwargs : tuple of dicts, optional
A tuple of dicts. Each dict in kwargs can have keys loc, scale, and
args to be passed to the respective distribution in dist. If not
provided, the distribution defaults are used.
Returns
-------
ndarray
Pdf of the mixture distribution evaluated at `x`.
Examples
--------
Say we want 5000 random variables from mixture of normals with two
distributions norm(-1,.5) and norm(1,.5) and we want to sample from the
first with probability .75 and the second with probability .25.
>>> import numpy as np
>>> from scipy import stats
>>> from statsmodels.distributions.mixture_rvs import MixtureDistribution
>>> x = np.arange(-4.0, 4.0, 0.01)
>>> prob = [.75,.25]
>>> mixture = MixtureDistribution()
>>> Y = mixture.pdf(x, prob, dist=[stats.norm, stats.norm],
... kwargs = (dict(loc=-1,scale=.5),dict(loc=1,scale=.5)))
"""
if len(prob) != len(dist):
raise ValueError("You must provide as many probabilities as distributions")
if not np.allclose(np.sum(prob), 1):
raise ValueError("prob does not sum to 1")
if kwargs is None:
kwargs = ({},) * len(prob)
for i in range(len(prob)):
loc = kwargs[i].get("loc", 0)
scale = kwargs[i].get("scale", 1)
args = kwargs[i].get("args", ())
if i == 0: # assume all broadcast the same as the first dist
pdf_ = prob[i] * dist[i].pdf(x, *args, loc=loc, scale=scale)
else:
pdf_ += prob[i] * dist[i].pdf(x, *args, loc=loc, scale=scale)
return pdf_
[docs]
def cdf(self, x, prob, dist, kwargs=None):
"""
cdf of a mixture of distributions.
Parameters
----------
x : array_like
Array containing locations where the CDF should be evaluated
prob : array_like
Probability of sampling from each distribution in dist
dist : array_like
An iterable of distributions objects from scipy.stats.
kwargs : tuple of dicts, optional
A tuple of dicts. Each dict in kwargs can have keys loc, scale, and
args to be passed to the respective distribution in dist. If not
provided, the distribution defaults are used.
Returns
-------
ndarray
Cdf of the mixture distribution evaluated at `x`.
Examples
--------
Say we want 5000 random variables from mixture of normals with two
distributions norm(-1,.5) and norm(1,.5) and we want to sample from the
first with probability .75 and the second with probability .25.
>>> import numpy as np
>>> from scipy import stats
>>> from statsmodels.distributions.mixture_rvs import MixtureDistribution
>>> x = np.arange(-4.0, 4.0, 0.01)
>>> prob = [.75,.25]
>>> mixture = MixtureDistribution()
>>> Y = mixture.cdf(x, prob, dist=[stats.norm, stats.norm],
... kwargs = (dict(loc=-1,scale=.5),dict(loc=1,scale=.5)))
"""
if len(prob) != len(dist):
raise ValueError("You must provide as many probabilities as distributions")
if not np.allclose(np.sum(prob), 1):
raise ValueError("prob does not sum to 1")
if kwargs is None:
kwargs = ({},) * len(prob)
for i in range(len(prob)):
loc = kwargs[i].get("loc", 0)
scale = kwargs[i].get("scale", 1)
args = kwargs[i].get("args", ())
if i == 0: # assume all broadcast the same as the first dist
cdf_ = prob[i] * dist[i].cdf(x, *args, loc=loc, scale=scale)
else:
cdf_ += prob[i] * dist[i].cdf(x, *args, loc=loc, scale=scale)
return cdf_
[docs]
def mv_mixture_rvs(prob, size, dist, nvars, rng=None, **kwargs):
"""
Sample from a mixture of multivariate distributions.
Parameters
----------
prob : array_like
Probability of sampling from each distribution in dist
size : int
The length of the returned sample.
dist : array_like
An iterable of distributions instances with callable method rvs.
nvars : int
dimension of the multivariate distribution, could be inferred instead
rng : int, array_like of int, numpy.random.Generator, or numpy.random.RandomState, optional
If `rng` is None, the legacy global (singleton) ``RandomState``
provided by ``numpy.random`` is used; this behavior is
deprecated and will change to creating a new ``Generator``
using fresh entropy from the operating system in a future
release. If `rng` is an int or array of ints, a new
``RandomState`` instance is created, seeded with `rng`. If `rng`
is already a ``Generator`` or ``RandomState`` instance, that
instance is used.
**kwargs
Ignored.
Returns
-------
ndarray
Sample from the mixture of multivariate distributions, with shape
(`size`, `nvars`).
Examples
--------
Say we want 2000 random variables from mixture of normals with two
multivariate normal distributions, and we want to sample from the
first with probability .4 and the second with probability .6.
import statsmodels.sandbox.distributions.mv_normal as mvd
cov3 = np.array([[ 1. , 0.5 , 0.75],
[ 0.5 , 1.5 , 0.6 ],
[ 0.75, 0.6 , 2. ]])
mu = np.array([-1, 0.0, 2.0])
mu2 = np.array([4, 2.0, 2.0])
mvn3 = mvd.MVNormal(mu, cov3)
mvn32 = mvd.MVNormal(mu2, cov3/2., 4)
rvs = mix.mv_mixture_rvs([0.4, 0.6], 2000, [mvn3, mvn32], 3)
"""
if len(prob) != len(dist):
raise ValueError("You must provide as many probabilities as distributions")
if not np.allclose(np.sum(prob), 1):
raise ValueError("prob does not sum to 1")
# if kwargs is None:
# kwargs = ({},) * len(prob)
idx = _make_index(prob, size, rng)
sample = np.empty((size, nvars))
for i in range(len(prob)):
sample_idx = idx[..., i]
sample_size = sample_idx.sum()
# loc = kwargs[i].get('loc',0)
# scale = kwargs[i].get('scale',1)
# args = kwargs[i].get('args',())
# use int to avoid numpy bug with np.random.multivariate_normal
try:
sample[sample_idx] = dist[i].rvs(size=int(sample_size), rng=rng)
except TypeError:
sample[sample_idx] = dist[i].rvs(size=int(sample_size), random_state=rng)
return sample
