Select Git revision
AdditionalPDFs.py
-
Antonello, Dr. Massimiliano authored
This reverts commit 16bb59d0
Antonello, Dr. Massimiliano authoredThis reverts commit 16bb59d0
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
AdditionalPDFs.py 3.60 KiB
import numpy as np
from scipy.stats import rv_discrete, rv_continuous, uniform
import scipy.special as sc
import matplotlib.pyplot as plt
from scipy.stats._distn_infrastructure import (
rv_discrete, _ncx2_pdf, _ncx2_cdf, get_distribution_names)
class gpd_gen(rv_discrete):
def _argcheck(self, mu, lbda):
return mu >= 0.0 and lbda >= 0.0 and lbda <= 1.0
def _rvs(self, mu, lbda):
population = np.asarray(
self._random_state.poisson(mu, self._size)
)
if population.shape == ():
population = population.reshape(-1)
offspring = population.copy()
while np.any(offspring > 0):
# probability dists are NOT ufuncs
# print("offspring", offspring)
offspring[:] = [
self._random_state.poisson(m)
for m in lbda*offspring
]
population += offspring
return population
def _pmf(self, k, mu, lbda):
return np.exp(self._logpmf(k, mu, lbda))
def _logpmf(self, k, mu, lbda):
mu_pls_klmb = mu + lbda*k
return np.log(mu) + sc.xlogy(k-1, mu_pls_klmb) - mu_pls_klmb - sc.gammaln(k+1)
def _munp(self, n, mu, lbda):
if n == 1:
return mu/(1-lbda)
elif n == 2:
return (mu/(1-lbda))**2+mu/(1-lbda)**3
gpoisson = gpd_gen(name='gpoisson')
class borel_gen(rv_discrete):
def _argcheck(self, mu):
return ((mu > 0) & (mu<1))
def _logpmf(self, k, mu):
n = k+1
Pk = sc.xlogy(n-1, mu*n) - sc.gammaln(n + 1) - mu*n
return Pk
def _pmf(self, k, mu):
return np.exp(self._logpmf(k, mu))
# def _rvs(self, mu, size=None, random_state=None):
# u = uniform.rvs(loc=0, scale = 1, size=size)
# cum = np.cumsum([self._pmf(_k, mu) for _k in range(0, 100)])
# print(cum)
# rnd = [ np.argmax( cum>=_u ) for _u in u ]
# return rnd
def _rvs(self, mu, size=None, random_state=None, epsilon=1e-4):
_u = uniform.rvs(loc=0, scale = 1-epsilon, size=size)
_sum = 0
_k=0
_elem = []
_max_u = np.max(_u)
while(_sum<_max_u):
_pmf = self._pmf(_k, mu)
_elem.append(_pmf)
_sum+=_pmf
_k+=1
_cum = np.cumsum(_elem)
_rnd = [ np.argmax( _cum>=__u ) for __u in _u ]
return _rnd
def _stats(self, mu):
_mu = 1/(1-mu)
_var = mu/(1-mu)**3
tmp = np.asarray(mu)
mu_nonzero = ((tmp > 0) & (tmp<1))
#g1 and g2: Lagrangian Probability Distributions, 978-0-8176-4365-2, page 159
g1 = scipy._lib._util._lazywhere(mu_nonzero, (tmp,), lambda x: (1+2*x)/scipy.sqrt(x*(1-x)), np.inf)
g2 = scipy._lib._util._lazywhere(mu_nonzero, (tmp,), lambda x: 3 + (1 + 8*x+6*x**2)/(x*(1-x)), np.inf)
return _mu, _var, g1, g2
borel= borel_gen(name='borel')
class erlang_gen(rv_discrete):
def _pdf(self, x, a):
# gamma.pdf(x, a) = x**(a-1) * exp(-x) / gamma(a)
return np.exp(self._logpdf(x, a))
def _logpdf(self, k, mu, nu):
return sc.xlogy(a-1.0, x) - x - sc.gammaln(a)
# def _rvs(self, mu, nu, size=None, random_state=None):
# u = scipy.stats.uniform.rvs(loc=0, scale = 1, size=size)
# cum = np.cumsum([self._pmf(_k, mu, nu) for _k in range(0, 100)])
# rnd = [ np.argmax( cum>=_u ) for _u in u ]
# return rnd
pairs = list(globals().items())
_distn_names, _distn_gen_names = get_distribution_names(pairs, rv_discrete)
__all__ = _distn_names + _distn_gen_names