import warnings
from abc import ABCMeta, abstractmethod
from collections import namedtuple
from functools import partial
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as sp_optimize
from scipy import interpolate, linalg, special, stats
from scipy.integrate import cumulative_trapezoid, quad
try:
from multiprocessing import cpu_count
import numexpr as ne
ne.set_num_threads(min(64, cpu_count()))
NE = True
except ImportError:
NE = False
import varwg
from varwg import helpers as my
from varwg.time_series_analysis import _kde as kde
# def _build_owens():
# if sys.platform != "win32":
# # import urllib.request, urllib.error, urllib.parse
# from urllib import request, error
# import socket
# socket.setdefaulttimeout(10)
# warn = False
# url = "http://people.sc.fsu.edu/~jburkardt/f_src/owens/owens.f90"
# try:
# src_dir = os.path.dirname(__file__)
# except NameError:
# src_dir = os.path.abspath(".")
# # if os.path.exists(os.path.join(src_dir, "owens.f90")):
# # warn = True
# if not warn:
# with my.chdir(src_dir):
# with open("owens.f90", "w") as owens_file:
# try:
# content_full = (request
# .urlopen(url)
# .read()
# .decode("utf-8"))
# # under python 3: "character ( len = 9 )" causes
# # f2py to fail, so we cut out the subroutine
# # timestep, which is not needed for owens-t
# content = []
# keep_line = True
# for line in content_full.split("\n"):
# if line.startswith("subroutine timestamp"):
# keep_line = False
# if keep_line:
# content += [line]
# if not keep_line and line == "end":
# keep_line = True
# content = os.linesep.join(content)
# owens_file.write(content + os.linesep)
# except error.URLError:
# warn = True
# try:
# subprocess.call("f2py -L$PREFIX$/lib -c -m owens owens.f90",
# shell=True)
# from varwg.time_series_analysis import owens
# warn = False
# except Exception:
# warn = True
# if warn:
# warnings.warn("""Could not import owens t function.
# Try (if on linux):
# wget {url}
# f2py -c -m owens owens.f90
# cp owens.so {src_dir}/
# """.format(url=url, src_dir=src_dir))
# return False
# return True
# else:
# return False
# try:
# from varwg.time_series_analysis import owens
# except ImportError:
# owens = _build_owens()
# if owens:
# # it might have been build freshly
# from varwg.time_series_analysis import owens
# owens_t = np.vectorize(lambda h, a: owens.t(h, a))
# some special functions vectorized to be able to handle arrays as input
gamma_func = np.vectorize(lambda x: special.gamma(x))
gammaln = np.vectorize(lambda x: special.gammaln(x))
gammainc = np.vectorize(lambda a, x: special.gammainc(a, x))
gammaincinv = np.vectorize(lambda a, qq: special.gammaincinv(a, qq))
digamma = np.vectorize(lambda a: special.digamma(a))
hyp1f1 = np.vectorize(lambda a, b, x: special.hyp1f1(a, b, x))
nctdtr = np.vectorize(lambda x1, x2, x3: special.nctdtr(x1, x2, x3))
nctdtrit = np.vectorize(lambda x1, x2, x3: special.nctdtrit(x1, x2, x3))
stdtr = np.vectorize(lambda df, x: special.stdtr(df, x))
stdtrit = np.vectorize(lambda df, qq: special.stdtrit(df, qq))
sqrt2 = np.sqrt(2)
[docs]
def max_likelihood(
density_func,
x0,
values=None,
opt_func=None,
disp=False,
weights=None,
method="BFGS",
bounds=None,
*args,
**kwds,
):
"""Fits parameters of a density function to data, using the log-maximum-
likelihood approach."""
if opt_func is None:
opt_func = sp_optimize.minimize
if weights is None and values is not None:
weights = np.ones_like(values)
# i hereby define the unlikelihood function as the negative
# log-likelihood function, so that i can maximize the likelihood
# by minimizing the unlikelihood
if values is None:
args = ()
def unlikelihood(params):
return -np.sum(np.log(weights * density_func(params) + 1e-12))
else:
args = (values,)
def unlikelihood(params, values):
# we have to express fixed values in terms of kwds,
# because there could be collisions otherwise...
# first come, first serve: interpret the first n elements
# of params as the first n variables in
# dist.parameter_names
p_kwds = {
par_name: value
for par_name, value in zip(
density_func.__self__.parameter_names, params
)
if par_name not in kwds
}
p_kwds.update(kwds)
dens = density_func(values, **p_kwds)
try:
inside = weights * dens
except FloatingPointError:
inside = dens
inside[~np.isfinite(inside)] = 1e-9
inside[inside <= 0] = 1e-9
return -np.sum(np.log(inside))
return opt_func(
unlikelihood,
x0,
args=args,
# bounds=density_func._bounds,
bounds=bounds,
# method=method,
# method=("L-BFGS-B" if density_func._bounds else method),
method=("L-BFGS-B" if bounds else method),
options={"maxiter": 1e4 * len(x0), "disp": disp},
)
[docs]
def min_ks(cdf, x0, values, opt_func=None, disp=False, *args, **kwds):
"""Fits parameters of a cumulative distribution function to data,
minimizing the Kolmogorov-Smirnof test statistic."""
if opt_func is None:
opt_func = sp_optimize.fmin
n = len(values)
values_sorted = np.sort(values)
ranks_plus = np.arange(0.0, n) / n
ranks_minus = np.arange(1.0, n + 1) / n
def ks(params):
# we have to express fixed values in terms of kwds, because there
# could be collisions otherwise...
# first come, first serve: interpret the first n elements of params
# as the first n variables in dist.parameter_names
p_kwds = {
par_name: value
for par_name, value in zip(cdf.__self__.parameter_names, params)
if par_name not in kwds
}
p_kwds.update(kwds)
cdf_values = cdf(values_sorted, **p_kwds)
cdf_values[np.isnan(cdf_values)] = np.inf
dmin_plus = np.abs(cdf_values - ranks_plus).max()
dmin_minus = np.abs(cdf_values - ranks_minus).max()
return max(dmin_plus, dmin_minus)
return opt_func(
ks,
x0,
args=args,
maxiter=1e4 * len(x0),
maxfun=1e4 * len(x0),
disp=disp,
)
[docs]
def min_fsum(cdf, x0, values, opt_func=None, disp=False, *args, **kwds):
if opt_func is None:
opt_func = sp_optimize.fmin
n = len(values)
values_sorted = np.sort(values)
ranks_plus = np.arange(0.0, n) / n
ranks_minus = np.arange(1.0, n + 1) / n
def ks(params):
# we have to express fixed values in terms of kwds, because there
# could be collisions otherwise...
# first come, first serve: interpret the first n elements of params
# as the first n variables in dist.parameter_names
p_kwds = {
par_name: value
for par_name, value in zip(cdf.__self__.parameter_names, params)
if par_name not in kwds
}
p_kwds.update(kwds)
cdf_values = cdf(values_sorted, **p_kwds)
cdf_values[np.isnan(cdf_values)] = np.inf
dmin_plus = ((cdf_values - ranks_plus) ** 2).sum()
dmin_minus = ((cdf_values - ranks_minus) ** 2).sum()
return dmin_plus + dmin_minus
return opt_func(
ks,
x0,
args=args,
maxiter=1e4 * len(x0),
maxfun=1e4 * len(x0),
disp=disp,
)
[docs]
class Dist(metaclass=DistMeta):
"""Mimics part of the interface of stats.rv_continuous. Comes with a few
extra goodies."""
# i would like to call the first item "solution", but in order to
# have a nice compatibility to the optime.minimize result object,
# it has to be "x"
Result = namedtuple("Result", ("x", "supplements", "success"))
isscipy = False
@abstractmethod
def _pdf(self):
pass
@abstractmethod
def _cdf(self):
pass
@abstractmethod
def _ppf(self):
pass
@abstractmethod
def _fit(self):
pass
@property
def scipy_(self):
return False
def _clean_kwds(self, kwds):
"""Return a copy with only the keywords that are also in
self.parameter_names."""
return {
key: value
for key, value in list(kwds.items())
if key in self.parameter_names
}
[docs]
def fit(self, *args, **kwds):
return self._fit(*args, **self._clean_kwds(kwds))
[docs]
def sample(self, *args, **kwds):
size = np.atleast_1d(args[0])
qq = varwg.rng.random(size)
return self.ppf(qq, *args[1:], **self._clean_kwds(kwds))
[docs]
@my.asscalar
def pdf(self, *args, **kwds):
densities = np.atleast_1d(self._pdf(*args, **self._clean_kwds(kwds)))
if "x" in kwds:
x = kwds.pop("x")
invalid_x = self._invalid_x(x, *args, **kwds)
else:
invalid_x = self._invalid_x(args[0], *args[1:], **kwds)
return np.where(invalid_x | np.isinf(densities), np.nan, densities)
[docs]
@my.asscalar
def cdf(self, *args, **kwds):
qq = np.atleast_1d(self._cdf(*args, **self._clean_kwds(kwds)))
if "x" in kwds:
x = kwds.pop("x")
invalid_x = self._invalid_x(x, *args, **kwds)
else:
invalid_x = self._invalid_x(args[0], *args[1:], **kwds)
return np.where(invalid_x | np.isinf(qq), np.nan, qq)
[docs]
@my.asscalar
def ppf(self, *args, **kwds):
if "qq" in kwds:
qq = kwds.pop("qq")
else:
qq = args[0]
qq = np.atleast_1d(qq)
finite_mask = np.isfinite(qq)
quantiles_finite = np.where(finite_mask, qq, -1)
x = np.atleast_1d(
self._ppf(quantiles_finite, *args[1:], **self._clean_kwds(kwds))
)
if quantiles_finite.size == 1:
quantiles_finite = np.full_like(x, quantiles_finite)
x = self._fix_x(x)
x[(quantiles_finite < 0) | (quantiles_finite > 1)] = np.nan
return x
[docs]
@my.asscalar
def mean(self, *args, **kwds):
"""Crude estimation of the expected value. Heavy tails break this!"""
epsilon = 1e-6
x_min, x_max = self.ppf([epsilon, 1 - epsilon], *args, **kwds)
result = quad(lambda x: x * self.pdf(x, *args, **kwds), x_min, x_max)
return result[0]
def __call__(self, *params):
return Frozen(self, *params)
[docs]
def fit_ml(self, values, x0=None, *args, **kwds):
if x0 is None:
x0 = (1,) * self.n_pars
result = max_likelihood(
self.pdf, x0, values, bounds=self._bounds, *args, **kwds
)
return self.Result(
x=result.x, supplements=None, success=result.success
)
# def fit_ml(self, values, x0=None, *args, **kwds):
# if x0 is None:
# x0 = (1, ) * len(self.parameter_names)
# if self.supplements_names is not None:
# x0 = list(x0)
# for name in self.supplements_names:
# x0[self.parameter_names.index(name)] = None
# # x0.pop(self.parameter_names.index(name))
# params = max_likelihood(self.pdf, x0, values, *args, **kwds)
# if self.supplements_names is not None:
# supplements = {name: params[self.parameter_names.index(name)]
# for name in self.supplements_names}
# for name in self.supplements_names:
# params.remove(name)
# else:
# supplements = self.supplements
# return self.Result(x=params, supplements=supplements)
[docs]
def fit_ks(
self, values, opt_func=sp_optimize.fmin, x0=None, *args, **kwds
):
if x0 is None:
x0 = (1,) * len(self.parameter_names)
return min_ks(self.cdf, x0, values, opt_func, *args, **kwds)
[docs]
def fit_fsum(
self, values, opt_func=sp_optimize.fmin, x0=None, *args, **kwds
):
if x0 is None:
x0 = (1,) * len(self.parameter_names)
return min_fsum(self.cdf, x0, values, opt_func, *args, **kwds)
def _constraints(self, *args, **kwds):
param = kwds[list(kwds.keys())[0]]
return np.full_like(param, False, dtype=bool)
def _invalid_x(self, x, *args, **kwds):
"""Returns a mask indicating where x < lower bound or x > upper bound."""
x = np.atleast_1d(x)
lower_key = (
"lc"
if "lc" in self.parameter_names
else "l" if "l" in self.parameter_names else None
)
lower = kwds.get(lower_key, False)
b_shape = np.broadcast(x, lower).shape
mask = np.atleast_1d(np.full(b_shape, False, dtype=bool))
if isinstance(lower, np.ndarray) or lower:
if lower is None:
lower_i = self.parameter_names.index(lower_key)
lower = np.atleast_1d(args[lower_i])
if x.size == 1 and isinstance(lower, np.ndarray):
x = np.full_like(lower, x)
mask = np.full_like(lower, mask)
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
mask[x < lower] = True
upper_key = (
"uc"
if "uc" in self.parameter_names
else "u" if "u" in self.parameter_names else None
)
upper = kwds.get(upper_key, False)
if isinstance(upper, np.ndarray) or upper:
if upper is None:
upper_i = self.parameter_names.index(upper_key)
upper = np.atleast_1d(args[upper_i])
if x.size == 1 and isinstance(upper, np.ndarray):
x = np.full_like(upper, x)
mask = np.full_like(upper, mask)
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
mask[x > upper] = True
return np.squeeze(mask)
def _fix_x(self, x, *args, **kwds):
"""The distributions might know what to replace invalid x-values with."""
return x
[docs]
class Frozen(object):
# this idea is stolen from the scipy.stats-module and expanded
def __init__(self, dist, *params):
self.dist = dist
self.params = params
self.parameter_names = dist.parameter_names
self.name = dist.name
[docs]
def pdf(self, x):
return self.dist.pdf(x, **self.parameter_dict)
[docs]
def cdf(self, x):
return self.dist.cdf(x, **self.parameter_dict)
[docs]
def ppf(self, qq):
return self.dist.ppf(qq, **self.parameter_dict)
[docs]
def sample(self, size):
return self.dist.sample(size, **self.parameter_dict)
@property
def mean(self):
"""Crude estimation of the expected value. Heavy tails break this!"""
return self.dist.mean(**self.parameter_dict)
@property
def median(self):
return self.ppf(0.5)
@property
def parameter_dict(self):
return {
name: value
for name, value in zip(self.parameter_names, self.params)
}
[docs]
def plot_fit(self, values=None, n_classes=30):
"""Display a combined plot of a histogram, fitted pdf, empirical cdf
and fitted cdf."""
if values is None and hasattr(self.dist, "x"):
values = self.dist.x
fig = plt.figure()
ax1 = fig.add_subplot(111)
# the histogram of the data
bins = ax1.hist(
values, n_classes, density=True, facecolor="green", alpha=0.75
)[1]
class_middles = 0.5 * (bins[1:] + bins[:-1])
density = self.pdf(class_middles)
ax1.plot(class_middles, density, "r--")
# the quantile part
ax2 = ax1.twinx()
# empirical cdf
values_sort = np.sort(values)
ranks_emp = (0.5 + np.arange(len(values))) / len(values)
ax2.plot(values_sort, ranks_emp)
# theoretical cdf
x = np.linspace(values.min(), values.max(), 5e2)
ranks_theory = self.cdf(x)
ax2.plot(x, ranks_theory, "r--")
ax2.grid()
plt.title(
", ".join(
" %s: %.3f" % (par_name, par)
for par_name, par in self.parameter_dict.items()
if par_name != "kernel_data"
)
)
fig.suptitle(self.dist.name)
return fig, (ax1, ax2)
[docs]
def plot_qq(self, values, *args, **kwds):
"""A qq-plot. Scatters theoretic over empirical ranks."""
ranks_emp = (stats.stats.rankdata(values) - 0.5) / len(values)
ranks_the = self.cdf(values)
fig = plt.figure()
plt.scatter(ranks_emp, ranks_the, marker="o", *args, **kwds)
plt.plot([0, 1], [0, 1])
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Empirical ranks")
plt.ylabel("Fitted ranks")
plt.grid()
return fig
[docs]
class Censored(Dist):
def __init__(self, distribution, lc=-np.inf, uc=np.inf):
self.dist = distribution
self.lc = lc
self.uc = uc
# read "lc" as "lower cut", "uc" as "upper cut"
# mind you: these are different than lower and upper bound in i.e. the
# beta distribution. we can also truncate the beta distribution and
# have upper bounds above the upper cut.
if hasattr(distribution, "parameter_names"):
self.parameter_names = tuple(
list(distribution.parameter_names) + ["lc", "uc"]
)
self.name = "censored " + distribution.name
def _pdf(self, x, *args, **kwds):
x = np.atleast_1d(x)
args = [np.asarray(arg) for arg in args]
lc = self.lc if "lc" not in kwds else kwds.pop("lc")
uc = self.uc if "uc" not in kwds else kwds.pop("uc")
density = self.dist.pdf(x, *args, **kwds)
qq = self.dist.cdf(x, *args, **kwds)
lc_above_ii = np.atleast_1d(x > lc)
if np.any(lc_above_ii):
try:
lc_index = np.argmin(x[lc_above_ii])
density[lc_index] += qq[lc_index]
except IndexError:
for jj, lc_above in enumerate(lc_above_ii):
lc_index = np.argmin(x[0, lc_above])
density[jj, lc_index] += qq[jj, lc_index]
uc_below_ii = np.atleast_1d(x < uc)
if np.any(uc_below_ii):
try:
uc_index = np.argmax(x[uc_below_ii])
density[uc_index] += 1 - qq[uc_index]
except IndexError:
for jj, uc_below in enumerate(uc_below_ii):
uc_index = np.argmax(x[0, uc_below])
density[jj, uc_index] += 1 - qq[jj, uc_index]
density[(x < lc) | (x > uc)] = 0
return density
def _cdf(self, x, *args, **kwds):
x = np.asarray(x)
args = [np.asarray(arg) for arg in args]
lc = self.lc if "lc" not in kwds else kwds.pop("lc")
uc = self.uc if "uc" not in kwds else kwds.pop("uc")
qq = np.atleast_1d(self.dist.cdf(x, *args, **kwds))
below_lc_ii = np.atleast_1d(x <= lc)
if np.any(below_lc_ii):
lc_quantile = self.dist.cdf(lc, *args, **kwds)
if np.isscalar(lc_quantile):
lc_quantile = np.full_like(x, lc_quantile)
qq[below_lc_ii] = lc_quantile[below_lc_ii]
above_uc_ii = np.atleast_1d(x >= uc)
if np.any(above_uc_ii):
uc_quantile = self.dist.cdf(uc, *args, **kwds)
if np.isscalar(uc_quantile):
uc_quantile = np.full_like(x, uc_quantile)
qq[above_uc_ii] = 1 - uc_quantile[above_uc_ii]
qq[qq < 0] = 0
qq[qq > 1] = 1
return np.squeeze(qq)
def _ppf(self, qq, *args, **kwds):
qq = np.asarray(qq)
args = [np.asarray(arg) for arg in args]
lc = self.lc if "lc" not in kwds else kwds["lc"] # kwds.pop("lc")
uc = self.uc if "uc" not in kwds else kwds["uc"] # kwds.pop("uc")
x = self.dist.ppf(qq, *args, **kwds)
lower_ii = x < lc
if np.any(lower_ii):
x[lower_ii] = (np.full_like(x, lc) if np.isscalar(lc) else lc)[
lower_ii
]
upper_ii = x > uc
if np.any(upper_ii):
x[upper_ii] = (np.full_like(x, uc) if np.isscalar(uc) else uc)[
upper_ii
]
return x
def _fit(self, x, *args, **kwds):
# try:
# x0 = self.dist._feasible_start
# except AttributeError:
# x0 = self.dist.fit(x, *args, **kwds)
# using the cdf to fit, because of the insanity inherent in the pdf
# (i.e. no maximum likelihood)
# return min_ks(self.cdf, x0, x, lc=self.lc, uc=self.uc, *args)
# fit on valid x
valid_x = x[~self._invalid_x(x, *args, **kwds)]
return self.dist.fit(
valid_x,
*args,
**{
key: val
for key, val in list(kwds.items())
if key != "lc" and key != "uc"
},
)
[docs]
def constraints(self, x, *args, **kwds):
return self.dist.constraints(x, *args, **kwds)
[docs]
class Truncated(Dist):
"""Wraps around any other distribution found here to truncate its upper
and/or lower tail."""
def __init__(self, distribution, lc=-np.inf, uc=np.inf):
self.dist = distribution
self.lc = lc
self.uc = uc
# read "lc" as "lower cut", "uc" as "upper cut"
# mind you: these are different than lower and upper bound in i.e. the
# beta distribution. we can also truncate the beta distribution and
# have upper bounds above the upper cut.
self.parameter_names = tuple(
list(distribution.parameter_names) + ["lc", "uc"]
)
self.name = "truncated " + distribution.name
def _pdf(self, x, *args, **kwds):
args = [np.asarray(arg) for arg in args]
if "uc" in kwds:
uc = kwds.pop("uc")
if len(args) > 0:
del args[-1]
else:
uc = self.uc
if "lc" in kwds:
lc = kwds.pop("lc")
if len(args) > 0:
del args[-1]
else:
lc = self.lc
if len(args) > len(self.dist.parameter_names):
args = args[: len(self.dist.parameter_names)]
# still sober
un_trunc = self.dist.pdf(x, *args, **kwds) / (
self.dist.cdf(uc, *args, **kwds) - self.dist.cdf(lc, *args, **kwds)
)
un_trunc[(x < lc) | (x > uc)] = 1e-9
# now totally drunkated
return un_trunc
def _cdf(self, x, *args, **kwds):
args = [np.asarray(arg) for arg in args]
if "uc" in kwds:
uc = kwds.pop("uc")
if len(args) > 0:
del args[-1]
else:
uc = self.uc
if "lc" in kwds:
lc = kwds.pop("lc")
if len(args) > 0:
del args[-1]
else:
lc = self.lc
if len(args) > len(self.dist.parameter_names):
args = args[: len(self.dist.parameter_names)]
qq = (
self.dist.cdf(x, *args, **kwds) - self.dist.cdf(lc, *args, **kwds)
) / (
self.dist.cdf(uc, *args, **kwds) - self.dist.cdf(lc, *args, **kwds)
)
qq[qq < 0] = 0
qq[qq > 1] = 1
return qq
def _ppf(self, qq, *args, **kwds):
args = [np.asarray(arg) for arg in args]
if "uc" in kwds:
uc = kwds.pop("uc")
if len(args) > 0:
del args[-1]
else:
uc = self.uc
if "lc" in kwds:
lc = kwds.pop("lc")
if len(args) > 0:
del args[-1]
else:
lc = self.lc
if len(args) > len(self.dist.parameter_names):
args = args[: len(self.dist.parameter_names)]
x = self.dist.ppf(
qq * self.dist.cdf(uc, *args, **kwds)
+ (1 - qq) * self.dist.cdf(lc, *args, **kwds),
*args,
**kwds,
)
lower_ii = x < lc
if np.sum(lower_ii) > 0:
x[lower_ii] = lc[lower_ii]
upper_ii = x > uc
if np.sum(upper_ii) > 1:
x[upper_ii] = uc[upper_ii]
return x
def _fit(self, x, *args, **kwds):
try:
x0 = self.dist._feasible_start
except AttributeError:
x0 = self.dist.fit(x, *args)
# using the cdf to fit, because of the insanity inherent in the pdf
# (i.e. no maximum likelihood)
if self.lc is None:
x0 = tuple(list(x0) + [np.min(x)])
kwds["lc"] = x0[-1]
else:
kwds["lc"] = self.lc
if self.uc is None:
x0 = tuple(list(x0) + [np.max(x)])
kwds["uc"] = x0[-1]
else:
kwds["uc"] = self.uc
solution = list(min_ks(self.cdf, x0, x, *args, **kwds))
if self.uc is None:
self.uc = solution.pop(-1)
if self.lc is None:
self.lc = solution.pop(-1)
if "lc" in kwds:
solution += [kwds["lc"]]
if "uc" in kwds:
solution += [kwds["uc"]]
return solution
[docs]
def constraints(self, x, *args, **kwds):
return self.dist.constraints(x, *args, **kwds)
[docs]
class Normal(Dist):
_feasible_start = (0, 2)
def _pdf(self, x, mu=0, sigma=1):
x, mu, sigma = np.atleast_1d(x, mu, sigma)
if NE:
_pi = np.pi
dens = ne.evaluate(
"(1. / (2 * _pi * sigma ** 2) ** .5 *"
+ "exp(-(x - mu) ** 2 / (2 * sigma ** 2)))"
)
else:
dens = (
1.0
/ (2 * np.pi * sigma**2) ** 0.5
* np.exp(-((x - mu) ** 2) / (2 * sigma**2))
)
return dens
def _cdf(self, x, mu=0, sigma=1):
x, mu, sigma = np.atleast_1d(x, mu, sigma)
qq = 0.5 * (1 + special.erf((x - mu) / (sigma * 2**0.5)))
return qq
def _ppf(self, qq, mu=0, sigma=1):
qq, mu, sigma = np.atleast_1d(qq, mu, sigma)
x = special.ndtri(qq) * sigma + mu
return x
def _fit(self, x, mu=None, sigma=None):
x_noninf = x[np.isfinite(x)]
if mu is None:
mu = x_noninf.mean()
if sigma is None:
sigma = x_noninf.std()
return mu, sigma
norm = Normal()
# class SkewNormal(Dist):
# """This is taken and adapted from
# http://promanthan.com/randomstuff/skewt-0.0.1.tgz
# which might end up in scipy.stats"""
# @staticmethod
# def _fui(h, i):
# return (h ** (2 * i)) / ((2 ** i) * gamma(i + 1))
#
# @staticmethod
# def _tInt(h, a, jmax, cutPoint):
# seriesL = np.empty(0)
# seriesH = np.empty(0)
# i = np.arange(0, jmax + 1)
# low = h <= cutPoint
# hL = h[low]
# hH = h[np.logical_not(low)]
# L = hL.size
# if L > 0:
# b = SkewNormal._fui(hL[:, np.newaxis], i)
# cumb = b.cumsum(axis=1) # transposed compared to R code
# b1 = np.exp(-0.5 * hL ** 2)[:, np.newaxis] * cumb
# matr = np.ones((jmax + 1, L)) - b1.transpose()
# jk = ([1.0, -1.0] * jmax)[0:jmax + 1] / (2 * i + 1)
# matr = np.inner((jk[:, np.newaxis] * matr).transpose(),
# a ** (2 * i + 1.0))
# seriesL = (np.arctan(a) - matr.flatten(1)) / (2 * np.pi)
# if hH.size > 0:
# seriesH = (np.arctan(a) *
# np.exp(-0.5 * (hH ** 2.0) * a / np.arctan(a))
# * (1 + 0.00868 * (hH ** 4.0) * a ** 4.0) /
# (2.0 * np.pi))
# series = np.empty(h.size)
# series[low] = seriesL
# series[np.logical_not(low)] = seriesH
# return series
#
# @staticmethod
# def _tOwen(h, a, jmax=50, cutPoint=6):
# aa = np.abs(a)
# ah = np.abs(h)
# if np.isnan(aa):
# raise ValueError("a is NaN")
# if np.isposinf(aa):
# return 0.5 * norm().cdf(-ah)
# if aa == 0.0:
# return np.zeros(h.size)
# na = np.isnan(h)
# inf = np.isposinf(ah)
# ah[np.logical_or(na, inf)] = 0
# if aa <= 1:
# owen = SkewNormal._tInt(ah, aa, jmax, cutPoint)
# else:
# owen = (0.5 * norm().cdf(ah) + norm().cdf(aa * ah)
# * (0.5 - norm().cdf(ah)) -
# SkewNormal._tInt(aa * ah, (1.0 / aa), jmax, cutPoint))
# owen[np.isposinf(owen)] = 0
# return owen * np.sign(a)
#
# def _pdf(self, x, zeta, omega, alpha):
# x, zeta, omega, alpha = np.atleast_1d(x, zeta, omega, alpha)
# return (1 / (omega * np.pi) *
# np.exp(-(x - zeta) ** 2 / (2 * omega ** 2))
# * norm.cdf(alpha * (x - zeta) / omega))
#
# def _cdf(self, x, zeta, omega, alpha):
# x, zeta, omega, alpha = np.atleast_1d(x, zeta, omega, alpha)
# return (norm.cdf((x - zeta) / omega) -
# 2 * SkewNormal._tOwen((x - zeta) / omega, alpha))
#
# def _ppfInternal(self, qq, shape):
# maxQ = np.sqrt(chi2.ppf(qq, 1))
# minQ = -np.sqrt(chi2.ppf(1 - qq, 1))
# if shape > 1e+5:
# return maxQ
# if shape < -1e+5:
# return minQ
# nan = np.isnan(qq) | (qq > 1) | (qq < 0)
# zero = qq == 0
# one = qq == 1
# qq[nan | zero | one] = 0.5
# cum = SkewNormal._cumulants(shape, 4)
# g1 = cum[2] / cum[1] ** (3 / 2.0)
# g2 = cum[3] / cum[1] ** 2
# x = norm().ppf(qq)
# x = (x + (x ** 2 - 1) * g1 / 6 + x * (x ** 2 - 3) * g2 / 24 -
# x * (2 * x ** 2 - 5) * g1 ** 2 / 36)
# x = cum[0] + np.sqrt(cum[1]) * x
# tol = 1e-8
# maxErr = 1
# while maxErr > tol:
# sn = skewnorm(shape)
# x1 = x - (sn.cdf(x) - qq) / (sn.pdf(x))
# x1 = np.minimum(x1, maxQ)
# x1 = np.maximum(x1, minQ)
# maxErr = np.amax(np.abs(x1 - x) / (1 + np.abs(x)))
# x = x1
# x[nan] = np.NaN
# x[zero] = -np.Infinity
# x[one] = np.Infinity
# return x
#
# def _ppf(self, qq, shape):
# if np.all(shape == shape[0]):
# return self._ppfInternal(qq, shape[0])
# else:
# vec = np.vectorize(lambda qq, shape:
# self._ppfInternal(np.array([qq]), shape))
# return vec(qq, shape)
#
# def _ppf(self, qq, zeta, omega, alpha):
# qq, zeta, omega, alpha = \
# np.atleast_1d(qq, zeta, omega, alpha)
# return self._ppf(qq, alpha) * omega + zeta
# class SkewNormal(Dist):
# _feasible_start = (0, 1, .5)
# def _pdf(self, x, zeta=0, omega=1, alpha=0):
# x, zeta, omega, alpha = np.atleast_1d(x, zeta, omega, alpha)
# t = (x - zeta) / omega
# dens = 2. / omega * norm.pdf(t) * norm.cdf(t * alpha)
# return dens
# def _cdf(self, x, zeta=0, omega=1, alpha=0):
# x, zeta, omega, alpha = np.atleast_1d(x, zeta, omega, alpha)
# t = (x - zeta) / omega
# qq = norm.cdf(t) - 2. * owens_t(t, alpha)
# return qq
# def _ppf(self, qq, zeta=0, omega=1, alpha=0):
# qq, zeta, omega, alpha = \
# np.atleast_1d(qq, zeta, omega, alpha)
# if len(zeta) == 1:
# zeta = np.full_like(qq, zeta)
# if len(omega) == 1:
# omega = np.full_like(qq, omega)
# if len(alpha) == 1:
# alpha = np.full_like(qq, alpha)
# qq0 = norm.cdf(qq) * omega + zeta
# x = np.empty_like(qq)
# for i in range(len(qq0)):
# q_exp = qq[i]
# def error(x):
# q_act = self.cdf(x, zeta[i], omega[i], alpha[i])
# return (q_act - q_exp) ** 2
# x[i] = sp_optimize.minimize(error, qq0[i],
# method="Nelder-Mead"
# )["x"]
# return x
# def _fit(self, x, skew_max=.9):
# skew = stats.skew(x)
# skew_23 = np.abs(skew)**(2. / 3)
# if not (-skew_max < skew < skew_max):
# delta = np.sqrt(np.pi / 2 * skew_max**(2. / 3) /
# (skew_max**(2. / 3) + (
# (4 - np.pi) / 2)**(2. / 3)))
# delta = np.copysign(delta, skew)
# # warnings.warn("Sample skew %.2f not in feasible range ~(-1,1)!" %
# # skew)
# else:
# delta = np.sqrt(np.pi / 2 * skew_23 / (skew_23 + (
# (4 - np.pi) / 2)**(2. / 3)))
# delta = np.copysign(delta, skew)
# alpha = delta / np.sqrt(1 - delta**2)
# omega = np.sqrt(np.var(x) / (1 - 2 * delta**2 / np.pi))
# zeta = np.mean(x) - omega * np.sqrt(2 / np.pi) * delta
# return self.fit_ml(x, x0=(zeta, omega, alpha)).x
# if not (-skew_max < skew < skew_max):
# zeta, omega, alpha = self.fit_ml(
# x, x0=(zeta, omega, alpha / 2))
# delta = alpha / np.sqrt(1 + alpha**2)
# return zeta, omega, alpha
# if owens:
# skewnorm = SkewNormal()
# class ExGauss(Dist):
# _feasible_start = (0, 1, 1)
# _bounds = [(-np.inf, np.inf),
# (1e-9, np.inf),
# (1e-9, np.inf)]
# def _pdf(self, x, mu, sigma, gamma):
# x, mu, sigma, gamma = np.atleast_1d(x, mu, sigma, gamma)
# gamma_sigma_sqr = gamma * sigma ** 2
# return (gamma / 2 *
# np.exp(gamma / 2 * (2 * mu + gamma_sigma_sqr - 2 * x)) *
# special.erfc((mu + gamma_sigma_sqr - x) / (sqrt2 * sigma)))
# def _cdf(self, x, mu, sigma, gamma):
# x, mu, sigma, gamma = np.atleast_1d(x, mu, sigma, gamma)
# return norm.cdf(x, mu, sigma) - self._pdf(x, mu, sigma, gamma) / gamma
# def _ppf(self, q, mu, sigma, gamma):
# q, mu, sigma, gamma = np.atleast_1d(q, mu, sigma, gamma)
[docs]
class TruncatedNormal(Dist):
_feasible_start = (0, 1, -1, 1)
_additional_kwds = {"lc": -1, "uc": 1}
_bounds = [(-np.inf, np.inf), (1e-9, np.inf)]
def _pdf(self, x, mu=0, sigma=1, lc=-np.inf, uc=np.inf):
x, mu, sigma, lc, uc = np.atleast_1d(x, mu, sigma, lc, uc)
if (x.shape != lc.shape) or (x.shape != uc.shape):
x, mu, sigma, lc, uc = np.broadcast_arrays(x, mu, sigma, lc, uc)
un_trunc = np.atleast_1d(
norm.pdf(x, mu, sigma)
/ (norm.cdf(uc, mu, sigma) - norm.cdf(lc, mu, sigma))
)
un_trunc[(x < lc) | (x > uc)] = 0
sigma_neg_mask = sigma <= 0
if np.any(sigma_neg_mask):
if sigma.size == 1:
un_trunc += sigma**2
else:
un_trunc[sigma_neg_mask] += sigma[sigma_neg_mask] ** 2
return un_trunc
def _cdf(self, x, mu=0, sigma=1, lc=-np.inf, uc=np.inf):
x, mu, sigma, lc, uc = np.atleast_1d(x, mu, sigma, lc, uc)
if (x.shape != lc.shape) or (x.shape != uc.shape):
x, mu, sigma, lc, uc = np.broadcast_arrays(x, mu, sigma, lc, uc)
qq = (norm.cdf(x, mu, sigma) - norm.cdf(lc, mu, sigma)) / (
norm.cdf(uc, mu, sigma) - norm.cdf(lc, mu, sigma)
)
qq = np.atleast_1d(qq)
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
qq[qq < 0] = 0
qq[qq > 1] = 1
return qq
def _ppf(self, qq, mu=0, sigma=1, lc=-np.inf, uc=np.inf):
qq, mu, sigma, lc, uc = np.atleast_1d(qq, mu, sigma, lc, uc)
if (qq.shape != lc.shape) or (qq.shape != uc.shape):
qq, mu, sigma, lc, uc = np.broadcast_arrays(qq, mu, sigma, lc, uc)
x = norm.ppf(
qq * norm.cdf(uc, mu, sigma) + (1 - qq) * norm.cdf(lc, mu, sigma),
mu,
sigma,
)
x = np.atleast_1d(x)
lower_ii = x < lc
if np.sum(lower_ii) > 0:
x[lower_ii] = lc[lower_ii]
upper_ii = x > uc
if np.sum(upper_ii) > 1:
x[upper_ii] = uc[upper_ii]
return x
def _fit(self, x, lc=-np.inf, uc=np.inf):
x0 = [np.nanmean(x), np.nanstd(x)]
kwds = {}
if lc is None: # or np.isneginf(lc):
x0 += [np.nanmin(x)]
else:
kwds["lc"] = lc
if uc is None: # or np.isinf(uc):
x0 += [np.nanmax(x)]
else:
kwds["uc"] = uc
return self.fit_ml(x, x0=x0, method="Nelder-Mead", **kwds).x
def _constraints(self, x, mu, sigma, lc=-np.inf, uc=np.inf, **kwds):
mask = (x > uc) | (x < lc) | (sigma <= 0)
return mask
truncnorm = TruncatedNormal()
[docs]
class LogNormal(Dist):
_feasible_start = (2.0, 0.5)
def _pdf(self, x, mu, sigma):
return norm.pdf(np.log(x), mu, sigma) / x
def _cdf(self, x, mu, sigma):
return norm.cdf(np.log(x), mu, sigma)
def _ppf(self, qq, mu, sigma):
return np.exp(norm.ppf(qq, mu, sigma))
def _fit(self, x):
x_noninf = x[(x > 0) & np.isfinite(x)]
mu = np.mean(np.log(x_noninf))
sigma = np.std(np.log(x_noninf))
return mu, sigma
def _constraints(self, x, mu, sigma):
mask = (x <= 0) | (sigma <= 0)
return mask
lognormal = LogNormal()
[docs]
class JohnsonSU(Dist):
_feasible_start = (1, 1, 0, 1)
def _pdf(self, x, a, b, loc, scale):
x, a, b, loc, scale = np.atleast_1d(x, a, b, loc, scale)
x = (x - loc) / scale
x2 = x * x
trm = norm.pdf(a + b * np.log(x + np.sqrt(x2 + 1)))
dens = b * 1.0 / np.sqrt(x2 + 1.0) * trm / scale
return dens
def _cdf(self, x, a, b, loc, scale):
x = np.atleast_1d(x - loc) / scale
qq = np.atleast_1d(norm.cdf(a + b * np.log(x + np.sqrt(x * x + 1))))
qq[np.isneginf(x)] = 0
return qq
def _ppf(self, qq, a, b, loc, scale):
qq, a, b, loc, scale = np.atleast_1d(qq, a, b, loc, scale)
z = np.sinh((norm.ppf(qq) - a) / b)
x = z * scale + loc
return x
def _fit(self, x):
x_noinf = x[np.isfinite(x)]
loc, scale = x_noinf.mean(), x_noinf.std()
return self.fit_fsum(x_noinf, x0=(1, 1, loc, scale))
def _constraints(self, x, a, b, loc, scale):
return b <= 0
johnsonsu = JohnsonSU()
[docs]
class Cauchy(Dist):
_feasible_start = (0, 1)
def _pdf(self, x, x0, gamma):
return 1 / (np.pi * gamma * (1 + ((x - x0) / gamma) ** 2))
def _cdf(self, x, x0, gamma):
return (1 / np.pi) * np.arctan2(x - x0, gamma) + 0.5
def _ppf(self, qq, x0, gamma):
return x0 + gamma * np.tan(np.pi * (qq - 0.5))
def _fit(self, x):
x0 = np.median(stats.trimboth(x, 0.38))
inner_quart = stats.trimboth(x, 0.25)
inner_range = inner_quart.max() - inner_quart.min()
return self.fit_ml(x, x0=(x0, 0.5 * inner_range)).x
def _constraints(self, x, x0, gamma):
return gamma < 0
cauchy = Cauchy()
[docs]
class StudentT(Dist):
_feasible_start = (0, 1)
def _pdf(self, x, mu=0, df=1):
x, mu, df = np.atleast_1d(x, mu, df)
Px = (
np.exp(gammaln((df + 1) / 2.0) - gammaln(df / 2.0))
/ np.sqrt(df * np.pi)
* (1 + (x - mu) ** 2 / df) ** (-(df + 1) / 2.0)
)
return Px
def _cdf(self, x, mu=0, df=1):
qq = stdtr(df, np.atleast_1d(x) - mu)
return qq
def _ppf(self, qq, mu=0, df=1):
qq = np.atleast_1d(qq)
x = stdtrit(df, qq) + mu
x[qq == 0] = -np.inf
x[qq == 1] = np.inf
return x
def _fit(self, x):
x_noninf = x[np.isfinite(x)]
sigma = np.var(x_noninf)
df = 2 * sigma / (sigma - 1)
if df > 6:
df = 6
return self.fit_ml(
x_noninf,
x0=(np.mean(x_noninf), df if df > 0 else 0.1),
method="Nelder-Mead",
).x
def _constraints(self, x, mu, df):
mask = (df <= 0) | (df > 8)
return mask
student_t = StudentT()
[docs]
class NoncentralT(Dist):
_feasible_start = (2, 3, 0)
# stolen from scipy.stats.distributions.nct
def _pdf(self, x, df, nc, mu=0):
x, df, nc = (np.atleast_1d(var).astype(float) for var in (x, df, nc))
pdf_x = lambda x: (
(df / x)
* (
self.cdf(x * (1 + 2 / df) ** 0.5, df + 2, nc)
- self.cdf(x, df, nc)
)
)
dens = np.where(
np.abs(x - mu) <= 0.5, self._pdf_old(x, df, nc, mu), pdf_x(x - mu)
)
dens[np.isinf(x)] = 0
return dens
def _pdf_old(self, x, df, nc, mu=0):
n = df * 1.0
nc = nc * 1.0
x2 = (x - mu) ** 2
ncx2 = nc * nc * x2
fac1 = n + x2
trm1 = (
n / 2.0 * np.log(n)
+ gammaln(n + 1)
- (
n * np.log(2)
+ nc * nc / 2.0
+ (n / 2.0) * np.log(fac1)
+ gammaln(n / 2.0)
)
)
Px = np.exp(trm1)
valF = ncx2 / (2 * fac1)
trm1 = (
sqrt2
* nc
* (x - mu)
* hyp1f1(n / 2 + 1, 1.5, valF)
/ (np.asarray(fac1 * gamma_func((n + 1) / 2)))
)
trm2 = hyp1f1((n + 1) / 2, 0.5, valF) / (
np.asarray(np.sqrt(fac1) * gamma_func(n / 2 + 1))
)
Px *= trm1 + trm2
return Px
def _cdf(self, x, df, nc, mu=0):
x, df, nc, mu = (
np.atleast_1d(var).astype(float) for var in (x, df, nc, mu)
)
qq = nctdtr(df, nc, x - mu)
return qq
def _ppf(self, qq, df, nc, mu=0):
qq, df, nc = (np.atleast_1d(var).astype(float) for var in (qq, df, nc))
x = nctdtrit(df, nc, qq) + mu
x[qq == 0] = -np.inf
x[qq == 1] = np.inf
return x
def _fit(self, x):
x_noninf = x[np.isfinite(x)]
mode = stats.mode(np.round(x_noninf, 1)).mode
x0 = (1, 2 * stats.skew(x_noninf), mode)
df, nc, mu = self.fit_ml(x_noninf, x0=x0, method="Nelder-Mead").x
# if df >= 30:
# # then we are normal anyway
# df = 30
return df, nc, mu
def _constraints(self, x, df, nc, mu):
mask = (df < 0.1) | (df > 8)
return mask
noncentral_t = NoncentralT()
# TODO: implement ppf
# class LogitNormal(Dist):
# def logit(self, x):
# return np.log(x / (1 - x))
#
# def _pdf(self, x, mu, sigma, l=0, u=1):
# x_normed = (x - l) / (u - l)
# logit = self.logit(x_normed)
# return norm_pdf(logit, mu, sigma) / (x_normed * (1 - x_normed))
#
# def _cdf(self, x, mu, sigma, l=0, u=1):
# x_normed = (x - l) / (u - l)
# logit = self.logit(x_normed)
# return norm_cdf(logit, mu, sigma)
#
# def _ppf(self, qq, mu, sigma, l=0, u=1):
# raise NotImplementedError
#
# def _fit(self, x, l=None, u=None):
# l = x.min() - 1e-9 if l is None else l
# u = x.max() + 1e-9 if u is None else u
# x_normed = (x - l) / (u - l)
# logit = self.logit(x_normed)
# mu, sigma = norm_fit(logit)
# return mu, sigma, l, u
#
# def _constraints(self, x, mu, sigma, l=0, u=1):
# x, mu, sigma, l, u = (np.asarray(var) for var in (x, mu, sigma, l, u))
# if np.any(x < l):
# return False
# if np.any(x > u):
# return False
# return True
# logitnormal = LogitNormal()
# class Rayleigh(Dist):
# def _pdf(self, x, sigma):
# return (x / sigma ** 2) * np.exp(-x ** 2 / (2 * sigma ** 2))
#
# def _cdf(self, x, sigma):
# return 1 - np.exp(-x ** 2 / (2 * sigma ** 2))
#
# def _fit(self, x):
# return ((.5 * np.average(x ** 2)) ** .5,)
#
# @abstractmethod
# def _ppf(self):
# raise NotImplementedError
#
# def _constraints(self, x, sigma):
# if np.any(sigma <= 0):
# return False
# return True
# rayleigh = Rayleigh()
# TODO: implement ppf
# class RayleighLU(Dist):
# def _pdf(self, x, sigma, l, u):
# x_normed = (x - l) / (u - l)
# return (x_normed / sigma ** 2) * \
# np.exp(-x_normed ** 2 / (2 * sigma ** 2))
#
# def _cdf(self, x, sigma, l, u):
# x_normed = (x - l) / (u - l)
# return 1 - np.exp(-x_normed ** 2 / (2 * sigma ** 2))
#
# def _fit(self, x, l=None, u=None):
# l = x.min() if l is None else l
# u = x.max() if u is None else u
# x_normed = (x - l) / (u - l)
# return (.5 * np.average(x_normed ** 2)) ** .5, l, u
#
# def _ppf(self):
# raise NotImplementedError
#
# def _constraints(self, x, sigma, l, u):
# if np.any(sigma <= 0):
# return False
# if np.any(l >= u):
# return False
# return True
# rayleighlu = RayleighLU()
# TODO: fit does not work well
# class LogNormalLU(Dist):
# def _pdf(self, x, mu, sigma, l, u):
# x_normed = (x - l) / (u - l)
# return norm_pdf(np.log(x_normed + 1e-9), mu, sigma) / (u - l)
#
# def _cdf(self, x, mu, sigma, l, u):
# x_normed = (x - l) / (u - l)
# return norm_cdf(np.log(x_normed + 1e-9), mu, sigma)
#
# def _ppf(self, qq, mu, sigma, l, u):
# return np.exp(norm_ppf(qq, mu, sigma)) * (u - l) + l
#
# def _fit(self, x, l=None, u=None):
# x_noninf = x[np.isfinite(x)]
# l = x_noninf.min() if l is None else l
# u = x_noninf.max() if u is None else u
# x_normed = (x_noninf - l) / (u - l)
# mu, sigma = norm_fit(np.log(x_normed + 1e-9))
# return self.fit_ml(x_noninf, x0=(mu, sigma, l, u))
#
# def _constraints(self, x, mu, sigma, l, u):
# if np.any(l >= u):
# return False
# return True
# lognormallu = LogNormalLU()
# lognormallu._feasible_start = (0, 1, -1, 1)
# class WeibullLU(Dist):
# def _pdf(self, x, alpha, beta, l=0, u=1):
# x_normed = (x - l) / (u - l)
# # alpha and beta should be positive
# densities = alpha * beta * x_normed ** (beta - 1) * \
# np.exp(-alpha * x_normed ** beta)
# return densities / (u - l)
#
# def _cdf(self, x, alpha, beta, l=0, u=1):
# x_normed = (x - l) / (u - l)
# return 1 - np.exp(-alpha * x_normed ** beta)
#
# def _ppf(self, qq, alpha, beta, l=0, u=1):
# x = (-np.log(1 - qq) / alpha) ** (1 / beta)
# return x * (u - l) + l
#
# def _fit(self, values, beta_start=10, *args, **kwds):
# """Implements a semi-analytical method of moments. beta is found by
# minimizing errors. After that, alpha can be estimated with the help of
# beta.
# See
# http://interstat.statjournals.net/YEAR/2000/articles/0010001.norm_pdf
# """
# values_noninf = values[np.isfinite(values)]
# l = values_noninf.min()
# u = values_noninf.max()
# values_noninf = (values_noninf - l) / (u - l)
# # we need the coefficient of variation now and the mean later. so do
# # not use stats.variation in order to not calculate the mean twice.
# mean = np.mean(values_noninf)
# cv = np.std(values_noninf) / mean
#
# def cv_error(beta):
# """Squared error between the theoretical and empirical coefficient
# of variation."""
# gamma1 = special.gamma(1 + 2 / beta)
# gamma2 = special.gamma(1 + 1 / beta)
# return (cv - (gamma1 - gamma2 ** 2) ** .5 / gamma2) ** 2
# beta = sp_optimize.fmin(cv_error, [beta_start], disp=False)[0]
# alpha = (mean / special.gamma(1 / beta + 1)) ** (-beta)
# return alpha, beta, l, u
#
# def _constraints(self, x, alpha, beta, l, u):
# x, alpha, beta = (np.asarray(var) for var in (x, alpha, beta))
# if np.any(alpha <= 0):
# return False
# if np.any(beta <= 0):
# return False
# if np.any(x < l):
# return False
# if np.any(x > u):
# return False
# return True
# weibulllu = WeibullLU()
# weibulllu._feasible_start = (1, 1, -1, 1)
[docs]
class Weibull(Dist):
_feasible_start = (0.5, 2.0)
def _pdf(self, x, alpha, beta):
# alpha and beta should be positive
densities = alpha * beta * x ** (beta - 1) * np.exp(-alpha * x**beta)
return densities
def _cdf(self, x, alpha, beta):
return 1 - np.exp(-alpha * x**beta)
def _ppf(self, qq, alpha, beta):
return (-np.log(1 - qq) / alpha) ** (1.0 / beta)
def _fit(self, values, beta_start=10, *args, **kwds):
"""Implements a semi-analytical method of moments. beta is found by
minimizing errors. After that, alpha can be estimated with the help of
beta.
See
http://interstat.statjournals.net/YEAR/2000/articles/0010001.norm_pdf
"""
# we need the coefficient of variation now and the mean later. so do
# not use stats.variation in order to not calculate the mean twice.
mean = np.nanmean(values)
cv = np.nanstd(values) / mean
def cv_error(beta):
"""Squared error between the theoretical and empirical coefficient
of variation."""
gamma1 = special.gamma(1 + 2 / beta)
gamma2 = special.gamma(1 + 1 / beta)
return (cv - (gamma1 - gamma2**2) ** 0.5 / gamma2) ** 2
beta = sp_optimize.fminbound(cv_error, -1, 1e6, disp=False)
alpha = (mean / special.gamma(1 / beta + 1)) ** (-beta)
return alpha, beta
def _constraints(self, x, alpha, beta):
mask = (x < 0) | (alpha <= 0) | (beta <= 0)
return mask
weibull = Weibull()
[docs]
class Kumaraswamy(Dist):
"""Resembles the Beta distribution, but does not need a transcendental
function."""
_feasible_start = (2, 2, 0, 1)
def _pdf(self, x, a, b, l=0, u=1):
x = np.atleast_1d(x).astype(float)
if NE:
x = ne.evaluate("(x - l) / (u - l)")
dens = ne.evaluate(
"a * b * x ** (a - 1) * (1 - x ** a) ** " + "(b - 1) / (u - l)"
)
else:
x = (x - l) / (u - l)
dens = a * b * x ** (a - 1) * (1 - x**a) ** (b - 1) / (u - l)
return dens
def _cdf(self, x, a, b, l=0, u=1):
x = np.atleast_1d(x).astype(float)
if np.any((x < l) | (x > u)):
warnings.warn("Some values below lower or above upper bounds.")
x = np.where(x > u, u, x)
x = np.where(x < l, l, x)
return 1 - (1 - ((x - l) / (u - l)) ** a) ** b
def _ppf(self, qq, a, b, l=0, u=1):
qq, a, b, l, u = np.atleast_1d(qq, a, b, l, u)
if NE:
x = ne.evaluate(
"((1 - (1 - qq) ** (1 / b)) ** " + "(1 / a)) * (u - l) + l"
)
else:
x = ((1 - (1 - qq) ** (1.0 / b)) ** (1.0 / a)) * (u - l) + l
return x
def _fit(self, x, l=None, u=None):
# return beta.fit(x, l, u)
x0 = beta.fit(x, l, u)
par_bounds = [(1e-9, np.inf), (1e-9, np.inf)]
# x0 will have entries for l and u if those where None
if l is None:
par_bounds += [(-np.inf, x0[-2 if u is None else -1])]
if u is None:
par_bounds += [(x0[-1], np.inf)]
return max_likelihood(self.pdf, x0, x, bounds=par_bounds).x
def _constraints(self, x, a, b, l=0, u=1):
mask = (a < 0) | (b < 0) | (x <= l) | (x > u)
return mask
kumaraswamy = Kumaraswamy()
[docs]
class Beta(Dist):
"""Beta distribution on the interval [l, u]."""
_feasible_start = (1, 1, 0, 1)
def _pdf(self, x, alpha, beta, l=0, u=1):
x, alpha, beta, l, u = list(map(np.asarray, (x, alpha, beta, l, u)))
# putting the values into [0, 1] according to [l, u]
# we avoid the value of exactly 0 or 1 to not get numerical problems
# x_normed = (x - l) / (u - l)
# x = np.copy(x_normed)
# ll, uu = sys.float_info.min, 1 - sys.float_info.min
# x_normed = np.where(x_normed >= uu, uu, x_normed)
# x_normed = np.where(x_normed <= ll, ll, x_normed)
if NE:
x = ne.evaluate("(x - l) / (u - l)")
return (
1
/ special.beta(alpha, beta)
* ne.evaluate(
"x ** (alpha - 1) * (1 - x) ** (beta - 1) / " + "(u - l)"
)
)
else:
x = (u - l) / (u - l)
return (
1
/ special.beta(alpha, beta)
* x ** (alpha - 1)
* (1 - x) ** (beta - 1)
/ (u - l)
)
def _cdf(self, x, alpha, beta, l=0, u=1):
if np.any((x < l) | (x > u)):
warnings.warn("Some values below lower or above upper bounds.")
# note: the betainc function below really delivers the regularized
# incomplete beta function
return special.betainc(alpha, beta, (x - l) / (u - l))
def _ppf(self, qq, alpha, beta, l=0, u=1):
return (u - l) * special.betaincinv(alpha, beta, qq) + l
def _fit(self, x, l=None, u=None):
"""stolen from
http://en.wikipedia.org/wiki/Beta_distribution#Parameter_estimation."""
# remember whether l and u where given. if they were, we do
# not return them
if l is None:
l = x.min()
return_l = True
else:
return_l = False
if u is None:
u = x.max()
return_u = True
else:
return_u = False
# have to make sure that the parameters all have the same lengths
try:
ll = np.atleast_1d(np.empty_like(u))
ll[:] = l
l = ll
except (TypeError, ValueError):
pass
try:
uu = np.empty_like(l)
uu[:] = u
u = uu
except (TypeError, ValueError):
pass
xmean = (x.mean() - l) / (u - l)
xvar = x.var()
if xvar == 0:
alpha = np.full_like(u, np.nan)
beta = np.full_like(u, np.nan)
else:
xvar /= (u - l) ** 2
alpha = xmean * (xmean * (1 - xmean) / xvar - 1)
beta = (1 - xmean) * alpha / xmean
return_list = [abs(alpha), abs(beta)]
if return_l:
return_list += [l]
if return_u:
return_list += [u]
return np.squeeze(return_list)
# return self.fit_ml(x, x0=(alpha, beta, l, u))
def _constraints(self, x, alpha, beta, l=0, u=1):
mask = (alpha <= 0) | (beta <= 0) | (x < l) | (x > u)
return mask
beta = Beta()
[docs]
class Gamma(Dist):
"""Gamma distribution."""
_feasible_start = (2, 2)
def _pdf(self, x, k, theta):
x = np.atleast_1d(x)
dens = (
x ** (k - 1)
* np.exp(-x / float(theta))
/ (special.gamma(k) * theta**k)
)
return dens
def _cdf(self, x, k, theta):
x = np.atleast_1d(x)
qq = gammainc(k, x / theta) / gamma_func(k)
# qq[np.isinf(x)] = 1
return qq
def _ppf(self, qq, k, theta):
qq = np.atleast_1d(qq)
x = theta * gammaincinv(k, qq * gamma_func(k))
# x[qq == 1] = np.inf
return x
def _fit(self, x, rel_change=1e-6):
"""A maximum-likelihood estimator. See
http://en.wikipedia.org/wiki/Gamma_distribution#Parameter_estimation
"""
def psi(k):
if k >= 8:
return np.log(k) - (
1.0
+ (1.0 - (0.1 - 1.0 / (21.0 * k**2)) / k**2) / (6.0 * k)
) / (2.0 * k)
else:
return psi(k + 1) - 1.0 / k
def psi_(k):
if k >= 8:
return (
1.0
+ (
1.0
+ (1.0 - (1.0 / 5.0 - 1.0 / (7.0 * k**2)) / k**2)
/ (3.0 * k)
)
/ (2.0 * k)
) / k
else:
return psi_(k + 1) + 1.0 / k**2
xmean = x[np.isfinite(x)].mean()
# we get problems here if there are 0's in x (log(x) = -inf, which
# happens for wind)
s = np.log(xmean) - np.mean(np.log(x[np.isfinite(x) & (x > 0)]))
k_old = (3 - s + ((s - 3) ** 2 + 24 * s) ** 0.5) / (12 * s)
k_rec = lambda k: k - (np.log(k) - psi(k) - s) / (1.0 / k - psi_(k))
k_new = k_rec(k_old)
while abs(k_old - k_new) / k_old > rel_change:
k_old = k_new
k_new = k_rec(k_old)
theta = xmean / k_new
return self.fit_fsum(x, x0=(k_new, theta))
def _constraints(self, x, k, theta):
mask = (k <= 0) | (theta <= 0) | (x < 0)
return mask
gamma = Gamma()
[docs]
class Gamma1(Dist):
_feasible_start = (1, 0.5, 2)
def _pdf(self, x, a, loc, scale):
x, a, loc, scale = np.atleast_1d(x, a, loc, scale)
z = (x - loc) / scale
dens = (np.exp((a - 1) * np.log(z) - z - gammaln(a))) / scale
return dens
def _cdf(self, x, a, loc, scale):
x, a, loc, scale = np.atleast_1d(x, a, loc, scale)
z = (x - loc) / scale
qq = gammainc(a, z)
qq[np.isinf(x)] = 1
return qq
def _ppf(self, qq, a, loc, scale):
qq, a, loc, scale = np.atleast_1d(qq, a, loc, scale)
x = gammaincinv(a, qq) * scale + loc
x[qq == 1] = np.inf
return x
def _fit(self, x):
return stats.gamma.fit(x[np.isfinite(x)])
# this is stolen from sp.stats.gamma-gen and not understood
# x_noninf = x[np.isfinite(x)]
# a = 4 / stats.skew(x_noninf) ** 2
# muhat = np.array(x_noninf).mean()
# mu2hat = np.array(x_noninf).var()
# Shat = (mu2hat / a) ** .5
# Lhat = muhat - Shat * a
# return self.fit_ml(x_noninf, x0=(a, Lhat, Shat))
gamma1 = Gamma1()
[docs]
class Expon(Dist):
_feasible_start = (1e-3, 1.0 / 5)
def _pdf(self, x, x0, lambd):
return np.where(x > x0, lambd * np.exp(-lambd * (x - x0)), 0)
def _cdf(self, x, x0, lambd):
return np.where(x > x0, 1 - np.exp(-lambd * (x - x0)), 0)
def _ppf(self, qq, x0, lambd):
return -np.log(1 - qq) / lambd + x0
def _fit(self, x):
x_noninf = x[np.isfinite(x)]
x0 = np.min(x_noninf)
lambd = 1 / np.mean(x_noninf)
return x0, lambd
def _constraints(self, x, x0, lambd):
return x <= x0
expon = Expon()
[docs]
class ExponTwo(Dist):
_feasible_start = (0, 0.5, 1, 1)
def _pdf(self, x, x0, q0, lambda1, lambda2):
x = np.atleast_1d(x)
dens = np.where(
x <= x0,
q0 * expon.pdf(x0 - x, 0, lambda1),
(1 - q0) * expon.pdf(x, x0, lambda2),
)
return dens
def _cdf(self, x, x0, q0, lambda1, lambda2):
x = np.atleast_1d(x)
qq = np.where(
x <= x0,
q0 * (1 - expon.cdf(x0 - x, 0, lambda1)),
q0 + (1 - q0) * expon.cdf(x, x0, lambda2),
)
return qq
def _ppf(self, qq, x0, q0, lambda1, lambda2):
qq = np.atleast_1d(qq)
x = np.where(
qq <= q0,
x0 + np.log(qq / q0) / lambda1,
x0 - np.log(1 - (qq - q0) / (1 - q0)) / lambda2,
)
return x
def _fit(self, x):
x_noninf = x[np.isfinite(x)]
x0 = stats.mode(np.round(x_noninf, 1)).mode
lambda1 = -1 / np.mean(x_noninf[x_noninf < x0] - x0)
lambda2 = 1 / np.mean(x_noninf[x_noninf > x0] - x0)
q0 = (np.argmin(np.abs(np.sort(x_noninf) - x0)) + 0.5) / len(x)
return self.fit_ml(
x_noninf, x0=(x0, q0, lambda1, lambda2), method="Nelder-Mead"
).x
def _constraints(self, x, x0, q0, lambda1, lambda2):
mask = (lambda1 <= 0) | (lambda2 <= 0)
return mask
expon_two = ExponTwo()
[docs]
class NoncentralLaplace(Dist):
_feasible_start = (1, 0.5, 0.5)
def _pdf(self, x, x0, q0, lambd):
lambda2 = lambd * q0 / (1 - q0)
x = np.atleast_1d(x)
dens = np.where(
x <= x0,
q0 * lambd * np.exp(-lambd * (x0 - x)),
(1 - q0) * lambda2 * np.exp(-lambda2 * (x - x0)),
)
dens[np.isinf(x)] = 0
return dens
def _cdf(self, x, x0, q0, lambd):
lambda2 = lambd * q0 / (1 - q0)
x = np.atleast_1d(x)
qq = np.where(
x <= x0,
q0 * (1 - expon.cdf(x0 - x, 0, lambd)),
q0 + (1 - q0) * expon.cdf(x, x0, lambda2),
)
return qq
def _ppf(self, qq, x0, q0, lambd):
lambda2 = lambd * q0 / (1 - q0)
qq = np.atleast_1d(qq)
x = np.where(
qq <= q0,
x0 + np.log(qq / q0) / lambd,
x0 - np.log(1 - (qq - q0) / (1 - q0)) / lambda2,
)
return x
def _fit(self, x):
x0 = stats.mode(np.round(x, 1)).mode
q0 = (np.argmin(np.abs(np.sort(x) - x0)) + 0.5) / len(x)
lambd = (1 - 2 * q0) / (q0 * np.mean(x[np.isfinite(x)]))
if lambd < 0:
# this actually means that x0 was not right, so provide a feasible
# lambd and let the ml estimation figure the rest out
lambd = -1.0 / np.mean(x[(x < x0) & np.isfinite(x)] - x0)
return self.fit_ks(x, x0=(x0, q0, lambd))
def _constraints(self, x, x0, q0, lambd):
mask = (lambd <= 0) | (q0 < 0) | (q0 > 1)
return mask
noncentral_laplace = NoncentralLaplace()
[docs]
class MDFt(object):
"""Multiple degrees of freedom t distribution.
The marginals have different degrees of freedom, therefore this is
not a "true" Multivariate t distribution.
See Serban et al 2007.
"""
@staticmethod
def _marginal_pdf(x, sigma, df):
x = np.atleast_1d(np.asarray(x))
dens = (
gamma_func((df + 1) / 2.0)
/ (np.sqrt(df * np.pi) * gamma_func(df / 2.0))
* (1 + x**2 / df) ** (-(df + 1) / 2.0)
)
return dens
[docs]
@staticmethod
def loglikelihood(data, sigma, df):
if np.any(~np.isfinite(sigma)):
return np.inf
sigma_inv_sqrt = np.linalg.inv(linalg.sqrtm(sigma))
yt = np.array(
[
np.squeeze(np.asarray(sigma_inv_sqrt[i] @ data))
/ np.sqrt((dfi - 2) / dfi)
for i, dfi in enumerate(df)
]
)
marginal_pdfs = [MDFt._marginal_pdf(yt, sigma, dfi) for dfi in df]
llh = (
np.sum(np.log(np.sqrt(df - 2) / df))
+ np.sum(list(map(np.log, marginal_pdfs)))
# is that needed? it should not change the location of
# the maximum?!
- 0.5 * np.log(linalg.det(sigma))
)
return -llh
[docs]
@staticmethod
def sample(size, sigma, df):
K = sigma.shape[0]
if isinstance(df, float):
df = np.array(K * [df])
yt = np.array([student_t.sample(size, df=dfi) for dfi in df])
zt = np.sqrt((df - 2.0) / df)[:, None] * yt
return np.dot(linalg.sqrtm(sigma), zt)
[docs]
@staticmethod
def fit(data):
K = data.shape[0]
def unlikelihood(params):
sigma = fill_lower(params[:-K])
df = params[-K:]
return -MDFt.loglikelihood(data, sigma, df)
sigma = np.cov(data)
sigma_upper = sigma[np.triu_indices_from(sigma)]
df = np.array([student_t.fit(values) for values in data])
x0 = np.concatenate((sigma_upper, df))
bounds = [
(0, None) if i == j else (None, None)
for i in range(K)
for j in range(i + 1)
]
bounds += K * [(5, None)]
result = sp_optimize.minimize(
unlikelihood, x0=x0, bounds=bounds, options=dict(disp=True)
)
sigma = fill_lower(result.x[:-K])
df = result.x[-K:]
return sigma, df
[docs]
def fill_lower(sequence):
K = int(np.ceil(np.sqrt(len(sequence))))
arr = np.empty((K, K))
arr[np.triu_indices_from(arr)] = sequence
upper_i = np.triu_indices_from(arr, k=1)
lower_i = np.tril_indices_from(arr, k=-1)
arr[lower_i] = arr[upper_i]
return arr
class _Rain(Dist):
def _mask_kwds(self, mask, kwds):
return {
name: (
vals[mask]
if (hasattr(vals, "shape") and vals.shape == mask.shape)
else vals
)
for name, vals in kwds.items()
}
def _pdf(self, meta, *args, **kwds):
pass
def _cdf(self, meta, *args, **kwds):
pass
def _ppf(self, meta, *args, **kwds):
pass
class _KDE(object):
# these are not part of the solution and should not be tested on
_cache_names = (
"kernel_data",
"x_eval",
"q_kde_eval",
) # "kernel_width"
def fit_kde(self, x, x0=None, bounds=None):
if bounds is not None:
bounds = np.squeeze(bounds)
# return kde.optimal_kernel_width(x, x0=None, bounds=bounds)
return kde.silvermans_rule(x)
def _kde_integral(
self, kernel_width, kernel_data, f_thresh, upper_eval, recalc=True
):
lower_eval = max(0.95 * f_thresh, 1e-9)
if upper_eval is None or np.any(upper_eval < kernel_data):
if len(kernel_data):
upper_eval = 1.75 * (kernel_data.max() + 1.75 * kernel_width)
else:
upper_eval = 1.025 * f_thresh
x_eval_log = np.log(np.linspace(lower_eval, upper_eval, 500))
if len(kernel_data) == 0:
x_eval_log = np.concatenate((x_eval_log, [np.log(upper_eval)]))
q_kde_eval = np.ones_like(x_eval_log)
q_kde_eval[0] = 0
return q_kde_eval, x_eval_log
if kernel_width is None:
kernel_width = kde.silvermans_rule(kernel_data)
if np.isnan(kernel_width):
kernel_data += (
upper_eval * 1e-6 * varwg.rng.normal(size=len(kernel_data))
)
kernel_width = kde.silvermans_rule(kernel_data)
dens_kde_eval = kde.kernel_density(
kernel_width,
np.log(kernel_data),
eval_points=x_eval_log,
recalc=recalc,
)
while np.isclose(dens_kde_eval.sum(), 0):
# increase number of evaluation points
n_eval_points = len(x_eval_log) * 2
x_eval_log = np.log(
np.linspace(lower_eval, upper_eval, n_eval_points)
)
dens_kde_eval = kde.kernel_density(
kernel_width,
np.log(kernel_data),
eval_points=x_eval_log,
recalc=recalc,
)
assert np.all(np.isfinite(dens_kde_eval))
x_eval = np.exp(x_eval_log)
dens_kde_eval /= x_eval
q_kde_eval = cumulative_trapezoid(y=dens_kde_eval, x=x_eval, initial=0)
assert np.all(np.isfinite(q_kde_eval))
# this is ugly, but I cannot find a better solution (now)
q_kde_eval /= q_kde_eval[-1]
assert np.all(np.isfinite(q_kde_eval))
# fill only the quantile range [q_thresh, 1]
q_kde_eval = self.q_thresh + (1.0 - self.q_thresh) * q_kde_eval
q_kde_eval = np.concatenate((q_kde_eval, [1.0]))
assert np.all(np.isfinite(q_kde_eval))
x_eval_log = np.concatenate((x_eval_log, [np.log(upper_eval)]))
return q_kde_eval, x_eval_log
def _kde_cdf(self, xx, q_kde_eval=None, x_eval=None):
xx = np.atleast_1d(xx)
qq_kde = np.empty_like(xx)
if isinstance(q_kde_eval, list):
q_kde_eval = np.asarray(q_kde_eval, dtype=object)
if isinstance(x_eval, list):
x_eval = np.asarray(x_eval, dtype=object)
for i, x in enumerate(xx):
if not (x_eval_ := x_eval[i]).ndim:
x_eval_ = x_eval
if x > x_eval_[-1]:
qq_kde[i] = 1.0
continue
if x < x_eval_[0]:
qq_kde[i] = 0.0
continue
if not (q_kde_eval_ := q_kde_eval[i]).ndim:
q_kde_eval_ = q_kde_eval
q_kde_interp = interpolate.interp1d(
x_eval_, q_kde_eval_, kind="linear", assume_sorted=True
)
qq_kde[i] = q_kde_interp(x)
return qq_kde
def _kde_ppf(self, qq, q_kde_eval=None, x_eval=None):
qq = np.atleast_1d(qq)
# assert np.all(qq >= self.q_thresh)
xx_kde = np.empty_like(qq)
if isinstance(q_kde_eval, list):
q_kde_eval = np.asarray(q_kde_eval, dtype=object)
if isinstance(x_eval, list):
x_eval = np.asarray(x_eval, dtype=object)
for i, q in enumerate(qq):
if not (x_eval_ := x_eval[i]).ndim:
x_eval_ = x_eval
if not (q_kde_eval_ := q_kde_eval[i]).ndim:
q_kde_eval_ = q_kde_eval
x_kde_interp = interpolate.interp1d(
q_kde_eval_, x_eval_, kind="linear", assume_sorted=True
)
xx_kde[i] = x_kde_interp(q)
return xx_kde
[docs]
class RainMix(_KDE, _Rain):
supplements_names = (
"q_thresh",
"f_thresh",
"kernel_data",
"q_kde_eval",
"x_eval",
)
def __init__(
self,
distribution,
threshold=0.0015,
q_thresh_lower=0.6,
q_thresh_upper=0.95,
):
"""Mixed Rain distribution with lower threshold and KDE above
q_threshold.
Requires self-generated sample data to be initialized because
of the KDE.
Parameters
----------
distribution : Dist
threshold : float, optional
Threshold rain intensity
q_threshold : float, optional
Quantile above which to estimate the distribution via KDE.
Not the quantile of the full distribution, but of the
'wet' values.
"""
self.dist = distribution
self.thresh = threshold
self.q_thresh_lower = q_thresh_lower
self.q_thresh_upper = q_thresh_upper
self.debug = False
# this enables testing
sample_data = self._gen_sample_data()
self.sample_data = sample_data
self._feasible_start = self.fit(sample_data)
# parameter_names exists because of the DistMeta class
self.parameter_names = (
self.parameter_names + distribution.parameter_names
)
# n_pars is set in class construction time with the Dist metaclass
self.n_pars += self.dist.n_pars
self.name = "rainmix " + distribution.name
def _gen_sample_data(self):
# the following solves problems that arise during testing
# as we have a mixed parametric, non-parametric distribution
# here, supplying a feasible starting solution is non-trivial.
# the data is part of the solution!
sample_quantiles = np.linspace(0.001, 0.999, 1000)
sample_data = self.dist.ppf(
sample_quantiles, *self.dist._feasible_start
)
sample_data[-50:] *= 2
sample_data = np.concatenate((np.zeros_like(sample_data), sample_data))
return sample_data
def _pdf(
self,
x,
q_thresh,
rain_prob,
kernel_width,
kernel_data,
f_thresh,
q_kde_eval=None,
x_eval=None,
*args,
**kwds,
):
(x, q_thresh, rain_prob, kernel_width, f_thresh) = np.atleast_1d(
x, q_thresh, rain_prob, kernel_width, f_thresh
)
assert np.all(np.isfinite(rain_prob))
if len(rain_prob) == 1:
rain_prob = np.broadcast_to(rain_prob, x.shape)
if len(kernel_width) == 1:
kernel_width = np.broadcast_to(kernel_width, x.shape)
if isinstance(kernel_data, np.ndarray):
if kernel_data.ndim == 1 and kernel_data.dtype != np.dtype("O"):
kernel_data = np.broadcast_to(
kernel_data, (x.size, kernel_data.size)
)
else:
kernel_data = np.asarray(kernel_data, dtype=object)
if len(q_thresh) == 1:
q_thresh = np.broadcast_to(q_thresh, x.shape)
if len(f_thresh) == 1:
f_thresh = np.broadcast_to(f_thresh, x.shape)
if x.ndim > 1 and rain_prob.ndim > 1:
(
x,
rain_prob,
f_thresh,
) = np.broadcast_arrays(x, rain_prob, f_thresh)
rain_mask = x > self.thresh
par_mask = (x < f_thresh) & rain_mask
kde_mask = x >= f_thresh
# the following tries to catch broadcast ambiguities, but might not
# be the most parsimonious formulation in terms of memory
dens = np.zeros_like(rain_prob + rain_mask)
def par_scale(x):
return (x - self.thresh) / (f_thresh[par_mask] - self.thresh)
def par_rescale(dens):
return dens / (f_thresh[par_mask] - self.thresh) * self.q_thresh
# parametric part
if np.any(par_mask):
x_par = x[par_mask]
x_par = par_scale(x_par)
kwds_par = self._mask_kwds(par_mask, kwds)
dens_par = self.dist.pdf(x_par, *args, **kwds_par)
dens_par = par_rescale(dens_par)
if dens.shape == par_mask.shape:
try:
dens[par_mask] = np.atleast_1d(dens_par)
except TypeError:
dens[par_mask] = np.atleast_1d(dens_par)[0]
else:
dens[:, par_mask[0]] = dens_par
def kde_single(dens, x_kde, mask, kernel_width, kernel_data):
if len(kernel_data.ravel()) == 0:
dens[mask] = 0
return
dens_kde = kde.kernel_density(
kernel_width, np.log(kernel_data), eval_points=np.log(x_kde)
)
# dividing by x_kde because of chain-rule differentiation
# of ln(x)
dens[mask] = (1.0 - self.q_thresh) * dens_kde / x_kde
if kernel_data.size:
dens[mask] /= kernel_data.size - 1
# non-parametric (KDE) part
if np.any(kde_mask):
if kde_mask.ndim == 2:
for tt, mask in enumerate(kde_mask):
kde_single(
dens[tt],
x[tt, mask],
mask,
kernel_width[tt],
kernel_data[tt],
)
else:
kde_single(
dens,
x[kde_mask],
kde_mask,
kernel_width[kde_mask],
kernel_data,
)
return dens
def _cdf(
self,
x,
q_thresh,
rain_prob,
kernel_width,
kernel_data,
f_thresh,
q_kde_eval=None,
x_eval=None,
*args,
**kwds,
):
(x, q_thresh, rain_prob, kernel_width, f_thresh) = np.atleast_1d(
x, q_thresh, rain_prob, kernel_width, f_thresh
)
# i want to manipulate kwds, but only here
kwds = {k: v for k, v in kwds.items()}
if len(rain_prob) == 1:
rain_prob = np.broadcast_to(rain_prob, x.shape)
if len(kernel_width) == 1:
kernel_width = np.broadcast_to(kernel_width, x.shape)
if isinstance(kernel_data, np.ndarray):
if kernel_data.ndim == 1 and kernel_data.dtype != np.dtype("O"):
kernel_data = np.broadcast_to(
kernel_data, (len(x), len(kernel_data))
)
else:
kernel_data = np.asarray(kernel_data, dtype=object)
if len(q_thresh) == 1:
q_thresh = np.broadcast_to(q_thresh, x.shape)
if len(f_thresh) == 1:
f_thresh = np.broadcast_to(f_thresh, x.shape)
upper_eval = np.atleast_1d(kwds.pop("u", [None]))
if len(upper_eval) == 1:
upper_eval = np.broadcast_to(upper_eval, x.shape)
if isinstance(q_kde_eval, list):
q_kde_eval = np.asarray(q_kde_eval, dtype=object)
if isinstance(x_eval, list):
x_eval = np.asarray(x_eval, dtype=object)
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)
rain_mask = x > self.thresh
kde_mask = x >= f_thresh
par_mask = ~kde_mask & rain_mask
# fill with lower part of south-of-p0 uniform random numbers
qq = np.zeros_like(x, dtype=float)
n_non_rain = (~rain_mask).sum()
if n_non_rain > 0:
qq[~rain_mask] = (1 - rain_prob[~rain_mask]) * varwg.rng.random(
n_non_rain
)
if self.debug:
fig, ax = plt.subplots(nrows=1, ncols=1)
ax.hist(qq[~rain_mask], 40, histtype="step", label="norain")
# parametric part
if np.any(par_mask):
x_par = x[par_mask]
kwds_par = self._mask_kwds(par_mask, kwds)
x_par = (x_par - self.thresh) / (f_thresh[par_mask] - self.thresh)
if "l" in kwds_par:
x_par = np.maximum(x_par, kwds_par["l"])
if self.debug:
assert np.all(x_par >= 0)
assert np.all(x_par <= 1)
if "l" not in kwds_par:
kwds_par["l"] = 0.0
if "u" not in kwds_par:
kwds_par["u"] = 1.0
q_par = self.dist.cdf(x_par, *args, **kwds_par)
p0 = (1 - rain_prob)[par_mask]
p1 = q_thresh[par_mask]
q_par = p0 + (p1 - p0) * q_par
qq[par_mask] = q_par
assert np.all(np.isfinite(q_par))
# non-parametric (KDE) part
if np.any(kde_mask):
assert np.all(kernel_width[kde_mask] > 0)
if isinstance(q_kde_eval, list):
q_kde_eval = [q_kde_eval[i] for i in np.where(kde_mask)[0]]
elif q_kde_eval.ndim == 2:
q_kde_eval = q_kde_eval[kde_mask]
if isinstance(x_eval, list):
x_eval = [x_eval[i] for i in np.where(kde_mask)[0]]
elif x_eval.ndim == 2:
x_eval = x_eval[kde_mask]
q_kde = self._kde_cdf(x[kde_mask], q_kde_eval, x_eval)
qq[kde_mask] = q_kde
if self.debug:
assert np.all(np.isfinite(q_kde))
if self.debug:
ax.hist(q_kde, 20, histtype="step", label="kde")
ax.legend()
assert np.all(np.isfinite(q_kde))
return qq
def _ppf(
self,
qq,
q_thresh,
rain_prob,
kernel_width,
kernel_data,
f_thresh,
q_kde_eval=None,
x_eval=None,
*args,
**kwds,
):
(qq, q_thresh, rain_prob, kernel_width, f_thresh) = np.atleast_1d(
qq, q_thresh, rain_prob, kernel_width, f_thresh
)
assert np.all(np.isfinite(qq))
# i want to manipulate kwds, but only here
kwds = {k: v for k, v in kwds.items()}
if len(rain_prob) == 1:
rain_prob = np.broadcast_to(rain_prob, qq.shape)
if len(kernel_width) == 1:
kernel_width = np.broadcast_to(kernel_width, qq.shape)
if isinstance(kernel_data, np.ndarray):
if kernel_data.ndim == 1 and kernel_data.dtype != np.dtype("O"):
kernel_data = np.broadcast_to(
kernel_data, (len(qq), len(kernel_data))
)
else:
kernel_data = np.asarray(kernel_data, dtype=object)
if len(q_thresh) == 1:
q_thresh = np.broadcast_to(q_thresh, qq.shape)
if len(f_thresh) == 1:
f_thresh = np.broadcast_to(f_thresh, qq.shape)
upper_eval = np.atleast_1d(kwds.pop("u", [None]))
if len(upper_eval) == 1:
upper_eval = np.broadcast_to(upper_eval, qq.shape)
if isinstance(q_kde_eval, list):
q_kde_eval = np.asarray(q_kde_eval, dtype=object)
if isinstance(x_eval, list):
x_eval = np.asarray(x_eval, dtype=object)
# we want to have zero rain where probability of rain is lower
# than rain_prob
x = np.zeros_like(qq, dtype=float)
rain_mask = qq > 1 - rain_prob
kde_mask = qq >= q_thresh
par_mask = rain_mask & ~kde_mask
if np.all(~rain_mask):
return x
# parametric part
if np.any(par_mask):
kwds_par = self._mask_kwds(par_mask, kwds)
q_par = qq[par_mask]
p0 = 1 - rain_prob[par_mask]
p1 = q_thresh[par_mask]
q_par = (q_par - p0) / (p1 - p0)
x_par = self.dist.ppf(q_par, *args, **kwds_par)
if self.debug:
assert np.all(q_par >= 0)
assert np.all(q_par <= 1)
x_par = (f_thresh[par_mask] - self.thresh) * x_par + self.thresh
x[par_mask] = x_par
assert np.all(np.isfinite(x_par))
# non-parametric (KDE) part
if np.any(kde_mask):
if (q_kde_eval.ndim == 2) or (
q_kde_eval.dtype == np.dtype("object")
):
q_kde_eval = q_kde_eval[kde_mask]
if (x_eval.ndim == 2) or (x_eval.dtype == np.dtype("object")):
x_eval = x_eval[kde_mask]
x_kde = self._kde_ppf(qq[kde_mask], q_kde_eval, x_eval)
x[kde_mask] = x_kde
assert np.all(np.isfinite(x_kde))
# # make sure there is no rain when rain_prob is very close to 0
# biblical_rain_mask = x[rain_prob < 1e-2] > 0
# if np.sum(biblical_rain_mask):
# # print(f"Capping biblical rain: {np.sum(biblical_rain_mask)} steps"
# # f"({100 * np.mean(biblical_rain_mask):.2f}%)")
# x[rain_prob < 1e-6] = 0
return x
[docs]
def mean(self, *args, **kwds):
f_thresh = kwds["f_thresh"]
epsilon = 1e-6
# integrate parametric and kde part separately
result_par = quad(
lambda x: x * self.pdf(x, *args, **kwds),
self.thresh,
f_thresh,
)
x_max = self.ppf(1 - epsilon, *args, **kwds)
result_kde = quad(
lambda x: x * self.pdf(x, *args, **kwds),
f_thresh,
x_max,
)
return 0.5 * (result_par[0] + result_kde[0])
def _fit_qthresh(
self, q_thresh, x, x0=None, kernel_bounds=None, *args, **kwds
):
rain_mask = x > self.thresh
rain_prob = np.mean(rain_mask)
self.q_thresh = q_thresh
f_thresh = max(np.percentile(x, 100 * self.q_thresh), self.thresh)
kde_mask = x >= f_thresh
par_mask = rain_mask & ~kde_mask
if np.sum(par_mask) > 1:
rain_kwds = self._mask_kwds(par_mask, kwds)
x_par = x[rain_mask & ~kde_mask]
x_par = (x_par - self.thresh) / (f_thresh - self.thresh)
par_names = self.dist.parameter_names
if "u" in par_names:
# the supplied upper limit is not meant for the
# parametric but for the kde part
rain_kwds["u"] = 1.0
# lower bound is not to be fitted
rain_kwds["l"] = 0.0
par_solution = list(self.dist.fit(x_par, **rain_kwds))
else:
par_solution = self.dist._feasible_start
self.success = True # haha!
kernel_data = x[kde_mask]
if np.any(par_mask):
kernel_data = np.concatenate(([np.max(x[par_mask])], kernel_data))
# if len(kernel_data) >= 5:
if len(kernel_data) > 2:
kernel_width = self.fit_kde(
np.log(kernel_data),
None if x0 is None else x0[2],
bounds=kernel_bounds,
)
if np.isnan(kernel_width):
kernel_width = None
else:
kernel_data = np.array([])
kernel_width = None
upper_eval = kwds.pop("u", None)
q_kde_eval, x_eval_log = self._kde_integral(
kernel_width, kernel_data, f_thresh, upper_eval
)
x_eval = np.exp(x_eval_log)
if kernel_width:
kernel_width = float(kernel_width)
self.fitted_pars = (
q_thresh,
rain_prob,
kernel_width,
kernel_data,
f_thresh,
q_kde_eval,
x_eval,
) + tuple(par_solution)
self.q_thresh = q_thresh
density = self.pdf(x, *self.fitted_pars)
mask = (density > 0) & np.isfinite(density)
lik = -np.nansum(np.log(density[mask]))
return lik
def _fit(self, x, x0=None, kernel_bounds=None, *args, **kwds):
result = sp_optimize.minimize_scalar(
partial(self._fit_qthresh, kernel_bounds=kernel_bounds, **kwds),
args=(x, x0) + args,
bounds=[self.q_thresh_lower, self.q_thresh_upper],
method="bounded",
)
return self.fitted_pars
[docs]
def fit_ml(self, x, x0=None, bounds=None, *args, **kwds):
"""This is a lie, we don't fit rainmix with ML."""
params = self._fit(x, x0, bounds=bounds, *args, **kwds)
params = list(params)
supplements = {
name: params[self.parameter_names.index(name)]
for name in self.supplements_names
}
fixed_names = [
name for name in kwds.keys() if name in self.parameter_names
]
params = [
param
for name, param in zip(self.parameter_names, params)
if (name not in self.supplements_names)
and (name not in fixed_names)
]
return self.Result(
x=params, supplements=supplements, success=self.success
)
def _constraints(
self,
x,
q_thresh,
rain_prob,
kernel_width,
kernel_data,
f_thresh,
q_kde_eval=None,
x_eval=None,
**kwds,
):
# (rain_prob, kernel_width, kernel_data, f_thresh) = map(
# np.asarray, (rain_prob, kernel_width, kernel_data, f_thresh)
# )
(rain_prob, kernel_width, f_thresh) = map(
np.asarray, (rain_prob, kernel_width, f_thresh)
)
if "kernel_bounds" in kwds:
kwds.pop("kernel_bounds")
mask = (
(q_thresh < self.q_thresh_lower)
| (q_thresh > self.q_thresh_upper)
| (rain_prob < 0)
| (rain_prob > 1)
| (kernel_width <= 0)
| (f_thresh < 0)
)
kde_mask = x >= f_thresh
par_mask = (x > self.thresh) & ~kde_mask
mask[par_mask] |= self.dist._constraints(
x[par_mask],
**{key: value[par_mask] for key, value in kwds.items()},
)
return mask
def _invalid_x(self, x, *args, **kwds):
return np.atleast_1d((np.full_like(x, False, dtype=bool)))
def _fix_x(self, x):
return np.where(x < 0, 0, x)
rainmix_kumaraswamy = RainMix(
kumaraswamy, threshold=0.001, q_thresh_lower=0.95, q_thresh_upper=0.99
)
# rainmix_expon = RainMix(expon, q_threshold=q_threshold)
# rainmix_weibull = RainMix(weibull, q_threshold=q_threshold)
# rainmix_gamma = RainMix(gamma, q_threshold=q_threshold)
# rainmix_lognormal = RainMix(lognormal, q_threshold=q_threshold)
# rainmix_gamma1 = RainMix(gamma1)
[docs]
class Rain(_Rain):
def __init__(self, distribution, threshold=0.0001):
self.dist = distribution
self.thresh = threshold
self.parameter_names = tuple(
["rain_prob"] + list(distribution.parameter_names)
)
self.n_pars = len(self.parameter_names)
self._feasible_start = (0.98,) + self.dist._feasible_start
self.name = "rain " + distribution.name
def _pdf(self, x, rain_prob, *args, **kwds):
x, rain_prob = np.atleast_1d(x, rain_prob)
if len(rain_prob) == 1:
rain_prob = np.broadcast_to(rain_prob, x.shape)
if x.ndim > 1 and rain_prob.ndim > 1:
x, rain_prob = np.broadcast_arrays(x, rain_prob)
# the following tries to catch broadcast ambiguities, but might not
# be the most parsimonious formulation in terms of memory
rain_mask = x > self.thresh
dens = np.empty_like(rain_prob + rain_mask, dtype=float)
dens[~rain_mask] = (1 - rain_prob[~rain_mask]) / self.thresh
kwds_rain = self._mask_kwds(rain_mask, kwds)
dens_rain = (
self.dist.pdf(x[rain_mask] - self.thresh, *args, **kwds_rain)
* rain_prob[rain_mask]
)
if dens.shape == rain_mask.shape:
try:
dens[rain_mask] = np.atleast_1d(dens_rain)
except TypeError:
dens[rain_mask] = np.atleast_1d(dens_rain)[0]
else:
dens[:, rain_mask[0]] = dens_rain
return dens
def _cdf(self, x, rain_prob, *args, **kwds):
x, rain_prob = np.atleast_1d(x, rain_prob)
qq = np.zeros_like(x, dtype=float)
if len(rain_prob) == 1:
rain_prob = np.full_like(x, rain_prob)
finite_mask = np.isfinite(x)
x_finite = np.where(finite_mask, x, 0)
rain_mask = x_finite > self.thresh
non_rain_mask = ~rain_mask
# subtracting self.thresh from x has the unwanted side-effect of
# shifting x to possibly out-of-bound values. so this should not be
# done with values where non_rain_mask.
# as we mask x, we also have to mask the parameters given in **kwds
kwds_rain = self._mask_kwds(rain_mask, kwds)
p0 = 1 - rain_prob[rain_mask]
qq[rain_mask] = p0 + rain_prob[rain_mask] * self.dist.cdf(
x[rain_mask] - self.thresh, *args, **kwds_rain
)
qq[non_rain_mask] = varwg.rng.uniform(size=non_rain_mask.sum()) * (
1 - rain_prob[non_rain_mask]
)
qq[~finite_mask] = np.nan
return qq
def _ppf(self, qq, rain_prob, *args, **kwds):
qq, rain_prob = np.atleast_1d(qq, rain_prob)
if len(rain_prob) == 1:
rain_prob = np.full_like(qq, rain_prob)
# we want to have zero rain, where probability of rain is lower than
# rain_prob
x = np.zeros_like(qq, dtype=float)
rain_mask = qq > 1 - rain_prob
if np.all(~rain_mask):
return x
# subtracting self.thresh from x has the unwanted side-effect of
# shifting x to possibly out-of-bound values. so this should not be
# done with values where ~rain_mask.
# as we mask x, we also have to mask the parameters given in **kwds
kwds_rain = self._mask_kwds(rain_mask, kwds)
x[rain_mask] = self.thresh + self.dist.ppf(
1 - (1 - qq[rain_mask]) / rain_prob[rain_mask], *args, **kwds_rain
)
return x
def _fit(self, x, **kwds):
rain_mask = x > self.thresh
rain_prob = np.mean(rain_mask)
if rain_prob:
solution = self.dist.fit(x[rain_mask] - self.thresh, **kwds)
else:
solution = np.full(self.dist.n_pars - len(kwds), np.nan)
return tuple([rain_prob] + list(solution))
rain_expon = Rain(expon, threshold=0.002)
rain_weibull = Rain(weibull, threshold=0.002)
if __name__ == "__main__":
# import config_konstanz_disag as conf
import config_konstanz as conf
import varwg
from varwg.core import base, plotting
varwg.set_conf(conf)
# times_hourly, met = varwg.read_met(varwg.conf.met_file)
# rain, times = vg.my.sumup(met["R"], 24, times_hourly)
met_vg = varwg.VG(
("R", "theta", "ILWR", "Qsw", "rh"), verbose=True, refit="R"
)
rain_dist, solution = met_vg.dist_sol["R"]
met_vg.simulate()
# rain = met_vg.data_raw[0]
# rain_dist = rainmix_expon
# rain_dist = vg.sd.SlidingDist(rain_dist, rain, met_vg.times, verbose=True)
# solution = rain_dist.fit()
# rain_dist.plot_seasonality_fit()
# rain_dist.plot_monthly_params()
rain_dist.plot_fourier_fit()
rain_dist.plot_monthly_fit(
solution,
# dists_alt=(Rain(gamma, threshold=.015),
# expon)
)
met_vg.plot_monthly_hists("R")
met_vg.plot_meteogramm_daily()
# rain_dist.plot_monthly_params()
quantiles = rain_dist.cdf(solution)
my.hist(quantiles, 20)
my.hist(met_vg.data_trans[0], 20, dist=norm)
# rain_dist.scatter_pdf(solution)
plt.show()