"""
Created on Wed Feb 17 15:35:23 2021
Author: Josef Perktold
License: BSD-3
"""
from statsmodels.compat.pandas import deprecate_kwarg
import numpy as np
from scipy import stats
from statsmodels.distributions.tools import (
_eval_bernstein_1d,
_eval_bernstein_2d,
_eval_bernstein_dd,
_Grid,
cdf2prob_grid,
prob2cdf_grid,
)
from statsmodels.tools._decorators import cache_readonly
from statsmodels.tools.rng_qrng import check_random_state
[docs]
class BernsteinDistribution:
"""Distribution based on Bernstein Polynomials on unit hypercube.
Parameters
----------
cdf_grid : array_like
cdf values on a equal spaced grid of the unit hypercube [0, 1]^d.
The dimension of the arrays define how many random variables are
included in the multivariate distribution.
Attributes
----------
cdf_grid : grid of cdf values
prob_grid : grid of cell or bin probabilities
k_dim : (int) number of components, dimension of random variable
k_grid : (tuple) shape of cdf_grid
k_grid_product : (int) total number of bins in grid
_grid : Grid instance with helper methods and attributes
"""
def __init__(self, cdf_grid):
self.cdf_grid = cdf_grid = np.asarray(cdf_grid)
self.k_dim = cdf_grid.ndim
self.k_grid = cdf_grid.shape
self.k_grid_product = np.prod([i - 1 for i in self.k_grid])
self._grid = _Grid(self.k_grid)
[docs]
@classmethod
def from_data(cls, data, k_bins):
"""Create distribution instance from data using histogram binning.
Classmethod to construct a distribution instance.
Parameters
----------
data : array_like
Data with observation in rows and random variables in columns.
Data can be 1-dimensional in the univariate case.
k_bins : int or list
Number or edges of bins to be used in numpy histogramdd.
If k_bins is a scalar int, then the number of bins of each
component will be equal to it.
Returns
-------
Instance of a Bernstein distribution
"""
data = np.asarray(data)
if np.any(data < 0) or np.any(data > 1):
raise ValueError("data needs to be in [0, 1]")
if data.ndim == 1:
data = data[:, None]
k_dim = data.shape[1]
if np.size(k_bins) == 1:
k_bins = [k_bins] * k_dim
bins = [np.linspace(-1 / ni, 1, ni + 2) for ni in k_bins]
c, e = np.histogramdd(data, bins=bins, density=False)
# TODO: check when we have zero observations, which bin?
# check bins start at 0 exept leading bin
assert all(ei[1] == 0 for ei in e)
c /= len(data)
cdf_grid = prob2cdf_grid(c)
return cls(cdf_grid)
[docs]
@cache_readonly
def prob_grid(self):
return cdf2prob_grid(self.cdf_grid, prepend=None)
[docs]
def cdf(self, x):
"""cdf values evaluated at x.
Parameters
----------
x : array_like
Points of multivariate random variable at which cdf is evaluated.
This can be a single point with length equal to the dimension of
the random variable, or two dimensional with points (observations)
in rows and random variables in columns.
In the univariate case, a 1-dimensional x will be interpreted as
different points for evaluation.
Returns
-------
pdf values
Notes
-----
Warning: 2-dim x with many points can be memory intensive because
currently the bernstein polynomials will be evaluated in a fully
vectorized computation.
"""
x = np.asarray(x)
if x.ndim == 1 and self.k_dim == 1:
x = x[:, None]
cdf_ = _eval_bernstein_dd(x, self.cdf_grid)
return cdf_
[docs]
def pdf(self, x):
"""pdf values evaluated at x.
Parameters
----------
x : array_like
Points of multivariate random variable at which pdf is evaluated.
This can be a single point with length equal to the dimension of
the random variable, or two dimensional with points (observations)
in rows and random variables in columns.
In the univariate case, a 1-dimensional x will be interpreted as
different points for evaluation.
Returns
-------
cdf values
Notes
-----
Warning: 2-dim x with many points can be memory intensive because
currently the bernstein polynomials will be evaluated in a fully
vectorized computation.
"""
x = np.asarray(x)
if x.ndim == 1 and self.k_dim == 1:
x = x[:, None]
# TODO: check usage of k_grid_product. Should this go into eval?
pdf_ = self.k_grid_product * _eval_bernstein_dd(x, self.prob_grid)
return pdf_
[docs]
def get_marginal(self, idx):
"""Get marginal BernsteinDistribution.
Parameters
----------
idx : int or list of int
Index or indices of the component for which the marginal
distribution is returned.
Returns
-------
BernsteinDistribution instance for the marginal distribution.
"""
# univariate
if self.k_dim == 1:
return self
sl = [-1] * self.k_dim
if np.shape(idx) == ():
idx = [idx]
for ii in idx:
sl[ii] = slice(None, None, None)
cdf_m = self.cdf_grid[tuple(sl)]
bpd_marginal = BernsteinDistribution(cdf_m)
return bpd_marginal
[docs]
@deprecate_kwarg("random_state", "rng")
def rvs(self, nobs, rng=None):
"""Generate random numbers from distribution.
Parameters
----------
nobs : int
Number of random observations to generate.
rng : {None, int, array_like[int], numpy.random.Generator, 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.
random_state : {None, int, array_like[int], numpy.random.Generator, 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.
"""
rng = check_random_state(rng)
rvs_mnl = rng.multinomial(nobs, self.prob_grid.flatten())
k_comp = self.k_dim
rvs_m = []
for i in range(len(rvs_mnl)):
if rvs_mnl[i] != 0:
idx = np.unravel_index(i, self.prob_grid.shape)
rvsi = []
for j in range(k_comp):
n = self.k_grid[j]
xgi = self._grid.x_marginal[j][idx[j]]
# Note: x_marginal starts at 0
# x_marginal ends with 1 but that is not used by idx
rvsi.append(
stats.beta.rvs(
n * xgi + 1,
n * (1 - xgi) + 0,
size=rvs_mnl[i],
random_state=rng,
)
)
rvs_m.append(np.column_stack(rvsi))
rvsm = np.concatenate(rvs_m)
return rvsm
[docs]
class BernsteinDistributionBV(BernsteinDistribution):
[docs]
def cdf(self, x):
cdf_ = _eval_bernstein_2d(x, self.cdf_grid)
return cdf_
[docs]
def pdf(self, x):
# TODO: check usage of k_grid_product. Should this go into eval?
pdf_ = self.k_grid_product * _eval_bernstein_2d(x, self.prob_grid)
return pdf_
[docs]
class BernsteinDistributionUV(BernsteinDistribution):
[docs]
def cdf(self, x, method="binom"):
cdf_ = _eval_bernstein_1d(x, self.cdf_grid, method=method)
return cdf_
[docs]
def pdf(self, x, method="binom"):
# TODO: check usage of k_grid_product. Should this go into eval?
pdf_ = self.k_grid_product * _eval_bernstein_1d(
x, self.prob_grid, method=method
)
return pdf_
