From 50473e412eed0893b29004b447b6c7aba10886d8 Mon Sep 17 00:00:00 2001 From: Nick Date: Thu, 26 Sep 2024 09:44:36 +0100 Subject: [PATCH] Updated code style to black --- .github/workflows/testing.yml | 2 +- pba/cbox.py | 72 +- pba/copula.py | 303 ++++-- pba/core.py | 175 ++-- pba/interval.py | 818 ++++++++-------- pba/logical.py | 226 ++--- pba/pbox.py | 1154 +++++++++++++---------- pba/pbox_constructors/distributions.py | 267 +++--- pba/pbox_constructors/non_parametric.py | 908 +++++++++--------- pba/pbox_constructors/parametric.py | 50 +- pyproject.toml | 2 +- 11 files changed, 2167 insertions(+), 1810 deletions(-) diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index eb9da83..c049ceb 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -1,4 +1,4 @@ -name: CI +name: Testing on: [push, pull_request] diff --git a/pba/cbox.py b/pba/cbox.py index f2ed730..342a0b1 100644 --- a/pba/cbox.py +++ b/pba/cbox.py @@ -2,51 +2,52 @@ from .pbox import Pbox else: from pbox import Pbox - + import numpy as np from matplotlib import pyplot as plt from typing import List, Tuple, Union -__all__ = ['Cbox'] +__all__ = ["Cbox"] + class Cbox(Pbox): """ Confidence boxes (c-boxes) are imprecise generalisations of traditional confidence distributions - They have a different interpretation to p-boxes but rely on the same underlying mathematics. As such in pba-for-python c-boxes inhert most of their methods from Pbox. + They have a different interpretation to p-boxes but rely on the same underlying mathematics. As such in pba-for-python c-boxes inhert most of their methods from Pbox. Args: Pbox (_type_): _description_ """ - - def __init__(self,*args,**kwargs): - if len(args) == 1 and isinstance(args[0],Pbox): - - super().__init__(**vars(args[0])) - + + def __init__(self, *args, **kwargs): + if len(args) == 1 and isinstance(args[0], Pbox): + + super().__init__(**vars(args[0])) + else: - - super().__init__(*args,**kwargs) - + + super().__init__(*args, **kwargs) + def __repr__(self): if self.mean_left == self.mean_right: - mean_text = f'{round(self.mean_left, 4)}' + mean_text = f"{round(self.mean_left, 4)}" else: - mean_text = f'[{round(self.mean_left, 4)}, {round(self.mean_right, 4)}]' + mean_text = f"[{round(self.mean_left, 4)}, {round(self.mean_right, 4)}]" if self.var_left == self.var_right: - var_text = f'{round(self.var_left, 4)}' + var_text = f"{round(self.var_left, 4)}" else: - var_text = f'[{round(self.var_left, 4)}, {round(self.var_right, 4)}]' + var_text = f"[{round(self.var_left, 4)}, {round(self.var_right, 4)}]" - range_text = f'[{round(np.min([self.left, self.right]), 4), round(np.max([self.left, self.right]), 4)}' + range_text = f"[{round(np.min([self.left, self.right]), 4), round(np.max([self.left, self.right]), 4)}" if self.shape is None: - shape_text = ' ' + shape_text = " " else: - shape_text = f' {self.shape}' # space to start; see below lacking space + shape_text = f" {self.shape}" # space to start; see below lacking space - return f'Cbox: ~ {shape_text} (range={range_text}, mean={mean_text}, var={var_text})' + return f"Cbox: ~ {shape_text} (range={range_text}, mean={mean_text}, var={var_text})" __str__ = __repr__ @@ -66,8 +67,8 @@ def __radd__(self, other): else: return Cbox(super().__radd__(other)) - def add(self,other, method = 'f'): - return Cbox(super().add(other, method = method)) + def add(self, other, method="f"): + return Cbox(super().add(other, method=method)) def __sub__(self, other): if isinstance(other, Cbox): @@ -85,8 +86,8 @@ def __rsub__(self, other): else: return Cbox(super().__rsub__(other)) - def sub(self,other, method = 'f'): - return Cbox(super().sub(other, method = method)) + def sub(self, other, method="f"): + return Cbox(super().sub(other, method=method)) def __mul__(self, other): if isinstance(other, Cbox): @@ -104,8 +105,8 @@ def __rmul__(self, other): else: return Cbox(super().__rmul__(other)) - def mul(self,other, method = 'f'): - return Cbox(super().mul(other, method = method)) + def mul(self, other, method="f"): + return Cbox(super().mul(other, method=method)) def __truediv__(self, other): if isinstance(other, Cbox): @@ -123,16 +124,21 @@ def __rtruediv__(self, other): else: return Cbox(super().__rtruediv__(other)) - def div(self,other, method = 'f'): - return Cbox(super().div(other, method = method)) + def div(self, other, method="f"): + return Cbox(super().div(other, method=method)) def __neg__(self): return Cbox(super().__neg__()) - + def recip(self): return Cbox(super().recip()) - -def singh(cboxes: List[Cbox], theta: Union[float, List[float]], figax: Tuple[plt.Figure, plt.Axes] = None) -> Tuple[plt.Figure, plt.Axes]: + + +def singh( + cboxes: List[Cbox], + theta: Union[float, List[float]], + figax: Tuple[plt.Figure, plt.Axes] = None, +) -> Tuple[plt.Figure, plt.Axes]: """ Generates a singh plot. Create s @@ -143,7 +149,7 @@ def singh(cboxes: List[Cbox], theta: Union[float, List[float]], figax: Tuple[plt Returns: tuple: A tuple containing the Matplotlib figure and axes objects used for plotting. - + """ x = np.linspace(0, 1, 1001) @@ -164,6 +170,6 @@ def singh(cboxes: List[Cbox], theta: Union[float, List[float]], figax: Tuple[plt ax.plot([0] + list(x) + [1], [0] + left + [1]) ax.plot([0] + list(x) + [1], [0] + right + [1]) - ax.plot([0, 1], [0, 1], 'k--', lw=2) + ax.plot([0, 1], [0, 1], "k--", lw=2) return fig, ax diff --git a/pba/copula.py b/pba/copula.py index 490d78f..636e5b5 100644 --- a/pba/copula.py +++ b/pba/copula.py @@ -12,7 +12,7 @@ from .interval import * else: from interval import * - + import numpy as np from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import Axes3D @@ -20,14 +20,34 @@ from scipy.stats import multivariate_normal as mvn from scipy.stats import norm -__all__ = ['Copula','ClaGen','ClaInv','FGen','FInv','indep','perf','opp','Cla','F','Gau','pi','M','W','Frank','Clayton','Gaussian'] +__all__ = [ + "Copula", + "ClaGen", + "ClaInv", + "FGen", + "FInv", + "indep", + "perf", + "opp", + "Cla", + "F", + "Gau", + "pi", + "M", + "W", + "Frank", + "Clayton", + "Gaussian", +] + + class Copula(object): - def __init__(self, cdf=None, func=None, param = None): + def __init__(self, cdf=None, func=None, param=None): - self.cdf = cdf # Could also include density, however expensive to compute - self.func = func # If the functional form is known - self.param = param # parameter for func + self.cdf = cdf # Could also include density, however expensive to compute + self.func = func # If the functional form is known + self.param = param # parameter for func def __repr__(self): @@ -36,42 +56,54 @@ def __repr__(self): if self.func is not None: func = self.func - if func == indep: func = "π" - if func == perf: func = "M" - if func == opp: func = "W" - if func == Gau: func = "Gau" - if func == Cla: func = "Clayton" - if func == F: func = "Frank" - statement1 = f'{func}' + if func == indep: + func = "π" + if func == perf: + func = "M" + if func == opp: + func = "W" + if func == Gau: + func = "Gau" + if func == Cla: + func = "Clayton" + if func == F: + func = "Frank" + statement1 = f"{func}" if self.param is not None: func = self.func parName = "par" - if func == Gau: parName = "r" - if func == F: parName = "s" - if func == Cla: parName = "t" - statement2 = f'{parName}={self.param}' + if func == Gau: + parName = "r" + if func == F: + parName = "s" + if func == Cla: + parName = "t" + statement2 = f"{parName}={self.param}" - return f'Copula ~ {statement1}({statement2})' + return f"Copula ~ {statement1}({statement2})" - def get_cdf(self, x, y): # x, y are points on the unit square + def get_cdf(self, x, y): # x, y are points on the unit square - if self.func is not None: # If function is known, return it - if self.param is not None: return self.func(x, y, self.param) - return self.func(x,y) + if self.func is not None: # If function is known, return it + if self.param is not None: + return self.func(x, y, self.param) + return self.func(x, y) - else: # Simple inner/outter interpolation. Returns interval for cdf value + else: # Simple inner/outter interpolation. Returns interval for cdf value xsize, ysize = self.cdf.shape - xIndexLower = int(np.floor(x * (xsize-1))) - yIndexLower = int(np.floor(y * (ysize-1))) + xIndexLower = int(np.floor(x * (xsize - 1))) + yIndexLower = int(np.floor(y * (ysize - 1))) - xIndexUpper = int(np.ceil(x * (xsize-1))) - yIndexUpper = int(np.ceil(y * (ysize-1))) + xIndexUpper = int(np.ceil(x * (xsize - 1))) + yIndexUpper = int(np.ceil(y * (ysize - 1))) - return Interval(self.cdf[xIndexLower, yIndexLower], self.cdf[xIndexUpper, yIndexUpper]) + return Interval( + self.cdf[xIndexLower, yIndexLower], self.cdf[xIndexUpper, yIndexUpper] + ) - def get_mass(self, x, y): # x, y are intervals on the unit square + def get_mass(self, x, y): # x, y are intervals on the unit square C22 = self.get_cdf(x.hi(), y.hi()) C21 = self.get_cdf(x.hi(), y.lo()) @@ -80,147 +112,209 @@ def get_mass(self, x, y): # x, y are intervals on the unit square return C22 - C21 - C12 + C11 - def show(self, pn = 50, fontsize = 20, cols = cm.RdGy): + def show(self, pn=50, fontsize=20, cols=cm.RdGy): ## # All the extra stuff is so that no more than 200 elements are plotted ## - A = self.cdf; m = len(A) - - if m < pn: pn = m - - x = y = np.linspace(0, 1, num = pn) - X, Y = np.meshgrid(x,y) - - nm = round(m/pn) - Z = A[::nm,::nm] # Skip over evelemts - - fig = plt.figure("SurfacPlots",figsize=(10,10)) - ax = fig.add_subplot(1,1,1,projection="3d") - ax.plot_surface(X, Y, Z, rstride=2,edgecolors="k", cstride=2, alpha=0.8, linewidth=0.25, cmap = cols) - plt.xlabel("X",fontsize = fontsize); plt.ylabel("Y", fontsize = fontsize) + A = self.cdf + m = len(A) + + if m < pn: + pn = m + + x = y = np.linspace(0, 1, num=pn) + X, Y = np.meshgrid(x, y) + + nm = round(m / pn) + Z = A[::nm, ::nm] # Skip over evelemts + + fig = plt.figure("SurfacPlots", figsize=(10, 10)) + ax = fig.add_subplot(1, 1, 1, projection="3d") + ax.plot_surface( + X, + Y, + Z, + rstride=2, + edgecolors="k", + cstride=2, + alpha=0.8, + linewidth=0.25, + cmap=cols, + ) + plt.xlabel("X", fontsize=fontsize) + plt.ylabel("Y", fontsize=fontsize) plt.show() - def showContour(self, fontsize = 20, cols = cm.coolwarm): + def showContour(self, fontsize=20, cols=cm.coolwarm): ## # All the extra stuff is so that no more than 200 elements are plotted ## - pn = 200 # Max plot number - A = self.cdf; m = len(A) - - if m < pn: pn = m - - x = y = np.linspace(0, 1, num = pn) - X, Y = np.meshgrid(x,y) - - nm = round(m/pn) - Z = A[::nm,::nm] # Skip over evelemts - - fig = plt.figure("SurfacPlots",figsize=(10,10)) - ax = fig.add_subplot(2,1,1,projection="3d") - ax.plot_surface(X, Y, Z, rstride=2,edgecolors="k", cstride=2, alpha=0.8, linewidth=0.25, cmap = cols) - plt.xlabel("X",fontsize = fontsize); plt.ylabel("Y", fontsize = fontsize) - plt.title("Surface Plot", fontsize = fontsize) - - ax1 = fig.add_subplot(2,1,2) - cp = ax1.contour(X, Y, Z, cmap = cols, levels = 15) + pn = 200 # Max plot number + A = self.cdf + m = len(A) + + if m < pn: + pn = m + + x = y = np.linspace(0, 1, num=pn) + X, Y = np.meshgrid(x, y) + + nm = round(m / pn) + Z = A[::nm, ::nm] # Skip over evelemts + + fig = plt.figure("SurfacPlots", figsize=(10, 10)) + ax = fig.add_subplot(2, 1, 1, projection="3d") + ax.plot_surface( + X, + Y, + Z, + rstride=2, + edgecolors="k", + cstride=2, + alpha=0.8, + linewidth=0.25, + cmap=cols, + ) + plt.xlabel("X", fontsize=fontsize) + plt.ylabel("Y", fontsize=fontsize) + plt.title("Surface Plot", fontsize=fontsize) + + ax1 = fig.add_subplot(2, 1, 2) + cp = ax1.contour(X, Y, Z, cmap=cols, levels=15) ax1.clabel(cp, inline=1, fontsize=10) - plt.xlabel("X", fontsize = fontsize); plt.ylabel("Y", fontsize = fontsize) - plt.title("Contour Plot", fontsize = fontsize) + plt.xlabel("X", fontsize=fontsize) + plt.ylabel("Y", fontsize=fontsize) + plt.title("Contour Plot", fontsize=fontsize) plt.tight_layout() plt.show() + ### # Copula generators for Archimedean (Frank and Clayton) copulas. Allows for easy and accurate copula generation # at any dimension in terms of univariate functions (generator and inverse generator). ### -def ClaGen(x, t = 1): return 1/t * (x **(-t) -1 ) # Clayton copula generator -def ClaInv(x, t = 1): return (1 + t * x) ** (-1/t) # Inverse generator -def FGen(x, s = 1): # Frank generator - X1 = np.exp( -(x * s) ) -1 - X2 = np.exp( -s ) -1 - return -np.log( X1 / X2) +def ClaGen(x, t=1): + return 1 / t * (x ** (-t) - 1) # Clayton copula generator + + +def ClaInv(x, t=1): + return (1 + t * x) ** (-1 / t) # Inverse generator + + +def FGen(x, s=1): # Frank generator + X1 = np.exp(-(x * s)) - 1 + X2 = np.exp(-s) - 1 + return -np.log(X1 / X2) -def FInv(x, s = 1): # Inverse - X1 = np.exp( -x ) * (np.exp(-s) - 1) + +def FInv(x, s=1): # Inverse + X1 = np.exp(-x) * (np.exp(-s) - 1) X2 = np.log(1 + X1) - return - (X2 / s) + return -(X2 / s) + ### # Copula functions and constructors ### -def indep(x,y): return x*y -def perf(x,y): return min(x,y) -def opp(x,y): return max(x+y-1,0) -def Cla(x, y, t = 1): return ClaInv( ClaGen(x, t) + ClaGen(y, t), t) -def F(x,y,s = 1): return FInv( FGen(x, s) + FGen(y, s), s) -def Gau(x,y,r=0): - if x == 0: return 0 - if y == 0: return 0 - return mvn.cdf([norm.ppf(x), norm.ppf(y)], mean = [0, 0], cov=[[1, r], [r, 1]]) +def indep(x, y): + return x * y + + +def perf(x, y): + return min(x, y) + + +def opp(x, y): + return max(x + y - 1, 0) + + +def Cla(x, y, t=1): + return ClaInv(ClaGen(x, t) + ClaGen(y, t), t) +def F(x, y, s=1): + return FInv(FGen(x, s) + FGen(y, s), s) + + +def Gau(x, y, r=0): + if x == 0: + return 0 + if y == 0: + return 0 + return mvn.cdf([norm.ppf(x), norm.ppf(y)], mean=[0, 0], cov=[[1, r], [r, 1]]) + # Copula constructors -def pi(steps = 200): +def pi(steps=200): x = y = np.linspace(0, 1, num=steps) cdf = np.array([[xs * ys for xs in x] for ys in y]) return Copula(cdf, indep) -def M(steps = 200): + +def M(steps=200): x = y = np.linspace(0, 1, num=steps) cdf = np.array([[perf(xs, ys) for xs in x] for ys in y]) return Copula(cdf, perf) -def W(steps = 200): + +def W(steps=200): x = y = np.linspace(0, 1, num=steps) cdf = np.array([[opp(xs, ys) for xs in x] for ys in y]) return Copula(cdf, opp) -def Frank(s = 0, steps = 200): # s is real; inf for perfect, 0 for indep, -inf for oposite - if s == float('-inf'): # Limit should be set earlier +def Frank( + s=0, steps=200 +): # s is real; inf for perfect, 0 for indep, -inf for oposite + + if s == float("-inf"): # Limit should be set earlier C = W() - return Copula(C.cdf, F, float('-inf')) + return Copula(C.cdf, F, float("-inf")) if s == 0: C = pi() return Copula(C.cdf, F, 1) - if s == float('inf'): + if s == float("inf"): C = M() - return Copula(C.cdf, F, float('inf')) + return Copula(C.cdf, F, float("inf")) x = y = np.linspace(0, 1, num=steps) cdf = np.array([[F(xs, ys, s) for xs in x] for ys in y]) return Copula(cdf, F, s) -def Clayton(t = 1, steps = 200): # t>-1; -1 for opposite, 0 for indep and inf for perfect + +def Clayton(t=1, steps=200): # t>-1; -1 for opposite, 0 for indep and inf for perfect if t == 0: C = pi() - return Copula(C.cdf, Cla , 0) + return Copula(C.cdf, Cla, 0) if t == -1: C = W() - return Copula(C.cdf, Cla ,-1) - if t == float('inf'): # Limit should be set earlier + return Copula(C.cdf, Cla, -1) + if t == float("inf"): # Limit should be set earlier C = M() - return Copula(C.cdf, Cla , float('inf')) + return Copula(C.cdf, Cla, float("inf")) x = np.linspace(0, 1, num=steps) - xx, yy = np.meshgrid(x,x,indexing='ij') + xx, yy = np.meshgrid(x, x, indexing="ij") X = np.array([xx.flatten(), yy.flatten()]) - cdf = Cla(X[0, ], X[1, ], t) + cdf = Cla(X[0,], X[1,], t) cdf = cdf.reshape(steps, steps) - cdf[0,] = 0; cdf[:,0] = 0 # Grounds C + cdf[0,] = 0 + cdf[:, 0] = 0 # Grounds C return Copula(cdf, Cla, t) -def Gaussian(r = 0, steps = 200): # -1 <= r <=1 ; -1 for opposite, 1 for indep and 1 for perfect + +def Gaussian( + r=0, steps=200 +): # -1 <= r <=1 ; -1 for opposite, 1 for indep and 1 for perfect if r == 0: C = pi() @@ -233,12 +327,13 @@ def Gaussian(r = 0, steps = 200): # -1 <= r <=1 ; -1 for opposite, 1 for ind return Copula(C.cdf, Gau, 1) x = np.linspace(0, 1, num=steps) - xx, yy = np.meshgrid(x,x,indexing='ij') + xx, yy = np.meshgrid(x, x, indexing="ij") X = np.array([xx.flatten(), yy.flatten()]) - cdf = mvn.cdf(norm.ppf(X.transpose()), mean=[0,0], cov=[[1, r], [r, 1]]) + cdf = mvn.cdf(norm.ppf(X.transpose()), mean=[0, 0], cov=[[1, r], [r, 1]]) cdf = cdf.reshape(steps, steps) - cdf[0,] = 0; cdf[:,0] = 0 # Grounds C + cdf[0,] = 0 + cdf[:, 0] = 0 # Grounds C return Copula(cdf, Gau, r) diff --git a/pba/core.py b/pba/core.py index 5c86339..c3ebc94 100644 --- a/pba/core.py +++ b/pba/core.py @@ -6,72 +6,77 @@ import numpy as np import warnings + def envelope(*args: Union[Interval, Pbox, float]) -> Union[Interval, Pbox]: - ''' + """ .. _core.envelope: - + Allows the envelope to be calculated for intervals and p-boxes. - + The envelope is the smallest interval/pbox that contains all values within the arguments. - + **Parameters**: ``*args``: The arguments for which the envelope needs to be calculated. The arguments can be intervals, p-boxes, or floats. - + **Returns**: ``Pbox|Interval``: The envelope of the given arguments, which can be an interval or a p-box. - + .. error:: - + ``ValueError``: If less than two arguments are given. - + ``TypeError``: If none of the arguments are intervals or p-boxes. - - ''' + + """ # args = list(*args) - + # Raise error if <2 arguments are given if len(args) < 2: - raise ValueError('At least two arguments are required') - + raise ValueError("At least two arguments are required") + # get the type of all arguments types = [arg.__class__.__name__ for arg in args] - + # check if all arguments are intervals or pboxes - if 'Interval' not in types and 'Pbox' not in types: - raise ValueError('At least one argument needs to be an Interval or Pbox') + if "Interval" not in types and "Pbox" not in types: + raise ValueError("At least one argument needs to be an Interval or Pbox") # check if there is a p-box in the arguments - elif 'Pbox' in types: + elif "Pbox" in types: # find first p-box - i = types.index('Pbox') + i = types.index("Pbox") # move previous values to the end args = args[i:] + args[:i] - + e = args[0].env(args[1]) for arg in args[2:]: - if not isinstance(arg,(Pbox,Interval)): - arg = Interval(arg,arg) + if not isinstance(arg, (Pbox, Interval)): + arg = Interval(arg, arg) e = e.env(arg) - - else: #Intervals only - - left = np.min([arg.left if isinstance(arg,Interval) else arg for arg in args]) - - right = np.max([arg.right if isinstance(arg,Interval) else arg for arg in args]) - - e = Interval(left,right) - + + else: # Intervals only + + left = np.min([arg.left if isinstance(arg, Interval) else arg for arg in args]) + + right = np.max( + [arg.right if isinstance(arg, Interval) else arg for arg in args] + ) + + e = Interval(left, right) + return e + def env(*args): - ''' + """ .. warning:: - + Deprecated function, use envelope() instead. - - ''' - warnings.warn('env() is deprecated, use envelope() instead', DeprecationWarning) + + """ + warnings.warn("env() is deprecated, use envelope() instead", DeprecationWarning) return envelope(*args) + # def min(x,y): # if x.__class__.__name__ == 'Pbox': # return x.min(y) @@ -88,101 +93,107 @@ def env(*args): # else: # raise NotImplementedError('At least one argument needs to be a Pbox') -def sum(*args: Union[list,tuple] ,method = 'f'): - ''' + +def sum(*args: Union[list, tuple], method="f"): + """ Allows the sum to be calculated for intervals and p-boxes - + **Parameters**: - + ``*args``: pboxes or intervals ``method`` (``f,i,o,p``): addition method to be used - + **Returns**: - + ``Interval | Pbox``: sum of interval or pbox objects within ``*args`` - + .. note:: - + If a list or tuple is given as the first argument, the elements of the list or tuple are used as arguments. If only one (non-list) argument is given, the argument is returned. - - - ''' - + + + """ + if len(args) == 1: - if isinstance(args[0],[list,tuple,np.ndarray]): + if isinstance(args[0], [list, tuple, np.ndarray]): args = args[0] else: return args[0] - + s = 0 for o in args: - if isinstance(o,[Interval,Pbox,Cbox]): - s = o.add(s,method = method) + if isinstance(o, [Interval, Pbox, Cbox]): + s = o.add(s, method=method) else: s += o return s -def mean(*args: Union[list,tuple] ,method = 'f'): - ''' + +def mean(*args: Union[list, tuple], method="f"): + """ Allows the mean to be calculated for intervals and p-boxes - + **Parameters**: l : list of pboxes or intervals - + method : pbox addition method to be used - + **Output**: ``Interval | Pbox``: mean of interval or pbox objects within ``*args`` - + .. important:: - + Implemented as - + >>> pba.sum(*args,method = method)/len(args) - ''' - s = sum(*args,method = method) - - return s/len(args) + """ + s = sum(*args, method=method) + + return s / len(args) + -def mul(*args, method = None): - for i,arg in enumerate(args): +def mul(*args, method=None): + for i, arg in enumerate(args): if i == 0: n = arg - elif n.__class__.__name__ == 'Interval': - if arg.__class__.__name__ == 'Interval': + elif n.__class__.__name__ == "Interval": + if arg.__class__.__name__ == "Interval": if method is None: n *= arg - elif method == 'p': + elif method == "p": n = n.pmul(arg) - elif method == 'o': + elif method == "o": n = n.omul(arg) else: - raise Exception(f"Method {method} unknown for Interval * Interval calculation") - elif arg.__class__.__name__ == 'Pbox': - n = arg.mul(n,method = method) + raise Exception( + f"Method {method} unknown for Interval * Interval calculation" + ) + elif arg.__class__.__name__ == "Pbox": + n = arg.mul(n, method=method) else: n *= arg - elif n.__class__.__name__ == 'Pbox': + elif n.__class__.__name__ == "Pbox": if method is None: n *= arg else: - n = n.mul(arg,method = method) + n = n.mul(arg, method=method) else: n *= arg return n - + + def sqrt(a): - if a.__class__.__name__ == 'Interval': - return Interval(np.sqrt(a.left),np.sqrt(a.right)) - elif a.__class__.__name__ == 'PBox': + if a.__class__.__name__ == "Interval": + return Interval(np.sqrt(a.left), np.sqrt(a.right)) + elif a.__class__.__name__ == "PBox": return Pbox( - left = np.sqrt(a.left), - right = np.sqrt(a.right), + left=np.sqrt(a.left), + right=np.sqrt(a.right), ) - else: + else: - return np.sqrt(a) + return np.sqrt(a) diff --git a/pba/interval.py b/pba/interval.py index 7b7870c..6ac5787 100644 --- a/pba/interval.py +++ b/pba/interval.py @@ -1,6 +1,7 @@ r""" """ + from typing import Union import numpy as np import random as r @@ -8,10 +9,11 @@ import warnings -__all__ = ['Interval','I','PM'] +__all__ = ["Interval", "I", "PM"] + class Interval: - r''' + r""" An interval is an uncertain number for which only the endpoints are known, :math:`x=[a,b]`. This is interpreted as :math:`x` being between :math:`a` and :math:`b` but with no more information about the value of :math:`x`. @@ -20,7 +22,7 @@ class Interval: Creation ________ - + Intervals can be created using either of the following: .. code-block:: python @@ -31,7 +33,7 @@ class Interval: Interval [2,3] .. tip:: - + The shorthand ``I`` is an alias for ``Interval`` Intervals can also be created from a single value ± half-width: @@ -65,24 +67,27 @@ class Interval: Alternative arithmetic methods are described in `interval.add`_, `interval.sub`_, `interval.mul`_, `interval.div`_. **Attributes**: - + ``left``: The left boundary of the interval. - + ``right``: The right boundary of the interval. - + .. admonition:: Default values - + If only 1 argument is given then the interval is assumed to be zero width around this value. - + If no arguments are given then the interval is assumed to be vaccous (i.e. :math:`[-\infty,\infty]`). This is implemented as ``Interval(-np.inf,np.inf)``. - ''' - def __init__(self,left = None, right = None): - + """ + + def __init__(self, left=None, right=None): + # disallow p-boxes - if left.__class__.__name__ == 'Pbox' or right.__class__.__name__ == 'Pbox': - raise ValueError("left and right must not be P-boxes. Use Pbox methods instead.") - + if left.__class__.__name__ == "Pbox" or right.__class__.__name__ == "Pbox": + raise ValueError( + "left and right must not be P-boxes. Use Pbox methods instead." + ) + # assume vaccous if no inputs if left is None and right is None: right = np.inf @@ -94,17 +99,16 @@ def __init__(self,left = None, right = None): elif left is not None and right is None: right = left - if hasattr(left, '__iter__') and not isinstance(left, Interval): - left = Interval(min(left),max(left)) - if hasattr(right, '__iter__') and not isinstance(right, Interval): - right = Interval(min(right),max(right)) - + if hasattr(left, "__iter__") and not isinstance(left, Interval): + left = Interval(min(left), max(left)) + if hasattr(right, "__iter__") and not isinstance(right, Interval): + right = Interval(min(right), max(right)) + if isinstance(left, Interval): left = left.left if isinstance(right, Interval): right = right.right - if left > right: LowerUpper = [left, right] left = min(LowerUpper) @@ -113,19 +117,22 @@ def __init__(self,left = None, right = None): self.left = left self.right = right - def __repr__(self) -> str: # return - return "Interval [%g, %g]"%(self.left,self.right) - - def __str__(self) -> str: # str(Interval) - return "[%g, %g]"%(self.left,self.right) + def __repr__(self) -> str: # return + return "Interval [%g, %g]" % (self.left, self.right) + def __str__(self) -> str: # str(Interval) + return "[%g, %g]" % (self.left, self.right) def __format__(self, format_spec: str) -> str: try: - return f'[{format(self.left, format_spec)},{format(self.right, format_spec)}]' + return ( + f"[{format(self.left, format_spec)},{format(self.right, format_spec)}]" + ) except: - raise ValueError(f'{format_spec} format specifier not understood for Interval object') - + raise ValueError( + f"{format_spec} format specifier not understood for Interval object" + ) + def __iter__(self): for bound in [self.left, self.right]: yield bound @@ -133,7 +140,7 @@ def __iter__(self): def __len__(self): return 2 - def __radd__(self,left): + def __radd__(self, left): return self.__add__(left) def __sub__(self, other): @@ -148,11 +155,11 @@ def __sub__(self, other): else: try: lo = self.left - other - hi = self.right - other + hi = self.right - other except: return NotImplemented - return Interval(lo,hi) + return Interval(lo, hi) def __rsub__(self, other): if other.__class__.__name__ == "Interval": @@ -171,12 +178,12 @@ def __rsub__(self, other): except: return NotImplemented - return Interval(lo,hi) + return Interval(lo, hi) def __neg__(self): return Interval(-self.right, -self.left) - - def __mul__(self,other): + + def __mul__(self, other): if other.__class__.__name__ == "Interval": b1 = self.lo() * other.lo() @@ -184,8 +191,8 @@ def __mul__(self,other): b3 = self.hi() * other.lo() b4 = self.hi() * other.hi() - lo = min(b1,b2,b3,b4) - hi = max(b1,b2,b3,b4) + lo = min(b1, b2, b3, b4) + hi = max(b1, b2, b3, b4) elif other.__class__.__name__ == "Pbox": @@ -202,12 +209,12 @@ def __mul__(self,other): return NotImplemented - return Interval(lo,hi) + return Interval(lo, hi) - def __rmul__(self,other): + def __rmul__(self, other): return self * other - def __truediv__(self,other): + def __truediv__(self, other): if other.__class__.__name__ == "Interval": @@ -227,147 +234,159 @@ def __truediv__(self,other): b3 = self.hi() / other.lo() b4 = self.hi() / other.hi() - lo = min(b1,b2,b3,b4) - hi = max(b1,b2,b3,b4) - + lo = min(b1, b2, b3, b4) + hi = max(b1, b2, b3, b4) + elif other.__class__.__name__ == "Pbox": return other.__rtruediv__(self) - + else: try: - lo = self.lo()/other - hi = self.hi()/other + lo = self.lo() / other + hi = self.hi() / other except: return NotImplemented - return Interval(lo,hi) + return Interval(lo, hi) + + def __rtruediv__(self, other): - def __rtruediv__(self,other): - if self.straddles_zero(): - + raise ZeroDivisionError() - + try: return other * self.recip() except: return NotImplemented - def __pow__(self,other): + def __pow__(self, other): if other.__class__.__name__ == "Interval": - pow1 = self.left ** other.left - pow2 = self.left ** other.right - pow3 = self.right ** other.left - pow4 = self.right ** other.right - powUp = max(pow1,pow2,pow3,pow4) - powLow = min(pow1,pow2,pow3,pow4) + pow1 = self.left**other.left + pow2 = self.left**other.right + pow3 = self.right**other.left + pow4 = self.right**other.right + powUp = max(pow1, pow2, pow3, pow4) + powLow = min(pow1, pow2, pow3, pow4) elif other.__class__.__name__ in ("int", "float"): - pow1 = self.left ** other - pow2 = self.right ** other - powUp = max(pow1,pow2) - powLow = min(pow1,pow2) + pow1 = self.left**other + pow2 = self.right**other + powUp = max(pow1, pow2) + powLow = min(pow1, pow2) if (self.right >= 0) and (self.left <= 0) and (other % 2 == 0): powLow = 0 - return Interval(powLow,powUp) + return Interval(powLow, powUp) - def __rpow__(self,left): + def __rpow__(self, left): if left.__class__.__name__ == "Interval": - pow1 = left.left ** self.left - pow2 = left.left ** self.right - pow3 = left.right ** self.left - pow4 = left.right ** self.right - powUp = max(pow1,pow2,pow3,pow4) - powLow = min(pow1,pow2,pow3,pow4) + pow1 = left.left**self.left + pow2 = left.left**self.right + pow3 = left.right**self.left + pow4 = left.right**self.right + powUp = max(pow1, pow2, pow3, pow4) + powLow = min(pow1, pow2, pow3, pow4) else: - pow1 = left ** self.left - pow2 = left ** self.right - powUp = max(pow1,pow2) - powLow = min(pow1,pow2) + pow1 = left**self.left + pow2 = left**self.right + powUp = max(pow1, pow2) + powLow = min(pow1, pow2) - return Interval(powLow,powUp) + return Interval(powLow, powUp) - def __lt__(self,other): + def __lt__(self, other): # < if not isinstance(other, Interval): try: other = Interval(other) except: - raise TypeError(f"'<' not supported between instances of Interval and {type(other)}") - + raise TypeError( + f"'<' not supported between instances of Interval and {type(other)}" + ) + if self.right < other.left: - return Interval(1,1).to_logical() + return Interval(1, 1).to_logical() elif self.left > other.right: - return Interval(0,0).to_logical() + return Interval(0, 0).to_logical() else: - return Interval(0,1).to_logical() + return Interval(0, 1).to_logical() def __rgt__(self, other): return self < other - - def __eq__(self,other): + + def __eq__(self, other): # == if not isinstance(other, Interval): try: other = Interval(other) except: - raise TypeError(f"'==' not supported between instances of Interval and {type(other)}") - + raise TypeError( + f"'==' not supported between instances of Interval and {type(other)}" + ) + if self.straddles(other.left) or self.straddles(other.right): - return Interval(0,1).to_logical() + return Interval(0, 1).to_logical() else: - return Interval(0,0).to_logical() + return Interval(0, 0).to_logical() - def __gt__(self,other): + def __gt__(self, other): # > try: lt = self < other except: - raise TypeError(f"'>' not supported between instances of Interval and {type(other)}") + raise TypeError( + f"'>' not supported between instances of Interval and {type(other)}" + ) return ~lt def __rlt__(self, other): return self > other - - def __ne__(self,other): + + def __ne__(self, other): # != try: eq = self == other except: - raise TypeError(f"'!=' not supported between instances of Interval and {type(other)}") + raise TypeError( + f"'!=' not supported between instances of Interval and {type(other)}" + ) return ~eq - - def __le__(self,other): + + def __le__(self, other): # <= if not isinstance(other, Interval): try: other = Interval(other) except: - raise TypeError(f"'<=' not supported between instances of Interval and {type(other)}") - + raise TypeError( + f"'<=' not supported between instances of Interval and {type(other)}" + ) + if self.right <= other.left: - return Interval(1,1).to_logical() + return Interval(1, 1).to_logical() elif self.left > other.right: - return Interval(0,0).to_logical() + return Interval(0, 0).to_logical() else: - return Interval(0,1).to_logical() - - def __ge__(self,other): + return Interval(0, 1).to_logical() + + def __ge__(self, other): if not isinstance(other, Interval): try: other = Interval(other) except: - raise TypeError(f"'>=' not supported between instances of Interval and {type(other)}") - + raise TypeError( + f"'>=' not supported between instances of Interval and {type(other)}" + ) + if self.right <= other.left: - return Interval(1,1).to_logical() + return Interval(1, 1).to_logical() elif self.left > other.right: - return Interval(0,0).to_logical() + return Interval(0, 0).to_logical() else: - return Interval(0,1).to_logical() + return Interval(0, 1).to_logical() def __bool__(self): @@ -377,323 +396,350 @@ def __bool__(self): else: return False except: - raise ValueError("Truth value of Interval %s is ambiguous" %self) - + raise ValueError("Truth value of Interval %s is ambiguous" % self) + def __abs__(self): if self.straddles_zero(): - return Interval(0, max(abs(self.left),abs(self.right))) + return Interval(0, max(abs(self.left), abs(self.right))) else: - return Interval(abs(self.left),abs(self.right)) + return Interval(abs(self.left), abs(self.right)) - def __contains__(self,other): - return self.straddles(other, endpoints = True) + def __contains__(self, other): + return self.straddles(other, endpoints=True) - def add(self,other,method = None): + def add(self, other, method=None): """ .. _interval.add: - + Adds the interval and another object together. - + **Args**: - + ``other``: The interval or numeric value to be added. This value must be transformable into an Interval object. - + **Methods**: - + p - perfect arithmetic :math:`[a,b]+[c,d] = [a + c, b + d]` - + o - opposite arithmetic :math:`[a,b]+[c,d] = [a + d, b + c]` - + None, i, f - Standard interval arithmetic is used. **Returns**: - + ``Interval`` """ if not isinstance(other, Interval): - if other.__class__.__name__ in ['Pbox','Cbox']: + if other.__class__.__name__ in ["Pbox", "Cbox"]: if method is None: - method = 'f' - return other.add(self, method = method) + method = "f" + return other.add(self, method=method) try: other = Interval(other) except: - raise TypeError(f"addition not supported between instances of Interval and {type(other)}") - - if method == 'p': + raise TypeError( + f"addition not supported between instances of Interval and {type(other)}" + ) + + if method == "p": Interval(self.left + other.left, self.right + other.right) - elif method == 'o': + elif method == "o": Interval(self.left + other.right, self.right + other.left) else: return self.__add__(other) - - def __add__(self,other): - if other.__class__.__name__ == 'Interval': + def __add__(self, other): + + if other.__class__.__name__ == "Interval": lo = self.left + other.left hi = self.right + other.right - elif other.__class__.__name__ == 'Pbox': + elif other.__class__.__name__ == "Pbox": # Perform Pbox addition assuming independance - return other.add(self, method = 'i') + return other.add(self, method="i") else: try: lo = self.left + other - hi = self.right + other + hi = self.right + other except: return NotImplemented - return Interval(lo,hi) - + return Interval(lo, hi) + def padd(self, other): """ .. warning:: This method is deprecated. Use add(other, method='p') instead. """ - warnings.warn("padd() is deprecated. Use add(other, method='p') instead.", DeprecationWarning) - return self.add(other, method='p') - - def oadd(self,other): + warnings.warn( + "padd() is deprecated. Use add(other, method='p') instead.", + DeprecationWarning, + ) + return self.add(other, method="p") + + def oadd(self, other): """ .. warning:: This method is deprecated. Use add(other, method='o') instead. """ - warnings.warn("oadd() is deprecated. Use add(other, method = 'o') instead.", DeprecationWarning) - return self.add(other, method='o') + warnings.warn( + "oadd() is deprecated. Use add(other, method = 'o') instead.", + DeprecationWarning, + ) + return self.add(other, method="o") - def sub(self, other, method = None): + def sub(self, other, method=None): """ .. _interval.sub: - + Subtracts other from self. - + **Args**: - - ``other``: The interval or numeric value to be subracted. This value must be transformable into an Interval object. - + + ``other``: The interval or numeric value to be subracted. This value must be transformable into an Interval object. + **Methods**: - + ``p``: perfect arithmetic :math:`a+b = [a.left - b.left, a.right - b.right]` - + ``o``: opposite arithmetic :math:`a+b = [a.left - b.right, a.right - b.left]` - + None, i, f - Standard interval arithmetic is used. **Returns**: - + ``Interval`` """ if not isinstance(other, Interval): - if other.__class__.__name__ in ['Pbox','Cbox']: + if other.__class__.__name__ in ["Pbox", "Cbox"]: if method is None: - method = 'f' - return other.rsub(self, method = method) + method = "f" + return other.rsub(self, method=method) try: other = Interval(other) except: - raise TypeError(f"Subtraction not supported between instances of Interval and {type(other)}") - - if method == 'p': + raise TypeError( + f"Subtraction not supported between instances of Interval and {type(other)}" + ) + + if method == "p": Interval(self.left - other.left, self.right - other.right) - elif method == 'o': + elif method == "o": Interval(self.left - other.right, self.right - other.left) else: return self.__sub__(other) - def psub(self,other): + def psub(self, other): """ .. warning:: Depreciated use self.sub(other, method = 'p') instead """ - warnings.warn("psub() is deprecated. Use sub(other, method = 'p') instead.", DeprecationWarning) + warnings.warn( + "psub() is deprecated. Use sub(other, method = 'p') instead.", + DeprecationWarning, + ) return Interval(self.left - other.left, self.right - other.right) - def osub(self,other): + def osub(self, other): """ .. warning:: Depreciated use self.sub(other, method = 'o') instead """ - warnings.warn("osub() is deprecated. Use sub(other, method = 'o') instead.", DeprecationWarning) + warnings.warn( + "osub() is deprecated. Use sub(other, method = 'o') instead.", + DeprecationWarning, + ) return Interval(self.left - other.right, self.right - other.left) def mul(self, other, method=None): """ .. _interval.mul: - + Multiplies self by other. **Args**: - + ``other``: The interval or numeric value to be multiplied. This value must be transformable into an Interval object. **Methods**: ``p``: perfect arithmetic :math:`[a,b],[c,d] = [a * c, b * d]` - + ``o``: opposite arithmetic :math:`[a,b],[c,d] = [a * d, b * c]` - + None, i, f - Standard interval arithmetic is used. **Returns**: Interval: The result of the multiplication. - + """ if not isinstance(other, Interval): - if other.__class__.__name__ in ['Pbox','Cbox']: + if other.__class__.__name__ in ["Pbox", "Cbox"]: if method is None: - method = 'f' - return other.mul(self, method = method) + method = "f" + return other.mul(self, method=method) try: other = Interval(other) except: - raise TypeError(f"Multiplication not supported between instances of Interval and {type(other)}") + raise TypeError( + f"Multiplication not supported between instances of Interval and {type(other)}" + ) - if method == 'p': + if method == "p": return Interval(self.left * other.left, self.right * other.right) - elif method == 'o': + elif method == "o": return Interval(self.left * other.right, self.right * other.left) else: return self.__mul__(other) - def pmul(self,other): + def pmul(self, other): """ .. warning:: Depreciated use self.mul(other, method = 'p') instead """ - warnings.warn("pmul() is deprecated. Use mul(other, method = 'p') instead.", DeprecationWarning) + warnings.warn( + "pmul() is deprecated. Use mul(other, method = 'p') instead.", + DeprecationWarning, + ) return Interval(self.left * other.left, self.right * other.right) - - - def omul(self,other): - """ + def omul(self, other): + """ .. warning:: Depreciated use self.mul(other, method = 'o') instead """ - warnings.warn("omul() is deprecated. Use mul(other, method = 'o') instead.", DeprecationWarning) + warnings.warn( + "omul() is deprecated. Use mul(other, method = 'o') instead.", + DeprecationWarning, + ) return Interval(self.left * other.right, self.right * other.left) - + def div(self, other, method=None): - ''' + """ .. _interval.div: - + Divides self by other - - + + If :math:`0 \\in other` it returns a division by zero error - + **Args**: - + ``other`` (Interval or numeric): The interval or numeric value to be multiplied. This value must be transformable into an Interval object. **Methods**: - + ``p``: perfect arithmetic :math:`[a,b],[c,d] = [a * 1/c, b * 1/d]` - + ``o``: opposite arithmetic :math:`[a,b],[c,d] = [a * 1/d, b * 1/c]` - + ``None``, ``i``, ``f`` - Standard interval arithmetic is used. - + .. admonition:: Implementation - + >>> self.add(1/other, method = method) - + .. error:: - + If :math:`0 \\in [a,b]` it returns a division by zero error - - ''' + + """ if not isinstance(other, Interval): - if other.__class__.__name__ in ['Pbox','Cbox']: + if other.__class__.__name__ in ["Pbox", "Cbox"]: if method is None: - return other.recip().mul(self, method = 'f') - elif method == 'o': - return other.recip().mul(self, method = 'p') - elif method == 'p': - return other.recip().mul(self, method = 'o') - else: - return other.recip().mul(self, method = method) + return other.recip().mul(self, method="f") + elif method == "o": + return other.recip().mul(self, method="p") + elif method == "p": + return other.recip().mul(self, method="o") + else: + return other.recip().mul(self, method=method) try: other = Interval(other) except: - raise TypeError(f"Division not supported between instances of Interval and {type(other)}") - - if method == 'o': + raise TypeError( + f"Division not supported between instances of Interval and {type(other)}" + ) + + if method == "o": return Interval(self.left / other.right, self.right / other.left) - elif method == 'p': + elif method == "p": return Interval(self.left / other.left, self.right / other.right) - - return self.mul(1/other, method=method) - - def pdiv(self,other): + + return self.mul(1 / other, method=method) + + def pdiv(self, other): """ .. warning:: Depreciated use self.div(other, method = 'p') instead """ - warnings.warn("pdiv() is deprecated. Use div(other, method = 'p') instead.", DeprecationWarning) - return Interval(self.left / other.left, self.right / other.right) + warnings.warn( + "pdiv() is deprecated. Use div(other, method = 'p') instead.", + DeprecationWarning, + ) + return Interval(self.left / other.left, self.right / other.right) - def odiv(self,other): + def odiv(self, other): """ .. warning:: Depreciated use self.div(other, method = 'o') instead - + """ - return Interval(self.left / other.right, self.right / other.left) + return Interval(self.left / other.right, self.right / other.left) def recip(self): """ Calculates the reciprocle of the interval. - + **Returns**: - + ``Interval``: Equal to :math:`[1/b,1/a]` - + **Example**: - + >>> pba.Interval(2,4).recip() Interval [0.25, 0.5] - + .. error:: If :math:`0 \\in [a,b]` it returns a division by zero error - + """ if self.straddles_zero(): # Cant divide by zero raise ZeroDivisionError() - elif 1/self.hi() < 1/self.lo(): - return Interval(1/self.hi(), 1/self.lo()) + elif 1 / self.hi() < 1 / self.lo(): + return Interval(1 / self.hi(), 1 / self.lo()) else: - return Interval(1/self.lo(), 1/self.hi()) + return Interval(1 / self.lo(), 1 / self.hi()) - def equiv(self,other: "Interval") -> bool: + def equiv(self, other: "Interval") -> bool: """ - Checks whether two intervals are equivalent. - + Checks whether two intervals are equivalent. + **Parameters**: - + ``other``: The interval to check against. - + **Returns** ``True`` **if**: - + ``self.left == other.right`` and ``self.right == other.right`` - + ``False`` otherwise. - + .. error:: - + ``TypeError``: If ``other`` is not an instance of ``Interval`` - + .. seealso:: :func:`~logical.is_same_as` - + **Examples**: - + >>> a = Interval(0,1) >>> b = Interval(0.5,1.5) >>> c = I(0,1) @@ -701,13 +747,13 @@ def equiv(self,other: "Interval") -> bool: False >>> a.equiv(c) True - + """ - if not isinstance(other,Interval): + if not isinstance(other, Interval): raise TypeError(f"Needs to be an instance of Interval not {type(other)}") - - return (self.left == other.left and self.right == other.right) - + + return self.left == other.left and self.right == other.right + def lo(self): """ **Returns**: @@ -731,145 +777,150 @@ def hi(self): """ return self.right - + def width(self) -> float: """ **Returns**: - + ``float``: The width of the interval, :math:`\mathrm{right} - \mathrm{left}` - + **Example**: - + >>> pba.Interval(0,3).width() 3 - + """ return self.right - self.left def halfwidth(self) -> float: """ **Returns**: - + ``float``: The half-width of the interval, :math:`(\mathrm{right} - \mathrm{left})/2` - + **Example**: - + >>> pba.Interval(0,3).halfwidth() 1.5 - + .. admonition:: Implementation - + >>> self.width()/2 - + """ - return self.width()/2 - + return self.width() / 2 + def midpoint(self) -> float: """ **Returns**: - + ``float``: The midpoint of the interval, :math:`(\mathrm{right} + \mathrm{left})/2` - + **Example**: - + >>> pba.Interval(0,2).midpoint() 1.0 """ - - return (self.left+self.right)/2 - + + return (self.left + self.right) / 2 + def to_logical(self): - ''' + """ Turns the interval into a logical interval, this is done by chacking the truth value of the ends of the interval - + **Returns**: - + ``Logical``: The logical interval - + .. admonition:: Implementation - + >>> left = self.left.__bool__() >>> right = self.right.__bool__() >>> Logical(left,right) - - - ''' + + + """ from .logical import Logical - return Logical(self.left.__bool__(),self.right.__bool__()) - + + return Logical(self.left.__bool__(), self.right.__bool__()) def env(self, other: Union[list, "Interval"]) -> "Interval": - ''' + """ Calculates the envelope between two intervals - + **Parameters**: ``other`` : Interval or list. The interval to envelope with self - + .. hint:: If other is a list then the envelope is calculated between self and each element of the list. In this case the envelope is calculated recursively and pba.envelope() may be more efficient. - + .. important:: If other is a Pbox then ``Pbox.env()`` is called - + .. seealso:: - + `pba.core.envelope`_ - + `pba.pbox.Pbox.env`_ - + **Returns**: ``Interval``: The envelope of self and other - - ''' - if isinstance(other, [list,tuple]): + + """ + if isinstance(other, [list, tuple]): e = self for o in other: e = e.env(o) return e - - if not isinstance(other,Interval): - if other.__class__.__name__ in ['Pbox','Cbox']: + + if not isinstance(other, Interval): + if other.__class__.__name__ in ["Pbox", "Cbox"]: return other.env(self) else: try: other = Interval(other) except: - raise TypeError(f"env() not supported between instances of Interval and {type(other)}") - - return Interval(min(self.left,other.left),max(self.right,other.right)) - - - def straddles(self,N: Union[int,float,'Interval'], endpoints:bool = True) -> bool: + raise TypeError( + f"env() not supported between instances of Interval and {type(other)}" + ) + + return Interval(min(self.left, other.left), max(self.right, other.right)) + + def straddles( + self, N: Union[int, float, "Interval"], endpoints: bool = True + ) -> bool: """ .. _interval.straddles: - + **Parameters**: - + ``N``: Number to check. If N is an interval checks whether the whole interval is within self. - + ``endpoints``: Whether to include the endpoints within the check **Returns** ``True`` **if**: - :math:`\mathrm{left} \leq N \leq \mathrm{right}` (Assuming ``endpoints=True``). - - For interval values. :math:`\mathrm{left} \leq N.left \leq \mathrm{right}` and :math:`\mathrm{left} \leq N.right \leq \mathrm{right}` (Assuming ``endpoints=True``). - + :math:`\mathrm{left} \leq N \leq \mathrm{right}` (Assuming ``endpoints=True``). + + For interval values. :math:`\mathrm{left} \leq N.left \leq \mathrm{right}` and :math:`\mathrm{left} \leq N.right \leq \mathrm{right}` (Assuming ``endpoints=True``). + ``False`` otherwise. - + .. tip:: ``N in self`` is equivalent to ``self.straddles(N)`` """ - if isinstance(N,Interval): - return self.straddles(N.left,endpoints) and self.straddles(N.right,endpoints) - + if isinstance(N, Interval): + return self.straddles(N.left, endpoints) and self.straddles( + N.right, endpoints + ) + if endpoints: if self.left <= N and self.right >= N: return True @@ -879,51 +930,58 @@ def straddles(self,N: Union[int,float,'Interval'], endpoints:bool = True) -> boo return False - def straddles_zero(self,endpoints = True): + def straddles_zero(self, endpoints=True): """ Checks whether :math:`0` is within the interval - + .. admonition:: Implementation - + Equivalent to ``self.straddles(0,endpoints)`` - + .. seealso:: interval.straddles_ - - """ - return self.straddles(0,endpoints) + """ + return self.straddles(0, endpoints) - def intersection(self, other: Union["Interval",list]) -> "Interval": - ''' + def intersection(self, other: Union["Interval", list]) -> "Interval": + """ Calculates the intersection between intervals - + **Parameters**: - ``other``: The interval to intersect with self. If an interval is not given will try to cast as an interval. If a list is given will calculate the intersection between self and each element of the list. + ``other``: The interval to intersect with self. If an interval is not given will try to cast as an interval. If a list is given will calculate the intersection between self and each element of the list. **Returns**: - + ``Interval``: The intersection of self and other. If no intersection is found returns ``None`` - + **Example**: >>> a = Interval(0,1) >>> b = Interval(0.5,1.5) >>> a.intersection(b) Interval [0.5, 1] - - - ''' + + + """ if isinstance(other, Interval): if self.straddles(other): - return I(max([x.left for x in [self, other]]), min([x.right for x in [self, other]])) + return I( + max([x.left for x in [self, other]]), + min([x.right for x in [self, other]]), + ) else: return None elif isinstance(other, list): if all([self.straddles(o) for o in other]): - assert all([isinstance(o, Interval) for o in other]), 'All intersected objects must be intervals' - return I(max([x.left for x in [self] + other]), min([x.right for x in [self] + other])) + assert all( + [isinstance(o, Interval) for o in other] + ), "All intersected objects must be intervals" + return I( + max([x.left for x in [self] + other]), + min([x.right for x in [self] + other]), + ) else: return None else: @@ -935,118 +993,121 @@ def intersection(self, other: Union["Interval",list]) -> "Interval": def exp(self): lo = np.exp(self.left) hi = np.exp(self.right) - return Interval(lo,hi) - + return Interval(lo, hi) + def log(self): lo = np.log(self.left) hi = np.log(self.right) - return Interval(lo,hi) - + return Interval(lo, hi) + def sqrt(self): if self.left >= 0: - return Interval(np.sqrt(self.left),np.sqrt(self.right)) + return Interval(np.sqrt(self.left), np.sqrt(self.right)) else: print("RuntimeWarning: invalid value encountered in sqrt") - return Interval(np.nan,np.sqrt(self.right)) - + return Interval(np.nan, np.sqrt(self.right)) + # def sin(self): # return Interval(np.sin(self.left),np.sin(self.right)) # def cos(self): # return Interval(np.cos(self.left),np.cos(self.right)) # def tan(self): # return Interval(np.tan(self.left),np.tan(self.right)) - + def sample(self, seed=None, numpy_rng: np.random.Generator = None) -> float: """ Generate a random sample within the interval. **Parameters**: - + ``seed`` (int, optional): Seed value for random number generation. Defaults to None. - + ``numpy_rng`` (numpy.random.Generator, optional): Numpy random number generator. Defaults to None. - + **Returns**: - + ``float``: Random sample within the interval. - + .. admonition:: Implementation - + If ``numpy_rng`` is given: - + >>> numpy_rng.uniform(self.left, self.right) - + Otherwise the following is used: - + >>> import random >>> random.seed(seed) >>> self.left + random.random() * self.width() - + **Examples**: >>> pba.Interval(0,1).sample() 0.6160988752201705 >>> pba.I(0,1).sample(seed = 1) 0.13436424411240122 - + If a numpy random number generator is given then it is used instead of the default python random number generator. It has to be initialised first. - + >>> import numpy as np >>> rng = np.random.default_rng(seed = 0) >>> pba.I(0,1).sample(numpy_rng = rng) 0.6369616873214543 - - + + """ - + if numpy_rng is not None: return numpy_rng.uniform(self.left, self.right) - + if seed is not None: r.seed(seed) return self.left + r.random() * self.width() - + + # Alias I = Interval + def PM(x, hw): """ Create an interval centered around x with a half-width of hw. **Parameters**: - + ``x`` (float): The center value of the interval. - + ``hw`` (float): The half-width of the interval. **Returns**: - + ``Interval``: An interval object with lower bound x-hw and upper bound x+hw. .. error:: - + ``ValueError``: If hw is less than 0. - + Example: >>> pba.pm(0, 1) Interval [-1, 1] - + """ if hw < 0: raise ValueError("hw must be greater than or equal to 0") - + return Interval(x - hw, x + hw) + def pm_repr(): - ''' + """ .. _interval.pm_repr: - + Modifies the interval class to display the interval in [midpoint ± half-width] format. - + **Example**: - + .. code-block:: python - + >>> import pba >>> pba.interval.pm_repr() >>> a = pba.Interval(0,1) # defined using left and right. This cannot be overriden. @@ -1055,31 +1116,32 @@ def pm_repr(): Interval [0.5 ± 0.5] >>> print(b) Interval [0 ± 1] - + .. seealso:: :func:`~pba.lr_interval_repr` - - ''' - + + """ + def new__repr__(self): - return f'Interval [{self.midpoint():g} ± {self.halfwidth():g}]' - + return f"Interval [{self.midpoint():g} ± {self.halfwidth():g}]" + Interval.__repr__ = new__repr__ Interval.__str__ = new__repr__ + def lr_repr(): - ''' + """ .. _interval.lr_repr: - + Modifies the interval class to display the interval in [left, right] format. - + .. note:: This function primarily exists to undo the effects of ``pm_interval_repr()``. By default the interval class displays in this format. - + **Example**: - + .. code-block:: python - + >>> import pba >>> pba.interval.pm_repr() >>> a = pba.Interval(0,1) # defined using left and right, this cannot be overriden. @@ -1089,11 +1151,11 @@ def lr_repr(): >>> b = pba.PM(0,1) # defined using midpoint and half-width >>> print(b) Interval [-1,1] - - ''' - + + """ + def new__repr__(self): - return f'Interval [{self.left:g}, {self.right:g}]' - + return f"Interval [{self.left:g}, {self.right:g}]" + Interval.__repr__ = new__repr__ - Interval.__str__ = new__repr__ \ No newline at end of file + Interval.__str__ = new__repr__ diff --git a/pba/logical.py b/pba/logical.py index f6b7e5e..1f0bd2b 100644 --- a/pba/logical.py +++ b/pba/logical.py @@ -5,22 +5,22 @@ from .pbox import Pbox from .interval import Interval - class Logical(Interval): - ''' + """ Represents a logical value that can be either True or False or dunno ([False,True]). Inherits from the Interval class. **Attributes**: - + ``left`` (``bool``): The left endpoint of the logical value. - + ``right`` (``bool``): The right endpoint of the logical value. - ''' - def __init__(self, left: bool ,right: bool = None): + """ + + def __init__(self, left: bool, right: bool = None): super().__init__(left, right) def __bool__(self): @@ -29,19 +29,21 @@ def __bool__(self): if self.left == 1 and self.right == 1: return True else: - print('WARNING: Truth value of Logical is ambiguous, use pba.sometime or pba.always') + print( + "WARNING: Truth value of Logical is ambiguous, use pba.sometime or pba.always" + ) return True def __repr__(self): if self.left == 0 and self.right == 0: - return 'False' + return "False" elif self.left == 1 and self.right == 1: - return 'True' + return "True" else: - return 'Dunno [False,True]' + return "Dunno [False,True]" __str__ = __repr__ - + def __invert__(self): if self.left == 0 and self.right == 0: return Logical(True, True) @@ -50,26 +52,32 @@ def __invert__(self): else: return self -def is_same_as(a: Union['Pbox', 'Interval'], b: Union['Pbox', 'Interval'], deep = False, exact_pbox = False): + +def is_same_as( + a: Union["Pbox", "Interval"], + b: Union["Pbox", "Interval"], + deep=False, + exact_pbox=False, +): """ Check if two objects of type 'Pbox' or 'Interval' are equal. **Parameters**: - + ``a``: The first object to be compared. - + ``b``: The second object to be compared. - + ``deep``: If True, performs a deep comparison, considering object identity. If False, performs a shallow comparison based on object attributes. Defaults to False. - + ``exact_pbox``: If True, performs a deep comparison of p-boxes, considering all attributes. If False, performs a shallow comparison of p-boxes, considering only the left and right attributes. Defaults to False. **Returns** ``True`` **if**: - + ``bool``: True if the objects have identical parameters. For Intervals this means that left and right are the same for both a and b. For p-boxes checks whether all p-box attributes are the same. If deep is True, checks whether the objects have the same id. - + **Examples**: - + >>> a = Interval(0, 2) >>> b = Interval(0, 2) >>> c = Interval(1, 3) @@ -77,9 +85,9 @@ def is_same_as(a: Union['Pbox', 'Interval'], b: Union['Pbox', 'Interval'], deep True >>> is_same_as(a, c) False - + For p-boxes: - + >>> a = pba.N([0,1],1) >>> b = pba.N([0,1],1) >>> c = pba.N([0,1],2) @@ -93,12 +101,12 @@ def is_same_as(a: Union['Pbox', 'Interval'], b: Union['Pbox', 'Interval'], deep False >>> is_same_as(e, f, exact_pbox = False) True - - + + """ - if not isinstance(a,(Interval, Pbox)) or not isinstance(b,(Interval, Pbox)): + if not isinstance(a, (Interval, Pbox)) or not isinstance(b, (Interval, Pbox)): return a == b - + if deep: if id(a) == id(b): return True @@ -107,45 +115,39 @@ def is_same_as(a: Union['Pbox', 'Interval'], b: Union['Pbox', 'Interval'], deep else: if a.__class__.__name__ != b.__class__.__name__: return False - elif isinstance(a,Pbox): - + elif isinstance(a, Pbox): + if exact_pbox: if ( - np.array_equal(a.left , b.left ) and - np.array_equal(a.right, b.right) and - a.steps == b.steps and - a.shape == b.shape and - a.mean.left == b.mean.left and - a.mean.right == b.mean.right and - a.var.left == b.var.left and - a.var.right == b.var.right - ): + np.array_equal(a.left, b.left) + and np.array_equal(a.right, b.right) + and a.steps == b.steps + and a.shape == b.shape + and a.mean.left == b.mean.left + and a.mean.right == b.mean.right + and a.var.left == b.var.left + and a.var.right == b.var.right + ): return True else: return False else: - if ( - np.array_equal(a.left, b.left) and - np.array_equal(a.right, b.right) - ): + if np.array_equal(a.left, b.left) and np.array_equal(a.right, b.right): return True else: return False - - - elif isinstance(a,Interval): - if ( - a.left == b.left and - a.right == b.right - ): + + elif isinstance(a, Interval): + if a.left == b.left and a.right == b.right: return True else: return False - -def always(logical: Union[Logical, Interval , Number, bool]) -> bool: + + +def always(logical: Union[Logical, Interval, Number, bool]) -> bool: """ Checks whether the logical value is always true. i.e. Every value from one interval or p-box is always greater than any other values from another. - + This function takes either a Logical object, an interval or a float as input and checks if both the left and right attributes of the Logical object are True. If an interval is provided, it checks that both the left and right attributes of the Logical object are 1. @@ -158,11 +160,11 @@ def always(logical: Union[Logical, Interval , Number, bool]) -> bool: ``bool``: True if both sides of the logical condition are True or if the float value is equal to 1, False otherwise. .. error:: - + ``TypeError``: If the input is not an instance of Interval, Logical or a numeric value. - + ``ValueError``: If the input float is not between 0 and 1 or the interval contains values outside of [0,1] - + **Examples**: >>> a = Interval(0, 2) >>> b = Interval(1, 3) @@ -171,10 +173,10 @@ def always(logical: Union[Logical, Interval , Number, bool]) -> bool: False >>> always(a < c) True - + """ - - if isinstance(logical,Logical): + + if isinstance(logical, Logical): if logical.left and logical.right: return True @@ -183,48 +185,52 @@ def always(logical: Union[Logical, Interval , Number, bool]) -> bool: elif isinstance(logical, Interval): if logical.left < 0 or logical.right > 1: - raise ValueError("If interval values needs to be between 0 and 1 (inclusive)") + raise ValueError( + "If interval values needs to be between 0 and 1 (inclusive)" + ) if logical.left == 1 and logical.right == 1: return True else: return False - elif isinstance(logical, (bool,Number)): - + elif isinstance(logical, (bool, Number)): + if logical < 0 or logical > 1: raise ValueError("If numeric input needs to be between 0 and 1 (inclusive)") if logical == 1: return True - else: False - + else: + False + else: raise TypeError("Input must be a Logical, Interval or a numeric value.") + def never(logical: Logical) -> bool: """ Checks whether the logical value is always true. i.e. Every value from one interval or p-box is always less than any other values from another. - + This function takes either a Logical object, an interval or a float as input and checks if both the left and right attributes of the Logical object are False. If an interval is provided, it checks that both the left and right attributes of the Logical object are 0. If a numeric value is provided, it checks if the is equal to 0. **Parameters**: - + ``logical`` (``Logical``, ``Interval`` , ``Number``): An object representing a logical condition with 'left' and 'right' attributes, or a number between 0 and 1. **Returns**: - + ``bool``: True if both sides of the logical condition are True or if the float value is equal to 0, False otherwise. .. error:: - + ``TypeError``: If the input is not an instance of Interval, Logical or a numeric value. - + ``ValueError``: If the input float is not between 0 and 1 or the interval contains values outside of [0,1] - + **Examples**: - + >>> a = Interval(0, 2) >>> b = Interval(1, 3) >>> c = Interval(4, 5) @@ -232,10 +238,10 @@ def never(logical: Logical) -> bool: False >>> never(a < c) True - + """ - - if isinstance(logical,Logical): + + if isinstance(logical, Logical): if not logical.left and not logical.right: return True @@ -244,48 +250,52 @@ def never(logical: Logical) -> bool: elif isinstance(logical, Interval): if logical.left < 0 or logical.right > 1: - raise ValueError("If interval values needs to be between 0 and 1 (inclusive)") + raise ValueError( + "If interval values needs to be between 0 and 1 (inclusive)" + ) if logical.left == 0 and logical.right == 0: return True else: return False - elif isinstance(logical, (bool,Number)): - + elif isinstance(logical, (bool, Number)): + if logical < 0 or logical > 1: raise ValueError("If numeric input needs to be between 0 and 1 (inclusive)") if logical == 0: return True - else: False - + else: + False + else: raise TypeError("Input must be a Logical, Interval or a numeric value.") + def sometimes(logical: Logical) -> bool: """ Checks whether the logical value is sometimes true. i.e. There exists one value from one interval or p-box is less than a values from another. - + This function takes either a Logical object, an interval or a float as input and checks if either the left and right attributes of the Logical object are True. If an interval is provided, it that both endpoints are not 0. If a numeric value is provided, it checks if the is not equal to 0. **Parameters**: - + ``logical`` (``Logical``, ``Interval`` , ``Number``): An object representing a logical condition with 'left' and 'right' attributes, or a number between 0 and 1. **Returns**: - + ``bool``: True if both sides of the logical condition are True or if the float value is equal to 0, False otherwise. .. error:: - + ``TypeError``: If the input is not an instance of Interval, Logical or a numeric value. - + ``ValueError``: If the input float is not between 0 and 1 or the interval contains values outside of [0,1] - + **Examples**: - + >>> a = pba.Interval(0, 2) >>> b = pba.Interval(1, 4) >>> c = pba.Interval(3, 5) @@ -295,10 +305,10 @@ def sometimes(logical: Logical) -> bool: True >>> pba.sometimes(c < b) True - + """ - - if isinstance(logical,Logical): + + if isinstance(logical, Logical): if not logical.left and not logical.right: return False @@ -307,48 +317,51 @@ def sometimes(logical: Logical) -> bool: elif isinstance(logical, Interval): if logical.left < 0 or logical.right > 1: - raise ValueError("If interval values needs to be between 0 and 1 (inclusive)") + raise ValueError( + "If interval values needs to be between 0 and 1 (inclusive)" + ) if logical.left != 0 or logical.right != 0: return True else: return False - elif isinstance(logical, (bool,Number)): - + elif isinstance(logical, (bool, Number)): + if 0 < logical <= 1: return True elif logical == 0: return False else: raise ValueError("If numeric input needs to be between 0 and 1 (inclusive)") - + else: raise TypeError("Input must be a Logical, Interval or a numeric value.") + def xtimes(logical: Logical) -> bool: """ Checks whether the logical value is exclusively sometimes true. i.e. There exists one value from one interval or p-box is less than a values from another but it is not always the case. - + This function takes either a Logical object, an interval or a float as input and checks that the left value is False and the right value is True If an interval is provided, it that both endpoints are not 0 or 1. If a numeric value is provided, it checks if the is not equal to 0 or 1. **Parameters**: - + ``logical`` (``Logical``, ``Interval`` , ``Number``): An object representing a logical condition with 'left' and 'right' attributes, or a number between 0 and 1. **Returns**: - + ``bool``: True if both sides of the logical condition are True or if the float value is equal to 0, False otherwise. .. error:: - + ``TypeError``: If the input is not an instance of Interval, Logical or a numeric value. - + ``ValueError``: If the input float is not between 0 and 1 or the interval contains values outside of [0,1] - + **Examples**: - + >>> a = pba.Interval(0, 2) >>> b = pba.Interval(2, 4) >>> c = pba.Interval(2.5,3.5) @@ -358,10 +371,10 @@ def xtimes(logical: Logical) -> bool: False >>> pba.xtimes(c < b) True - + """ - - if isinstance(logical,Logical): + + if isinstance(logical, Logical): if not logical.left and logical.right: return True @@ -370,21 +383,22 @@ def xtimes(logical: Logical) -> bool: elif isinstance(logical, Interval): if logical.left < 0 or logical.right > 1: - raise ValueError("If interval values needs to be between 0 and 1 (inclusive)") + raise ValueError( + "If interval values needs to be between 0 and 1 (inclusive)" + ) if logical.left != 0 and logical.right != 1: return True else: return False - elif isinstance(logical, (bool,Number)): - - if Interval(0,1).straddles(logical, endpoints=False): + elif isinstance(logical, (bool, Number)): + + if Interval(0, 1).straddles(logical, endpoints=False): return True elif logical == 0 or logical == 1: return False else: raise ValueError("If numeric input needs to be between 0 and 1 (inclusive)") - + else: raise TypeError("Input must be a Logical, Interval or a numeric value.") - diff --git a/pba/pbox.py b/pba/pbox.py index 0596b42..f128900 100644 --- a/pba/pbox.py +++ b/pba/pbox.py @@ -10,84 +10,100 @@ from .interval import Interval from .copula import Copula -__all__ = [ - 'Pbox', - 'mixture', - 'truncate', - 'imposition', - 'NotIncreasingError' -] +__all__ = ["Pbox", "mixture", "truncate", "imposition", "NotIncreasingError"] ### Local functions ### def _check_increasing(arr): return np.all(np.diff(arr) >= 0) + def _interpolate(left, right, steps, method): - if method == 'linear': - nleft = np.interp(np.linspace(0,1,steps), np.linspace(0,1,len(left)), left) - nright = np.interp(np.linspace(0,1,steps), np.linspace(0,1,len(right)), right) - elif method == 'cubicspline': - nleft = spi.CubicSpline(np.linspace(0,1,len(left)), left)(np.linspace(0,1,steps)) - nright = spi.CubicSpline(np.linspace(0,1,len(right)), right)(np.linspace(0,1,steps)) - elif method == 'step': - percentiles = {i:Interval(j,k) for i,j,k in zip(np.linspace(0,1,len(left)), left, right)} + if method == "linear": + nleft = np.interp(np.linspace(0, 1, steps), np.linspace(0, 1, len(left)), left) + nright = np.interp( + np.linspace(0, 1, steps), np.linspace(0, 1, len(right)), right + ) + elif method == "cubicspline": + nleft = spi.CubicSpline(np.linspace(0, 1, len(left)), left)( + np.linspace(0, 1, steps) + ) + nright = spi.CubicSpline(np.linspace(0, 1, len(right)), right)( + np.linspace(0, 1, steps) + ) + elif method == "step": + percentiles = { + i: Interval(j, k) + for i, j, k in zip(np.linspace(0, 1, len(left)), left, right) + } nleft = [] nright = [] - for i in np.linspace(0,1,steps): + for i in np.linspace(0, 1, steps): left_key = max(key for key in percentiles.keys() if key <= i) - right_key = min(key for key in percentiles.keys() if key >= i) + right_key = min(key for key in percentiles.keys() if key >= i) nleft.append(percentiles[left_key].left) nright.append(percentiles[right_key].right) else: - raise ValueError("Invalid interpolation method. Must be one of: linear, cubicspline, step") - + raise ValueError( + "Invalid interpolation method. Must be one of: linear, cubicspline, step" + ) + return nleft, nright + class NotIncreasingError(Exception): pass + def _check_div_by_zero(pbox): if 0 <= min(pbox.left) and 0 >= max(pbox.right): raise DivisionByZero("Pbox contains 0") -def _interval_list_to_array(l, left = True): + +def _interval_list_to_array(l, left=True): if left: f = lambda x: x.left if isinstance(x, Interval) else x - else: # must be right + else: # must be right f = lambda x: x.right if isinstance(x, Interval) else x - + return np.array([f(i) for i in l]) - + + def _check_steps(a, b): if a.steps > b.steps: - warn("Pboxes have different number of steps. Interpolating {b.__name__} to {a.steps} steps") - b.interpolate(a.steps, inplace = True) + warn( + "Pboxes have different number of steps. Interpolating {b.__name__} to {a.steps} steps" + ) + b.interpolate(a.steps, inplace=True) elif a.steps < b.steps: - warn("Pboxes have different number of steps. Interpolating {a.__name__} to {b.steps} steps") - a.interpolate(b.steps, inplace = True) - + warn( + "Pboxes have different number of steps. Interpolating {a.__name__} to {b.steps} steps" + ) + a.interpolate(b.steps, inplace=True) + return a, b - + + def _check_moments(left, right, steps, mean, var): - - def _sideVariance(w,mu): - if not isinstance(w, np.ndarray): w = np.array(w) - if mu is None: mu = np.mean(w) + + def _sideVariance(w, mu): + if not isinstance(w, np.ndarray): + w = np.array(w) + if mu is None: + mu = np.mean(w) return max(0, np.mean((w - mu) ** 2)) - + cmean = Interval(np.mean(left), np.mean(right)) - if not mean.equiv(cmean) and not mean.equiv(Interval(-np.inf,np.inf)): + if not mean.equiv(cmean) and not mean.equiv(Interval(-np.inf, np.inf)): warn("Mean specified does not match calculated mean. Using calculated mean") print(f"Specified mean: {mean}, calculated mean: {cmean}") mean = cmean if np.any(np.isinf(left)) or np.any(np.isinf(right)): - cvar = Interval(0, np.inf) + cvar = Interval(0, np.inf) if np.all(right[0] == right) and np.all(left[0] == left): - cvar = Interval(0, (right[0] - left[0]) ** (2 / 4)) - + cvar = Interval(0, (right[0] - left[0]) ** (2 / 4)) vr = _sideVariance(left, np.mean(left)) w = np.copy(left) @@ -113,7 +129,7 @@ def _sideVariance(w,mu): here = w[i] if 1 < i: - for j in reversed(range(i-1)): + for j in reversed(range(i - 1)): if w[i] < w[j]: w[j] = here @@ -126,214 +142,217 @@ def _sideVariance(w,mu): cvar = Interval(vl, vr) - if not var.equiv(cvar) and not var.equiv(Interval(0,np.inf)): - warn("Variance specified does not match calculated variance. Using calculated variance") + if not var.equiv(cvar) and not var.equiv(Interval(0, np.inf)): + warn( + "Variance specified does not match calculated variance. Using calculated variance" + ) print(f"Specified variance: {var}, calculated variance: {cvar}") var = cvar return cmean, cvar + ### Arithmetic Functions ### -#** This is prefered as means can only test once for +, -, *, /, **# -def _arithmetic(a, b, method, op, enforce_steps = True, interpolation_method = 'linear'): - #* If enforce_steps is True, the number of steps in the returned p-box is the maximum of the number of steps in a and b. +# ** This is prefered as means can only test once for +, -, *, /, **# +def _arithmetic(a, b, method, op, enforce_steps=True, interpolation_method="linear"): + # * If enforce_steps is True, the number of steps in the returned p-box is the maximum of the number of steps in a and b. - if b.__class__.__name__ == 'Interval': - other = Pbox(other, steps = a.steps) + if b.__class__.__name__ == "Interval": + other = Pbox(other, steps=a.steps) - if b.__class__.__name__ == 'Pbox': + if b.__class__.__name__ == "Pbox": a, b = _check_steps(a, b) - if method == 'f': + if method == "f": with catch_warnings(): - simplefilter('ignore') + simplefilter("ignore") nleft, nright = _f_arithmetic(a, b, op) - elif method == 'p': + elif method == "p": with catch_warnings(): - simplefilter('ignore') + simplefilter("ignore") nleft, nright = _p_arithmetic(a, b, op) - elif method == 'o': + elif method == "o": with catch_warnings(): - simplefilter('ignore') + simplefilter("ignore") nleft, nright = _o_arithmetic(a, b, op) - elif method == 'i': + elif method == "i": with catch_warnings(): - simplefilter('ignore') + simplefilter("ignore") nleft, nright = _i_arithmetic(a, b, op) - else: raise ArithmeticError("Calculation method unkown") - + else: + raise ArithmeticError("Calculation method unkown") + nleft.sort() nright.sort() if enforce_steps: # Steps needs to match the maximum number of steps in a return Pbox( - left = nleft, - right = nright, - steps = max(a.steps, b.steps), - interpolation=interpolation_method + left=nleft, + right=nright, + steps=max(a.steps, b.steps), + interpolation=interpolation_method, ) else: - return Pbox( - left = nleft, - right = nright - ) + return Pbox(left=nleft, right=nright) else: try: # Try adding constant - if a.shape in ['uniform','normal','cauchy','triangular','skew-normal']: + if a.shape in ["uniform", "normal", "cauchy", "triangular", "skew-normal"]: s = a.shape else: - s = '' + s = "" return Pbox( - left = op(a.left,other), - right = op(a.right,other), - shape = s, - steps = a.steps + left=op(a.left, other), right=op(a.right, other), shape=s, steps=a.steps ) except: return NotImplemented -def _f_arithmetic(a,b,op): - nleft = np.empty(a.steps) +def _f_arithmetic(a, b, op): + nleft = np.empty(a.steps) nright = np.empty(a.steps) - for i in range(0,a.steps): + for i in range(0, a.steps): j = np.array(range(i, a.steps)) - k = np.array(range(a.steps - 1, i-1, -1)) + k = np.array(range(a.steps - 1, i - 1, -1)) - nright[i] = np.min(op(a.right[j],b.right[k])) + nright[i] = np.min(op(a.right[j], b.right[k])) jj = np.array(range(0, i + 1)) - kk = np.array(range(i, -1 , -1)) + kk = np.array(range(i, -1, -1)) + + nleft[i] = np.max(op(a.left[jj], b.left[kk])) - nleft[i] = np.max(op(a.left[jj],b.left[kk])) - return nleft, nright -def _i_arithmetic(a,b,op): - nleft = np.array([op(i,j) for i,j in itertools.product(a.left,b.left)]) - nright = np.array([op(i,j) for i,j in itertools.product(a.right,b.right)]) - +def _i_arithmetic(a, b, op): + + nleft = np.array([op(i, j) for i, j in itertools.product(a.left, b.left)]) + nright = np.array([op(i, j) for i, j in itertools.product(a.right, b.right)]) + return nleft, nright -def _o_arithmetic(a,b,op): - nleft = op(a.left,np.flip(b.right)) - nright = op(a.right,np.flip(b.left)) + +def _o_arithmetic(a, b, op): + nleft = op(a.left, np.flip(b.right)) + nright = op(a.right, np.flip(b.left)) return nleft, nright -def _p_arithmetic(a,b,op): - nleft = op(a.left,b.left) - nright = op(a.right,b.right) + +def _p_arithmetic(a, b, op): + nleft = op(a.left, b.left) + nright = op(a.right, b.right) return nleft, nright ### Pbox Class ### class Pbox: - r''' + r""" A probability distribution is a mathematical function that gives the probabilities of occurrence for different possible values of a variable. Probability boxes (p-boxes) represent interval bounds on probability distributions. The simplest kind of p-box can be expressed mathematically as - + .. math:: - + \mathcal{F}(x) = [\underline{F}(x),\overline{F}(x)], \ \underline{F}(x)\geq \overline{F}(x)\ \forall x \in \mathbb{R} - where :math:`\underline{F}(x)` is the function that defines the left bound of the p-box and :math:`\overline{F}(x)` defines the right bound of the p-box. In PBA the left and right bounds are each stored as a NumPy array containing the percent point function (the inverse of the cumulative distribution function) for `steps` evenly spaced values between 0 and 1. P-boxes can be defined using all the probability distributions that are available through SciPy's statistics library, + where :math:`\underline{F}(x)` is the function that defines the left bound of the p-box and :math:`\overline{F}(x)` defines the right bound of the p-box. In PBA the left and right bounds are each stored as a NumPy array containing the percent point function (the inverse of the cumulative distribution function) for `steps` evenly spaced values between 0 and 1. P-boxes can be defined using all the probability distributions that are available through SciPy's statistics library, - Naturally, precise probability distributions can be defined in PBA by defining a p-box with precise inputs. This means that within probability bounds analysis probability distributions are considered a special case of a p-box with zero width. Resultantly, all methodology that applies to p-boxes can also be applied to probability distributions. + Naturally, precise probability distributions can be defined in PBA by defining a p-box with precise inputs. This means that within probability bounds analysis probability distributions are considered a special case of a p-box with zero width. Resultantly, all methodology that applies to p-boxes can also be applied to probability distributions. Distribution-free p-boxes can also be generated when the underlying distribution is unknown but parameters such as the mean, variance or minimum/maximum bounds are known. Such p-boxes make no assumption about the shape of the distribution and instead return bounds expressing all possible distributions that are valid given the known information. Such p-boxes can be constructed making use of Chebyshev, Markov and Cantelli inequalities from probability theory. .. attention:: - + It is usually better to define p-boxes using distributions or non-parametric methods (see ). This constructor is provided for completeness and for the construction of p-boxes with precise inputs. - + :arg left: Left bound of the p-box. Can be a list, NumPy array, Interval or numeric type. :arg right: Right bound of the p-box. Can be a list, NumPy array, Interval or numeric type. - :arg steps: Number of steps to discretize the p-box into. Default is None. + :arg steps: Number of steps to discretize the p-box into. Default is None. :arg shape: The shap eof the distribution used to construct the p-box. The shape is defined by the p-box constructor. Default is `None`. :arg mean: Interval containing the mean of the p-box. Default is `Interval(-np.inf,np.inf)`. :arg var: Interval containing the variance of the p-box. Default is `Interval(-np.inf,np.inf)`. :arg interpolation: Interpolation method to use. See interpolation for more details. Default is `linear`. :arg check_moments: If `True`, the mean and variance of the p-box are checked and recalculated if necessary. Default is `True`. - + .. important:: - - If steps is not specified, left and right must be arrays of the same length and steps is set at that value. - + + If steps is not specified, left and right must be arrays of the same length and steps is set at that value. + If steps is specified, both left and right are interpolated to the specified number of steps using the specified interpolation method (see interpolation_). In this case if steps is less than the length of left or right, a warning is raised and steps is set to the length of left or right. .. warning:: - + The statistic values be specified by constructor function are also calculated automatically. If specified values differ from the calculated values, the calculated values are used and a warning is raised. - + If check_moments is set to ``True`` and mean and/or var are specified, if the calculated values differ from the specified values, the calculated values are used and a warning is raised. - + .. error:: - + If the left and right bounds are not increasing, a `NotIncreasingError` is raised. ``ValueError`` is raised if left and right are not the same length. - - - ''' - + + + """ + STEPS = 200 - def __init__(self, - left: Union[list,np.ndarray], - right: Union[list,np.ndarray]=None, - steps: int =None, - shape: str="", - mean: Interval = Interval(-np.inf, np.inf), - var: Interval = Interval(0, np.inf), - mean_left: float = None, - mean_right: float = None, - var_left: float = None, - var_right: float = None, - interpolation: str ='linear', - check_moments: bool = True): - + def __init__( + self, + left: Union[list, np.ndarray], + right: Union[list, np.ndarray] = None, + steps: int = None, + shape: str = "", + mean: Interval = Interval(-np.inf, np.inf), + var: Interval = Interval(0, np.inf), + mean_left: float = None, + mean_right: float = None, + var_left: float = None, + var_right: float = None, + interpolation: str = "linear", + check_moments: bool = True, + ): + #!! TO BE DEPRECATED !!# if mean_left is not None or mean_right is not None: mean = Interval(mean_left, mean_right) if var_left is not None or var_right is not None: var = Interval(var_left, var_right) - + if right is None: right = left - - if not hasattr(left, '__iter__'): + + if not hasattr(left, "__iter__"): raise ValueError("left must be an iterable") - if not hasattr(right, '__iter__'): + if not hasattr(right, "__iter__"): raise ValueError("right must be an iterable") if isinstance(left, Interval): if steps is None: raise ValueError("steps must be specified if left is an Interval") - left = np.array([left.left]*steps) - + left = np.array([left.left] * steps) + if isinstance(right, Interval): if steps is None: raise ValueError("steps must be specified if right is an Interval") - right = np.array([right.right]*steps) + right = np.array([right.right] * steps) if len(left) != len(right): raise ValueError("left and right must be the same number of steps") - + if not _check_increasing(left) or not _check_increasing(right): raise NotIncreasingError("Left and right arrays must be increasing") - - if steps is None: + + if steps is None: steps = len(left) if isinstance(left, list): @@ -342,421 +361,457 @@ def __init__(self, left = np.array(left) if isinstance(right, list): - right = _interval_list_to_array(right, left = False) + right = _interval_list_to_array(right, left=False) elif not isinstance(right, np.ndarray): right = np.array(right) if steps != len(left): - warn(f"Number of steps does not match length of left. Interpolating left/right to {steps} steps") + warn( + f"Number of steps does not match length of left. Interpolating left/right to {steps} steps" + ) left, right = _interpolate(left, right, steps, interpolation) - - l,r = zip(*[(min(i),max(i)) for i in zip(left,right)]) + + l, r = zip(*[(min(i), max(i)) for i in zip(left, right)]) self.left = np.array(l) self.right = np.array(r) self.steps = steps self.shape = shape - - if mean.equiv(Interval(-np.inf,np.inf)) or var.equiv(Interval(-np.inf,np.inf)) or check_moments: - self.mean, self.var = _check_moments(l,r,steps,mean,var) - + if ( + mean.equiv(Interval(-np.inf, np.inf)) + or var.equiv(Interval(-np.inf, np.inf)) + or check_moments + ): + + self.mean, self.var = _check_moments(l, r, steps, mean, var) + else: self.mean = mean self.var = var def __repr__(self): - return f'Pbox: ~ {self.shape}(range={Interval(self.lo(),self.hi()):g}, mean={self.mean:g}, var={self.var:g})' + return f"Pbox: ~ {self.shape}(range={Interval(self.lo(),self.hi()):g}, mean={self.mean:g}, var={self.var:g})" __str__ = __repr__ def __iter__(self): - for val in np.array([self.left,self.right]).flatten(): + for val in np.array([self.left, self.right]).flatten(): yield val - def __neg__(self): - if self.shape in ['uniform','normal','cauchy','triangular','skew-normal']: + if self.shape in ["uniform", "normal", "cauchy", "triangular", "skew-normal"]: s = self.shape else: - s = '' + s = "" return Pbox( - left = sorted(-np.flip(self.right)), - right = sorted(-np.flip(self.left)), - steps = len(self.left), - shape = s, - mean = -self.mean, - var = self.var, - check_moments = False + left=sorted(-np.flip(self.right)), + right=sorted(-np.flip(self.left)), + steps=len(self.left), + shape=s, + mean=-self.mean, + var=self.var, + check_moments=False, ) - def __lt__(self,other): - return self.lt(other, method = 'f') + def __lt__(self, other): + return self.lt(other, method="f") - def __rlt__(self,other): - return self.ge(other, method = 'f') + def __rlt__(self, other): + return self.ge(other, method="f") - def __le__(self,other): - return self.le(other, method = 'f') + def __le__(self, other): + return self.le(other, method="f") - def __rle__(self,other): - return self.gt(other, method = 'f') + def __rle__(self, other): + return self.gt(other, method="f") - def __gt__(self,other): - return self.gt(other, method = 'f') + def __gt__(self, other): + return self.gt(other, method="f") - def __rgt__(self,other): - return self.le(other, method = 'f') + def __rgt__(self, other): + return self.le(other, method="f") - def __ge__(self,other): - return self.ge(other, method = 'f') + def __ge__(self, other): + return self.ge(other, method="f") - def __rge__(self,other): - return self.lt(other, method = 'f') + def __rge__(self, other): + return self.lt(other, method="f") def __and__(self, other): - return self.logicaland(other, method = 'f') + return self.logicaland(other, method="f") - def __rand__(self,other): - return self.logicaland(other, method = 'f') + def __rand__(self, other): + return self.logicaland(other, method="f") def __or__(self, other): - return self.logicalor(other, method = 'f') + return self.logicalor(other, method="f") - def __ror__(self,other): - return self.logicalor(other, method = 'f') + def __ror__(self, other): + return self.logicalor(other, method="f") def __add__(self, other): - return self.add(other, method = 'f') + return self.add(other, method="f") - def __radd__(self,other): - return self.add(other, method = 'f') + def __radd__(self, other): + return self.add(other, method="f") - def __sub__(self,other): - return self.sub(other, method = 'f') + def __sub__(self, other): + return self.sub(other, method="f") - def __rsub__(self,other): - self = - self - return self.add(other, method = 'f') + def __rsub__(self, other): + self = -self + return self.add(other, method="f") - def __mul__(self,other): - return self.mul(other, method = 'f') + def __mul__(self, other): + return self.mul(other, method="f") - def __rmul__(self,other): - return self.mul(other, method = 'f') + def __rmul__(self, other): + return self.mul(other, method="f") - def __pow__(self,other): - return self.pow(other, method='f') - - def __rpow__(self,other): - if not hasattr(other, '__iter__'): + def __pow__(self, other): + return self.pow(other, method="f") + + def __rpow__(self, other): + if not hasattr(other, "__iter__"): other = np.array((other)) - + b = Pbox(other) - return b.pow(self, method='f') - + return b.pow(self, method="f") + def __truediv__(self, other): - return self.div(other, method = 'f') + return self.div(other, method="f") - def __rtruediv__(self,other): + def __rtruediv__(self, other): try: return other * self.recip() except: return NotImplemented - - def lo(self): - ''' + + def lo(self): + """ Returns the left-most value in the interval - ''' + """ return self.left[0] - - def hi(self): - ''' + + def hi(self): + """ Returns the right-most value in the interval - ''' - return self.right[-1] - + """ + return self.right[-1] + def unary(self, func, *args, **kwargs): - ''' + """ Allows for unary operations to be performed on a p-box. This is acheived by applying the function to each interval in the p-box. - + **Arguments:** ``func`` (``function``): Function to apply to each interval in the p-box. ``args`` (``tuple``): Arguments to pass to the function. ``kwargs`` (``dict``): Keyword arguments to pass to the function. - + .. important:: - + The function must accept an Interval object as its first argument and return an Interval object. ``args`` and ``kwargs`` are passed to the function as additional arguments. - + >>> func(Interval(l,r),*args, **kwargs) - + The function must return an Interval object. Behaviour may be unpredictable if the endpoints of the inputted interval do not correspond to the endpoints of the outputted p-box. - - ''' - ints = [func(Interval(l,r),*args, **kwargs) for l,r in zip(self.left,self.right)] + + """ + ints = [ + func(Interval(l, r), *args, **kwargs) for l, r in zip(self.left, self.right) + ] return Pbox( - left = np.array([i.left for i in ints]), - right = np.array([i.right for i in ints]) + left=np.array([i.left for i in ints]), + right=np.array([i.right for i in ints]), ) - - def interpolate(self, steps: int, method: str = 'linear', inplace = True) -> None: - ''' + + def interpolate(self, steps: int, method: str = "linear", inplace=True) -> None: + """ Function to interpolate a p-box to a new number of steps. - + **Arguments:** ``steps`` (``int``): Number of steps to interpolate to. ``method`` (``str``): Interpolation method to use. Must be one of ``linear``, ``cubicspline`` or ``step``. ``inplace`` (``bool``): If ``True``, the p-box is interpolated in place. If ``False``, a new p-box is returned. - + .. note:: - + ``method = linear`` uses ``numpy.interp`` - ``method = cubicspline`` uses ``scipy.interpolate.CubicSpline`` + ``method = cubicspline`` uses ``scipy.interpolate.CubicSpline`` ``method = step`` uses a step interpolation method. - + .. example:: - + .. image:: https://github.com/Institute-for-Risk-and-Uncertainty/pba-for-python/blob/master/docs/images/interpolation.png?raw=true - - ''' + + """ if steps < self.steps: - raise ValueError("New number of steps must be greater than current number of steps") - + raise ValueError( + "New number of steps must be greater than current number of steps" + ) + new_left, new_right = _interpolate(self.left, self.right, steps, method) - + if inplace: self.left = new_left self.right = new_right self.steps = steps - + else: - return Pbox(left = new_left, right = new_right, steps = steps,shape = self.shape) - + return Pbox(left=new_left, right=new_right, steps=steps, shape=self.shape) + ### Arithmetic Functions - def add(self, other: Union["Pbox",Interval,float,int], method = 'f', enforce_steps = True) -> "Pbox": - ''' + def add( + self, other: Union["Pbox", Interval, float, int], method="f", enforce_steps=True + ) -> "Pbox": + """ Adds other to Pbox to other using the defined dependency method - ''' + """ try: - return _arithmetic(self, other, method, op = lambda x,y: x+y, enforce_steps = enforce_steps) + return _arithmetic( + self, other, method, op=lambda x, y: x + y, enforce_steps=enforce_steps + ) except NotImplementedError: - raise NotImplementedError(f"Addition of {other.__class__.__name__} to Pbox not implemented") + raise NotImplementedError( + f"Addition of {other.__class__.__name__} to Pbox not implemented" + ) except: raise Exception(f"Addition of {other.__class__.__name__} to Pbox failed") - def sub(self, other, method = 'f'): - ''' + def sub(self, other, method="f"): + """ Subtracts other from Pbox using the defined dependency method - ''' - if method == 'o': - method = 'p' - elif method == 'p': - method = 'o' + """ + if method == "o": + method = "p" + elif method == "p": + method = "o" return self.add(-other, method) - def mul(self, other, method = 'f'): - ''' + def mul(self, other, method="f"): + """ Multiplies other and Pbox using the defined dependency method - ''' + """ try: - return _arithmetic(self, other, method, op = lambda x,y: x*y) + return _arithmetic(self, other, method, op=lambda x, y: x * y) except NotImplementedError: - raise NotImplementedError(f"Multiplication of {other.__class__.__name__} and Pbox not implemented") + raise NotImplementedError( + f"Multiplication of {other.__class__.__name__} and Pbox not implemented" + ) except: - raise Exception(f"Multiplication of {other.__class__.__name__} to Pbox failed") + raise Exception( + f"Multiplication of {other.__class__.__name__} to Pbox failed" + ) - def div(self, other, method = 'f'): - ''' + def div(self, other, method="f"): + """ Divides Pbox by other using the defined dependency method - ''' - if method == 'o': - method = 'p' - elif method == 'p': - method = 'o' + """ + if method == "o": + method = "p" + elif method == "p": + method = "o" - if isinstance(other, (Interval,Pbox)): + if isinstance(other, (Interval, Pbox)): return self.mul(other.recip(), method) else: - return self.mul(1/other, method) + return self.mul(1 / other, method) - def pow(self, other: Union["Pbox",Interval,float,int], method = 'f') -> "Pbox": - ''' + def pow(self, other: Union["Pbox", Interval, float, int], method="f") -> "Pbox": + """ Raises a p-box to the power of other using the defined dependency method - ''' + """ try: - return _arithmetic(self, other, method, op = lambda x,y: x**y) + return _arithmetic(self, other, method, op=lambda x, y: x**y) except NotImplementedError: - raise NotImplementedError(f"Power of {other.__class__.__name__} to Pbox not implemented") + raise NotImplementedError( + f"Power of {other.__class__.__name__} to Pbox not implemented" + ) except: raise Exception(f"Power of {other.__class__.__name__} to Pbox failed") - def exp(self): - return self.unary(function = lambda x: x.exp()) - + def exp(self): + return self.unary(function=lambda x: x.exp()) + def sqrt(self): - return self.unary(function = lambda x: x.sqrt()) - + return self.unary(function=lambda x: x.sqrt()) + def recip(self): - ''' + """ Calculates the reciprocal of a p-box. - + .. error:: - + ``DivisionByZero`` is raised if the p-box contains 0. - - ''' + + """ _check_div_by_zero(self) return Pbox( - left = 1 / np.flip(self.right), - right = 1 / np.flip(self.left), - steps = self.steps + left=1 / np.flip(self.right), right=1 / np.flip(self.left), steps=self.steps ) - - def lt(self, other, method = 'f'): + def lt(self, other, method="f"): b = self.add(-other, method) - return b.get_probability(0) # return (self.add(-other, method)).get_probability(0) + return b.get_probability( + 0 + ) # return (self.add(-other, method)).get_probability(0) - def le(self, other, method = 'f'): + def le(self, other, method="f"): b = self.add(-other, method) - return b.get_probability(0) # how is the "or equal to" affecting the calculation? + return b.get_probability( + 0 + ) # how is the "or equal to" affecting the calculation? - def gt(self, other, method = 'f'): - self = - self + def gt(self, other, method="f"): + self = -self b = self.add(other, method) - return b.get_probability(0) # maybe 1-prob ? + return b.get_probability(0) # maybe 1-prob ? - def ge(self, other, method = 'f'): - self = - self + def ge(self, other, method="f"): + self = -self b = self.add(other, method) return b.get_probability(0) - def min(self, other, method='f'): + def min(self, other, method="f"): """ Returns a new Pbox object that represents the element-wise minimum of two Pboxes. **Arguments**: - + ``other``: Another Pbox, Interval or a numeric value. ``method``: Calculation method to determine the minimum. Can be one of 'f', 'p', 'o', 'i'. **Returns**: - + ``Pbox`` - + """ - if isinstance(other, (Interval,Pbox)): - if method == 'f': - return _arithmetic(self, other, method, op = lambda x,y: min(list(x)+list(y))) + if isinstance(other, (Interval, Pbox)): + if method == "f": + return _arithmetic( + self, other, method, op=lambda x, y: min(list(x) + list(y)) + ) else: - return _arithmetic(self, other, method, op = lambda x,y: np.minimum(x,y)) + return _arithmetic( + self, other, method, op=lambda x, y: np.minimum(x, y) + ) else: try: return Pbox( - left = np.array([i if i < other else other for i in self.left]), - right = np.array([i if i < other else other for i in self.right]), + left=np.array([i if i < other else other for i in self.left]), + right=np.array([i if i < other else other for i in self.right]), ) except: - return NotImplemented(f"Minimum of {other.__class__.__name__} and Pbox not implemented") + return NotImplemented( + f"Minimum of {other.__class__.__name__} and Pbox not implemented" + ) - def max(self, other, method = 'f'): + def max(self, other, method="f"): """ Returns a new Pbox object that represents the element-wise minimum of two Pboxes. **Arguments**: - + ``other``: Another Pbox, Interval or a numeric value. ``method``: Calculation method to determine the minimum. Can be one of 'f', 'p', 'o', 'i'. **Returns**: - + ``Pbox`` - + """ - if isinstance(other, (Interval,Pbox)): - if method == 'f': - return _arithmetic(self, other, method, op = lambda x,y: max(list(x)+list(y))) + if isinstance(other, (Interval, Pbox)): + if method == "f": + return _arithmetic( + self, other, method, op=lambda x, y: max(list(x) + list(y)) + ) else: - return _arithmetic(self, other, method, op = lambda x,y: np.maximum(x,y)) + return _arithmetic( + self, other, method, op=lambda x, y: np.maximum(x, y) + ) else: try: return Pbox( - left = np.array([i if i > other else other for i in self.left]), - right = np.array([i if i > other else other for i in self.right]), + left=np.array([i if i > other else other for i in self.left]), + right=np.array([i if i > other else other for i in self.right]), ) except: - return NotImplemented(f"Minimum of {other.__class__.__name__} and Pbox not implemented") + return NotImplemented( + f"Minimum of {other.__class__.__name__} and Pbox not implemented" + ) - def truncate(self, a: Union[Interval,float,int], b: Union[float,int] = None, method = 'f'): - ''' + def truncate( + self, a: Union[Interval, float, int], b: Union[float, int] = None, method="f" + ): + """ Truncates a p-box to the interval [a,b], or a if b is not specified and a is an Interval. - + **Arguments:** ``a`` (``Interval``, ``float``, ``int``): The lower bound of the truncation interval. ``b`` (``float``, ``int``): The upper bound of the truncation interval. If not specified, the upper bound of ``a`` is used. ``method`` (``str``): The dependency method to use. Can be one of ``f``, ``p``, ``o``, ``i``. - + .. admonition:: Implementation - + >>> self.min(a).max(b) - + .. error:: - + ``ValueError`` is raised if ``b`` is not specified and ``a`` is not an ``Interval``. - - ''' + + """ if isinstance(a, Interval): a, b = a.left, a.right else: if b is None: raise ValueError("b must be specified if a is not an Interval") - - return self.min(a,method=method).max(b,method=method) - - def logicaland(self, other, method = 'f'): # conjunction - if method=='i': - return(self.mul(other,method)) # independence a * b - elif method=='p': - return(self.min(other,method)) # perfect min(a, b) - elif method=='o': - return(max(self.add(other,method)-1, 0)) # opposite max(a + b – 1, 0) - elif method=='+': - return(self.min(other,method)) # positive env(a * b, min(a, b)) - elif method=='-': - return(self.min(other,method)) # negative env(max(a + b – 1, 0), a * b) + + return self.min(a, method=method).max(b, method=method) + + def logicaland(self, other, method="f"): # conjunction + if method == "i": + return self.mul(other, method) # independence a * b + elif method == "p": + return self.min(other, method) # perfect min(a, b) + elif method == "o": + return max(self.add(other, method) - 1, 0) # opposite max(a + b – 1, 0) + elif method == "+": + return self.min(other, method) # positive env(a * b, min(a, b)) + elif method == "-": + return self.min(other, method) # negative env(max(a + b – 1, 0), a * b) else: - return self.min(other,method).env(max(0, self.add(other,method) - 1)) - - def logicalor(self, other, method = 'f'): # disjunction - if method=='i': - return(1 - (1-self) * (1-other)) # independent 1 – (1 – a) * (1 – b) - elif method=='p': - return(self.max(other,method)) # perfect max(a, b) - elif method=='o': - return(min(self.add(other,method),1)) # opposite min(1, a + b) + return self.min(other, method).env(max(0, self.add(other, method) - 1)) + + def logicalor(self, other, method="f"): # disjunction + if method == "i": + return 1 - (1 - self) * (1 - other) # independent 1 – (1 – a) * (1 – b) + elif method == "p": + return self.max(other, method) # perfect max(a, b) + elif method == "o": + return min(self.add(other, method), 1) # opposite min(1, a + b) # elif method=='+': # return(env(,min(self.add(other,method),1)) # positive env(max(a, b), 1 – (1 – a) * (1 – b)) # elif method=='-': # return() # negative env(1 – (1 – a) * (1 – b), min(1, a + b)) else: - return(self.env(self.max(other,method), min(self.add(other,method),1))) + return self.env(self.max(other, method), min(self.add(other, method), 1)) def env(self, other): """ .. _pbox.env: - + Computes the envelope of two Pboxes. **Arguments**: - + ``other``: Another Pbox, Interval or a numeric value. **Returns**: @@ -764,32 +819,39 @@ def env(self, other): ``Pbox`` .. error:: - + ``NotImplementedError`` is raised if ``other`` is not a Pbox. Imputation of other needs to be done manually - + """ - if other.__class__.__name__ == 'Pbox': + if other.__class__.__name__ == "Pbox": if self.steps != other.steps: raise ArithmeticError("Both Pboxes must have the same number of steps") else: - other = Pbox(other, steps = self.steps) - - nleft = np.minimum(self.left, other.left) - nright = np.maximum(self.right, other.right) + other = Pbox(other, steps=self.steps) - return Pbox( - left = nleft, - right = nright, - steps = self.steps - ) + nleft = np.minimum(self.left, other.left) + nright = np.maximum(self.right, other.right) - def show(self,figax = None, now = True, title = '', xlabel = 'x', ylabel = r'$\Pr(x \leq X)$',left_col = 'red',right_col = 'black', label=None, **kwargs): - ''' + return Pbox(left=nleft, right=nright, steps=self.steps) + + def show( + self, + figax=None, + now=True, + title="", + xlabel="x", + ylabel=r"$\Pr(x \leq X)$", + left_col="red", + right_col="black", + label=None, + **kwargs, + ): + """ Plots the p-box - + **Arguments:** - + ``figax`` (``tuple``): Tuple containing a matplotlib figure and axis object. If not specified, a new figure and axis object are created. ``now`` (``bool``): If ``True``, the figure is shown. If ``False``, the figure is returned. ``title`` (``str``): Title of the plot. @@ -799,21 +861,21 @@ def show(self,figax = None, now = True, title = '', xlabel = 'x', ylabel = r'$\P ``left_col`` (``str``): Colour of the left bound of the p-box. ``right_col`` (``str``): Colour of the right bound of the p-box. ``kwargs`` (``dict``): Additional keyword arguments to pass to ``matplotlib.pyplot.plot``. - + **Example:** .. code-block:: python - + >>> p = pba.N([-1,1],1) >>> fig, ax = plt.subplots() >>> p.show(figax = (fig,ax), now = True, title = 'Example', xlabel = 'x', ylabel = r'$\Pr(x \leq X)$',left_col = 'red',right_col = 'black') - ''' + """ if figax is None: fig, ax = plt.subplots() else: fig, ax = figax - + # now respects discretization L = self.left R = self.right @@ -821,24 +883,35 @@ def show(self,figax = None, now = True, title = '', xlabel = 'x', ylabel = r'$\P LL = np.concatenate((L, L, np.array([R[-1]]))) RR = np.concatenate((np.array([L[0]]), R, R)) - ii = np.concatenate((np.arange(steps), np.arange(1, steps + 1), np.array([steps]))) / steps - jj = np.concatenate((np.array([0]),np.arange(steps + 1), np.arange(1, steps))) / steps + ii = ( + np.concatenate( + (np.arange(steps), np.arange(1, steps + 1), np.array([steps])) + ) + / steps + ) + jj = ( + np.concatenate((np.array([0]), np.arange(steps + 1), np.arange(1, steps))) + / steps + ) - ii.sort(); jj.sort(); LL.sort(); RR.sort() + ii.sort() + jj.sort() + LL.sort() + RR.sort() - if 'color' in kwargs.keys(): - ax.plot(LL,ii,label = label, **kwargs) - ax.plot(RR,jj,**kwargs) + if "color" in kwargs.keys(): + ax.plot(LL, ii, label=label, **kwargs) + ax.plot(RR, jj, **kwargs) else: - ax.plot(LL,ii,color=left_col,label = label, **kwargs) - ax.plot(RR,jj,color=right_col,**kwargs) - - if title != '' : - ax.set_title(title,**kwargs) - + ax.plot(LL, ii, color=left_col, label=label, **kwargs) + ax.plot(RR, jj, color=right_col, **kwargs) + + if title != "": + ax.set_title(title, **kwargs) + ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) - + if now: fig.show() else: @@ -852,10 +925,10 @@ def get_interval(self, *args) -> Interval: if args[0] == 1: # asking for whole pbox bounds - return Interval(min(self.left),max(self.right)) + return Interval(min(self.left), max(self.right)) - p1 = (1-args[0])/2 - p2 = 1-p1 + p1 = (1 - args[0]) / 2 + p2 = 1 - p1 elif len(args) == 2: @@ -863,34 +936,33 @@ def get_interval(self, *args) -> Interval: p2 = args[1] else: - raise Exception('Too many inputs') + raise Exception("Too many inputs") - y = np.append(np.insert(np.linspace(0,1,self.steps),0,0),1) + y = np.append(np.insert(np.linspace(0, 1, self.steps), 0, 0), 1) y1 = 0 while y[y1] < p1: y1 += 1 - y2 = len(y)-1 + y2 = len(y) - 1 while y[y2] > p2: y2 -= 1 x1 = self.left[y1] x2 = self.right[y2] - return Interval(x1,x2) + return Interval(x1, x2) def get_probability(self, val) -> Interval: - ''' + """ Returns the interval - ''' - - p = np.append(np.insert(np.linspace(0,1,self.steps),0,0),1) + """ + + p = np.append(np.insert(np.linspace(0, 1, self.steps), 0, 0), 1) i = 0 while i < self.steps and self.left[i] < val: i += 1 - ub = p[i] j = 0 @@ -898,148 +970,152 @@ def get_probability(self, val) -> Interval: while j < self.steps and self.right[j] < val: j += 1 - lb = p[j] - return Interval(lb,ub) + return Interval(lb, ub) def summary(self) -> str: - ''' + """ Returns a summary of the p-box - ''' - s = 'Pbox Summary\n' - s += '------------\n' - if self.shape != '': s += f'Shape: {self.shape}\n' - s += f'Range: {self.support()}\n' - s += f'Mean: {self.mean}\n' - s += f'Variance: {self.var}\n' - s += f'Steps: {self.steps}\n' + """ + s = "Pbox Summary\n" + s += "------------\n" + if self.shape != "": + s += f"Shape: {self.shape}\n" + s += f"Range: {self.support()}\n" + s += f"Mean: {self.mean}\n" + s += f"Variance: {self.var}\n" + s += f"Steps: {self.steps}\n" return s def support(self) -> Interval: - ''' + """ Returns the range of the pbox - + .. admonition:: Implementation - + >>> Interval(self.lo(),self.hi()) - - ''' - return Interval(min(self.left),max(self.right)) + + """ + return Interval(min(self.left), max(self.right)) def mean(self) -> Interval: - ''' + """ Returns the mean of the pbox - ''' + """ return self.mean def median(self) -> Interval: - ''' + """ Returns the median of the distribution - ''' - return Interval(np.median(self.left),np.median(self.right)) - + """ + return Interval(np.median(self.left), np.median(self.right)) - def straddles(self,N, endpoints = True) -> bool: + def straddles(self, N, endpoints=True) -> bool: """ Checks whether a number is within the p-box's support - + **Arguments:** - + ``N`` (``float``): Number to check ``endpoints`` (``bool``): If ``True``, the endpoints of the p-box are included in the check. - + **Returns:** - + ``bool`` """ - return self.support().straddles(N,endpoints) + return self.support().straddles(N, endpoints) - def straddles_zero(self,endpoints = True) -> bool: + def straddles_zero(self, endpoints=True) -> bool: """ Checks whether :math:`0` is within the p-box """ - return self.straddles(0,endpoints) + return self.straddles(0, endpoints) def imp(self, other): - ''' + """ Returns the imposition of self with other - ''' - if other.__class__.__name__ != 'Pbox': + """ + if other.__class__.__name__ != "Pbox": try: pbox = Pbox(pbox) except: - raise TypeError("Unable to convert %s object (%s) to Pbox" %(type(pbox),pbox)) - + raise TypeError( + "Unable to convert %s object (%s) to Pbox" % (type(pbox), pbox) + ) + u = [] d = [] - + assert self.steps == other.steps - - for sL,sR,oL,oR in zip(self.left,self.right,other.left,other.right): - - if max(sL,oL) > min(sR,oR): + + for sL, sR, oL, oR in zip(self.left, self.right, other.left, other.right): + + if max(sL, oL) > min(sR, oR): raise Exception("Imposition does not exist") - - u.append(max(sL,oL)) - d.append(min(sR,oR)) - - return Pbox( - left = u, - right = d - ) + + u.append(max(sL, oL)) + d.append(min(sR, oR)) + + return Pbox(left=u, right=d) def left_list(implist, verbose=False): - if not hasattr(implist,"__iter__"): + if not hasattr(implist, "__iter__"): return np.array(implist) return np.array([imp.lo() for imp in implist]) + def right_list(implist, verbose=False): - if not hasattr(implist,"__iter__"): + if not hasattr(implist, "__iter__"): return np.array(implist) return np.array([imp.hi() for imp in implist]) - -def truncate(pbox,min,max): - return pbox.truncate(min,max) -def imposition(*args: Union[Pbox,Interval,float,int]): - ''' + +def truncate(pbox, min, max): + return pbox.truncate(min, max) + + +def imposition(*args: Union[Pbox, Interval, float, int]): + """ Returns the imposition of the p-boxes in *args - + Parameters ---------- *args : Number of p-boxes or objects to be mixed - + Returns ---------- - Pbox - ''' + Pbox + """ x = [] for pbox in args: - if pbox.__class__.__name__ != 'Pbox': + if pbox.__class__.__name__ != "Pbox": try: pbox = Pbox(pbox) except: - raise TypeError("Unable to convert %s object (%s) to Pbox" %(type(pbox),pbox)) + raise TypeError( + "Unable to convert %s object (%s) to Pbox" % (type(pbox), pbox) + ) x.append(pbox) p = x[0] - - for i in range(1,len(x)): + + for i in range(1, len(x)): p.imp(x[i]) - + return p - + + def mixture( - *args: Union[Pbox,Interval,float,int], - weights: List[Union[float,int]] = [], - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + *args: Union[Pbox, Interval, float, int], + weights: List[Union[float, int]] = [], + steps: int = Pbox.STEPS, +) -> Pbox: + """ Mixes the pboxes in *args Parameters ---------- @@ -1047,36 +1123,36 @@ def mixture( Number of p-boxes or objects to be mixed weights: Right side of box - + Returns ---------- Pbox - ''' - #TODO: IMPROVE READBILITY + """ + # TODO: IMPROVE READBILITY x = [] for pbox in args: - if pbox.__class__.__name__ != 'Pbox': + if pbox.__class__.__name__ != "Pbox": try: pbox = Pbox(pbox) except: - raise TypeError("Unable to convert %s object (%s) to Pbox" %(type(pbox),pbox)) + raise TypeError( + "Unable to convert %s object (%s) to Pbox" % (type(pbox), pbox) + ) x.append(pbox) k = len(x) if weights == []: weights = [1] * k - # temporary hack # k = 2 # x = [self, x] # w = [1,1] - if k != len(weights): - return('Need same number of weights as arguments for mixture') - weights = [i/sum(weights) for i in weights] # w = w / sum(w) + return "Need same number of weights as arguments for mixture" + weights = [i / sum(weights) for i in weights] # w = w / sum(w) u = [] d = [] n = [] @@ -1086,10 +1162,12 @@ def mixture( vl = [] vh = [] v = [] - for i in range(k) : + for i in range(k): u = u + list(x[i].left) - d = np.append(d,x[i].right) - n = n + [weights[i] / x[i].steps] * x[i].steps # w[i]*rep(1/x[i].steps,x[i].steps)) + d = np.append(d, x[i].right) + n = ( + n + [weights[i] / x[i].steps] * x[i].steps + ) # w[i]*rep(1/x[i].steps,x[i].steps)) # mu = mean(x[i]) # ml = ml + [mu.left()] @@ -1104,99 +1182,137 @@ def mixture( MR = x[i].mean.right VL = x[i].var.left VR = x[i].var.right - m = m + [Interval(ML,MR)] - v = v + [Interval(VL,VR)] + m = m + [Interval(ML, MR)] + v = v + [Interval(VL, VR)] ml = ml + [ML] mh = mh + [MR] vl = vl + [VL] vh = vh + [VR] - n = [_/sum(n) for _ in n] # n = n / sum(n) + n = [_ / sum(n) for _ in n] # n = n / sum(n) su = sorted(u) su = [su[0]] + su - pu = [0] + list(np.cumsum([n[i] for i in np.argsort(u)])) # pu = c(0,cumsum(n[order(u)])) - sd = sorted(d); sd = sd + [sd[-1]] - pd = list(np.cumsum([n[i] for i in np.argsort(d)])) + [1] # pd = c(cumsum(n[order(d)]),1) - u = []; d = [] + pu = [0] + list( + np.cumsum([n[i] for i in np.argsort(u)]) + ) # pu = c(0,cumsum(n[order(u)])) + sd = sorted(d) + sd = sd + [sd[-1]] + pd = list(np.cumsum([n[i] for i in np.argsort(d)])) + [ + 1 + ] # pd = c(cumsum(n[order(d)]),1) + u = [] + d = [] j = len(pu) - 1 - for p in reversed(np.arange(steps)/steps) : # ii = np.arange(steps))/steps # ii = 0: (Pbox$steps-1) / Pbox$steps - while p < pu[j] : j = j - 1 # repeat {if (pu[j] <= p) break; j = j - 1} + for p in reversed( + np.arange(steps) / steps + ): # ii = np.arange(steps))/steps # ii = 0: (Pbox$steps-1) / Pbox$steps + while p < pu[j]: + j = j - 1 # repeat {if (pu[j] <= p) break; j = j - 1} u = [su[j]] + u j = 0 - for p in (np.arange(steps)+1)/steps : # jj = (np.arange(steps)+1)/steps # jj = 1: Pbox$steps / Pbox$steps - while pd[j] < p : j = j + 1 # repeat {if (p <= pu[j]) break; j = j + 1} + for p in ( + np.arange(steps) + 1 + ) / steps: # jj = (np.arange(steps)+1)/steps # jj = 1: Pbox$steps / Pbox$steps + while pd[j] < p: + j = j + 1 # repeat {if (p <= pu[j]) break; j = j + 1} d = d + [sd[j]] - mu = Interval(np.sum([W * M for M,W in zip(weights,ml)]), np.sum([W * M for M,W in zip(weights,mh)])) + mu = Interval( + np.sum([W * M for M, W in zip(weights, ml)]), + np.sum([W * M for M, W in zip(weights, mh)]), + ) s2 = 0 - for i in range(k) : s2 = s2 + weights[i] * (v[i] + m[i]**2) + for i in range(k): + s2 = s2 + weights[i] * (v[i] + m[i] ** 2) s2 = s2 - mu**2 - return Pbox(np.array(u),np.array(d), mean = Interval(mu.left, mu.right), var = Interval(s2.left, s2.right), steps = steps, check_moments=False) + return Pbox( + np.array(u), + np.array(d), + mean=Interval(mu.left, mu.right), + var=Interval(s2.left, s2.right), + steps=steps, + check_moments=False, + ) + def change_default_arithmetic_method(method): - ''' + """ Changes the default arithmetic method for p-boxes - + **Arguments:** ``method`` (``str``): Method to use. Must be one of ``f``, ``p``, ``o``, ``i``. - - ''' - if method not in ['f','p','o','i']: + """ + + if method not in ["f", "p", "o", "i"]: raise ValueError("Method must be one of 'f', 'p', 'o', 'i'") - - def nadd(self, other: Union["Pbox",Interval,float,int], method = method, enforce_steps = True) -> "Pbox": + + def nadd( + self, + other: Union["Pbox", Interval, float, int], + method=method, + enforce_steps=True, + ) -> "Pbox": try: - return _arithmetic(self, other, method, op = lambda x,y: x+y, enforce_steps = enforce_steps) + return _arithmetic( + self, other, method, op=lambda x, y: x + y, enforce_steps=enforce_steps + ) except NotImplementedError: - raise NotImplementedError(f"Addition of {other.__class__.__name__} to Pbox not implemented") + raise NotImplementedError( + f"Addition of {other.__class__.__name__} to Pbox not implemented" + ) except: raise Exception(f"Addition of {other.__class__.__name__} to Pbox failed") - def nsub(self, other, method = method): + def nsub(self, other, method=method): - if method == 'o': - method = 'p' - elif method == 'p': - method = 'o' + if method == "o": + method = "p" + elif method == "p": + method = "o" return self.add(-other, method) - def nmul(self, other, method = method): + def nmul(self, other, method=method): try: - return _arithmetic(self, other, method, op = lambda x,y: x*y) + return _arithmetic(self, other, method, op=lambda x, y: x * y) except NotImplementedError: - raise NotImplementedError(f"Multiplication of {other.__class__.__name__} and Pbox not implemented") + raise NotImplementedError( + f"Multiplication of {other.__class__.__name__} and Pbox not implemented" + ) except: - raise Exception(f"Multiplication of {other.__class__.__name__} to Pbox failed") + raise Exception( + f"Multiplication of {other.__class__.__name__} to Pbox failed" + ) - def ndiv(self, other, method = method): + def ndiv(self, other, method=method): - if method == 'o': - method = 'p' - elif method == 'p': - method = 'o' + if method == "o": + method = "p" + elif method == "p": + method = "o" - if isinstance(other, (Interval,Pbox)): + if isinstance(other, (Interval, Pbox)): return self.mul(other.recip(), method) else: - return self.mul(1/other, method) + return self.mul(1 / other, method) - def npow(self, other: Union["Pbox",Interval,float,int], method = method) -> "Pbox": + def npow(self, other: Union["Pbox", Interval, float, int], method=method) -> "Pbox": try: - return _arithmetic(self, other, method, op = lambda x,y: x**y) + return _arithmetic(self, other, method, op=lambda x, y: x**y) except NotImplementedError: - raise NotImplementedError(f"Power of {other.__class__.__name__} to Pbox not implemented") + raise NotImplementedError( + f"Power of {other.__class__.__name__} to Pbox not implemented" + ) except: raise Exception(f"Power of {other.__class__.__name__} to Pbox failed") - + Pbox.__add__ = nadd Pbox.__sub__ = nsub Pbox.__mul__ = nmul Pbox.__truediv__ = ndiv Pbox.__pow__ = npow - \ No newline at end of file diff --git a/pba/pbox_constructors/distributions.py b/pba/pbox_constructors/distributions.py index 23697c4..3bbffa5 100644 --- a/pba/pbox_constructors/distributions.py +++ b/pba/pbox_constructors/distributions.py @@ -1,21 +1,21 @@ -''' +""" This module contains functions to create p-boxes for various distributions. These functions have been hand-coded to ensure accuracy. All other distributions can be found in the :mod:`pba.pbox_constructors.parametric` module. -''' +""" __all__ = [ - 'beta', - 'foldnorm', - 'normal', - 'N', - 'unif', - 'U', - 'KN', - 'KM', - 'lognormal', - 'lognorm', - 'trapz', - 'uniform', - 'weibull' + "beta", + "foldnorm", + "normal", + "N", + "unif", + "U", + "KN", + "KM", + "lognormal", + "lognorm", + "trapz", + "uniform", + "weibull", ] @@ -32,106 +32,103 @@ from .parametric import parametric, norm, weibull_max, weibull_min -def lognormal(mean, var, steps = 200): - ''' +def lognormal(mean, var, steps=200): + """ Creates a p-box for the lognormal distribution - + *Note: the parameters used are the mean and variance of the lognormal distribution not the mean and variance of the underlying normal* See: `[1]` `[2]` - - + + Parameters ---------- - mean : + mean : mean of the lognormal distribution var : variance of the lognormal distribution - + Returns ---------- Pbox - - ''' + + """ if steps > 1000: - x = np.linspace(1/steps,1-1/steps,steps) + x = np.linspace(1 / steps, 1 - 1 / steps, steps) else: - x = np.linspace(0.001,0.999,steps) + x = np.linspace(0.001, 0.999, steps) + + if mean.__class__.__name__ != "Interval": + mean = Interval(mean, mean) + if var.__class__.__name__ != "Interval": + var = Interval(var, var) - if mean.__class__.__name__ != 'Interval': - mean = Interval(mean,mean) - if var.__class__.__name__ != 'Interval': - var = Interval(var,var) - def __lognorm(mean, var): - - sigma = np.sqrt(np.log1p(var/mean**2)) - mu = np.log(mean) - 0.5*sigma*sigma - - return stats.lognorm(sigma, loc = 0,scale = np.exp(mu)) - - - bounds = np.array([ - __lognorm(mean.left,var.left).ppf(x), - __lognorm(mean.right,var.left).ppf(x), - __lognorm(mean.left,var.right).ppf(x), - __lognorm(mean.right,var.right).ppf(x) - ]) + + sigma = np.sqrt(np.log1p(var / mean**2)) + mu = np.log(mean) - 0.5 * sigma * sigma + + return stats.lognorm(sigma, loc=0, scale=np.exp(mu)) + + bounds = np.array( + [ + __lognorm(mean.left, var.left).ppf(x), + __lognorm(mean.right, var.left).ppf(x), + __lognorm(mean.left, var.right).ppf(x), + __lognorm(mean.right, var.right).ppf(x), + ] + ) print(bounds) - Left = np.min(bounds,axis = 0) - Right = np.max(bounds,axis = 0) + Left = np.min(bounds, axis=0) + Right = np.max(bounds, axis=0) Left = np.array(Left) Right = np.array(Right) - return Pbox( - Left, - Right, - steps = steps, - shape='lognormal' - ) - - + return Pbox(Left, Right, steps=steps, shape="lognormal") + + lognorm = lognormal -def beta(a,b, steps = 200): - ''' - Beta distribution. - ''' + +def beta(a, b, steps=200): + """ + Beta distribution. + """ args = list(args) - if not isinstance(a,Interval): + if not isinstance(a, Interval): a = Interval(a) - if not isinstance(b,Interval): + if not isinstance(b, Interval): b = Interval(b) - + if sometimes(a < 0): - raise ValueError('a must be always greater than 0') + raise ValueError("a must be always greater than 0") if sometimes(b < 0): - raise ValueError('b must be always greater than 0') - + raise ValueError("b must be always greater than 0") + if a.left == 0: a.left = 1e-5 if b.left == 0: b.left = 1e-5 - - return parametric('beta').make_pbox(a,b,steps = steps, support=Interval(0,1)) -def foldnorm(mu,s, steps = 200): + return parametric("beta").make_pbox(a, b, steps=steps, support=Interval(0, 1)) + - x = np.linspace(0.0001,0.9999,steps) - if mu.__class__.__name__ != 'Interval': +def foldnorm(mu, s, steps=200): + + x = np.linspace(0.0001, 0.9999, steps) + if mu.__class__.__name__ != "Interval": mu = Interval(mu) - if s.__class__.__name__ != 'Interval': + if s.__class__.__name__ != "Interval": s = Interval(s) new_args = [ - [mu.lo()/s.lo(),0,s.lo()], - [mu.hi()/s.lo(),0,s.lo()], - [mu.lo()/s.hi(),0,s.hi()], - [mu.hi()/s.hi(),0,s.hi()] + [mu.lo() / s.lo(), 0, s.lo()], + [mu.hi() / s.lo(), 0, s.lo()], + [mu.lo() / s.hi(), 0, s.hi()], + [mu.hi() / s.hi(), 0, s.hi()], ] - bounds = [] mean_hi = -np.inf @@ -141,8 +138,8 @@ def foldnorm(mu,s, steps = 200): for a in new_args: - bounds.append(stats.foldnorm.ppf(x,*a)) - bmean, bvar = stats.foldnorm.stats(*a, moments = 'mv') + bounds.append(stats.foldnorm.ppf(x, *a)) + bmean, bvar = stats.foldnorm.stats(*a, moments="mv") if bmean < mean_lo: mean_lo = bmean @@ -153,26 +150,26 @@ def foldnorm(mu,s, steps = 200): if bvar < var_lo: var_lo = bvar - Left = [min([b[i] for b in bounds]) for i in range(steps)] Right = [max([b[i] for b in bounds]) for i in range(steps)] - var = Interval(np.float64(var_lo),np.float64(var_hi)) - mean = Interval(np.float64(mean_lo),np.float64(mean_hi)) + var = Interval(np.float64(var_lo), np.float64(var_hi)) + mean = Interval(np.float64(mean_lo), np.float64(mean_hi)) Left = np.array(Left) Right = np.array(Right) return Pbox( - Left, - Right, - steps = steps, - shape = 'foldnorm', - mean_left = mean.left, - mean_right = mean.right, - var_left = var.left, - var_right = var.right - ) + Left, + Right, + steps=steps, + shape="foldnorm", + mean_left=mean.left, + mean_right=mean.right, + var_left=var.left, + var_right=var.right, + ) + # def frechet_r(*args, steps = 200): # args = list(args) @@ -212,78 +209,76 @@ def foldnorm(mu,s, steps = 200): # var_right = var.right # ) -def trapz(a,b,c,d , steps = 200): - if a.__class__.__name__ != 'Interval': + +def trapz(a, b, c, d, steps=200): + if a.__class__.__name__ != "Interval": a = Interval(a) - if b.__class__.__name__ != 'Interval': + if b.__class__.__name__ != "Interval": b = Interval(b) - if c.__class__.__name__ != 'Interval': + if c.__class__.__name__ != "Interval": c = Interval(c) - if d.__class__.__name__ != 'Interval': + if d.__class__.__name__ != "Interval": d = Interval(d) - x = np.linspace(0.0001,0.9999,steps) - left = stats.trapz.ppf(x,*sorted([b.lo()/d.lo(),c.lo()/d.lo(),a.lo(),d.lo()-a.lo()])) - right = stats.trapz.ppf(x,*sorted([b.hi()/d.hi(),c.hi()/d.hi(),a.hi(),d.hi()-a.hi()])) + x = np.linspace(0.0001, 0.9999, steps) + left = stats.trapz.ppf( + x, *sorted([b.lo() / d.lo(), c.lo() / d.lo(), a.lo(), d.lo() - a.lo()]) + ) + right = stats.trapz.ppf( + x, *sorted([b.hi() / d.hi(), c.hi() / d.hi(), a.hi(), d.hi() - a.hi()]) + ) - return Pbox( - left, - right, - steps = steps, - shape = 'trapz' - ) + return Pbox(left, right, steps=steps, shape="trapz") -def uniform(a, b, steps = 200): - x = np.linspace(0,1,steps) +def uniform(a, b, steps=200): - if a.__class__.__name__ != 'Interval': - a = Interval(a,a) - if b.__class__.__name__ != 'Interval': - b = Interval(b,b) + x = np.linspace(0, 1, steps) - Left = np.linspace(a.left,b.left) - Right = np.linspace(a.right,b.right) + if a.__class__.__name__ != "Interval": + a = Interval(a, a) + if b.__class__.__name__ != "Interval": + b = Interval(b, b) - mean = 0.5 * (a+b) - var = ((b-a)**2 )/12 + Left = np.linspace(a.left, b.left) + Right = np.linspace(a.right, b.right) + + mean = 0.5 * (a + b) + var = ((b - a) ** 2) / 12 return Pbox( - Left, - Right, - steps = steps, - shape = 'uniform', - mean_left = mean.left, - mean_right = mean.right, - var_left = var.left, - var_right = var.right - ) - -def weibull(*args, steps = 200): - - + Left, + Right, + steps=steps, + shape="uniform", + mean_left=mean.left, + mean_right=mean.right, + var_left=var.left, + var_right=var.right, + ) + + +def weibull(*args, steps=200): + wm = weibull_max(*args) wl = weibull_min(*args) - - return Pbox( - left = wl.left, - right = wm.right - ) + return Pbox(left=wl.left, right=wm.right) ### Other distributions -def KM(k,m,steps = 200): +def KM(k, m, steps=200): with catch_warnings(): simplefilter("ignore") - return beta(Interval(k,k+1),Interval(m,m+1),steps = steps) + return beta(Interval(k, k + 1), Interval(m, m + 1), steps=steps) + + +def KN(k, n, steps=200): + return KM(k, n - k, steps=steps) -def KN(k,n,steps = 200): - return KM(k,n-k,steps=steps) ### Alternate names normal = norm N = normal unif = uniform U = uniform - diff --git a/pba/pbox_constructors/non_parametric.py b/pba/pbox_constructors/non_parametric.py index 6b7ee20..766c33c 100644 --- a/pba/pbox_constructors/non_parametric.py +++ b/pba/pbox_constructors/non_parametric.py @@ -1,22 +1,23 @@ -''' +""" Functions that can be used to generate non-parametric p-boxes. These functions are used to generate p-boxes based upon the minimum, maximum, mean, median, mode, standard deviation, variance, and coefficient of variation of the variable. -''' +""" + __all__ = [ - 'what_I_know', - 'box', - 'min_max', - 'min_max_mean', - 'min_max_mean_std', - 'min_max_mean_var', - 'min_max_mode', - 'min_max_median', - 'min_max_median_is_mode', - 'mean_std', - 'mean_var', - 'pos_mean_std', - 'symmetric_mean_std', - 'from_percentiles' - ] + "what_I_know", + "box", + "min_max", + "min_max_mean", + "min_max_mean_std", + "min_max_mean_var", + "min_max_mode", + "min_max_median", + "min_max_median_is_mode", + "mean_std", + "mean_var", + "pos_mean_std", + "symmetric_mean_std", + "from_percentiles", +] from ..pbox import Pbox, imposition, NotIncreasingError from ..interval import Interval @@ -29,28 +30,29 @@ from warnings import warn import numpy as np + def what_I_know( - minimum: Optional[Union[Interval,float,int]] = None, - maximum: Optional[Union[Interval,float,int]] = None, - mean: Optional[Union[Interval,float,int]] = None, - median: Optional[Union[Interval,float,int]] = None, - mode: Optional[Union[Interval,float,int]] = None, - std: Optional[Union[Interval,float,int]] = None, - var: Optional[Union[Interval,float,int]] = None, - cv: Optional[Union[Interval,float,int]] = None, - percentiles: Optional[dict[Union[Interval,float,int]]] = None, - # coverages: Optional[Union[Interval,float,int]] = None, - # shape: Optional[Literal['unimodal', 'symmetric', 'positive', 'nonnegative', 'concave', 'convex', 'increasinghazard', 'decreasinghazard', 'discrete', 'integervalued', 'continuous', '', 'normal', 'lognormal']] = None, - # data: Optional[list] = None, - # confidence: Optional[float] = 0.95, - debug: bool = False, - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + minimum: Optional[Union[Interval, float, int]] = None, + maximum: Optional[Union[Interval, float, int]] = None, + mean: Optional[Union[Interval, float, int]] = None, + median: Optional[Union[Interval, float, int]] = None, + mode: Optional[Union[Interval, float, int]] = None, + std: Optional[Union[Interval, float, int]] = None, + var: Optional[Union[Interval, float, int]] = None, + cv: Optional[Union[Interval, float, int]] = None, + percentiles: Optional[dict[Union[Interval, float, int]]] = None, + # coverages: Optional[Union[Interval,float,int]] = None, + # shape: Optional[Literal['unimodal', 'symmetric', 'positive', 'nonnegative', 'concave', 'convex', 'increasinghazard', 'decreasinghazard', 'discrete', 'integervalued', 'continuous', '', 'normal', 'lognormal']] = None, + # data: Optional[list] = None, + # confidence: Optional[float] = 0.95, + debug: bool = False, + steps: int = Pbox.STEPS, +) -> Pbox: + """ Generates a distribution free p-box based upon the information given. This function works by calculating every possible non-parametric p-box that can be generated using the information provided. The returned p-box is the intersection of these p-boxes. - + **Parameters**: - + ``minimum``: Minimum value of the variable ``maximum``: Maximum value of the variable ``mean``: Mean value of the variable @@ -61,34 +63,36 @@ def what_I_know( ``cv``: Coefficient of variation of the variable ``percentiles``: Dictionary of percentiles and their values (e.g. {0.1: 1, 0.5: 2, 0.9: Interval(3,4)}) ``steps``: Number of steps to use in the p-box - + .. error:: - + ``ValueError``: If any of the arguments are not consistent with each other. (i.e. if ``std`` and ``var`` are both given, but ``std != sqrt(var)``) - + **Returns**: - + ``Pbox``: Imposition of possible p-boxes - - ''' - - def _print_debug(skk): print("\033[93m {}\033[00m" .format(skk),end=' ') - - def _get_pbox(func,*args,steps = steps,debug=False): - if debug: _print_debug(func.__name__) - try: - return func(*args,steps=steps) + + """ + + def _print_debug(skk): + print("\033[93m {}\033[00m".format(skk), end=" ") + + def _get_pbox(func, *args, steps=steps, debug=False): + if debug: + _print_debug(func.__name__) + try: + return func(*args, steps=steps) except: - raise Exception(f'Unable to generate {func.__name__} pbox') - + raise Exception(f"Unable to generate {func.__name__} pbox") + # if 'positive' in shape: # if minimum is None: # minimum = 0 # else: # minimum = max(0,minimum) - + # if debug: _print_debug("Shape is positive") - + # if 'negative' in shape: # if maximum is None: # maximum = 0 @@ -96,88 +100,114 @@ def _get_pbox(func,*args,steps = steps,debug=False): # maximum = min(0,maximum) # if debug: _print_debug("Shape is negative") - + if std is not None and var is not None: if std != np.sqrt(var): raise ValueError("std and var are not consistent") - + imp = [] - - if minimum is not None and maximum is not None: imp += _get_pbox(min_max,minimum,maximum,debug=debug) - - if minimum is not None and mean is not None: imp += _get_pbox(min_mean,minimum,mean,debug=debug) - - if maximum is not None and mean is not None: imp += _get_pbox(max_mean,maximum,mean,debug=debug) - - if minimum is not None and maximum is not None and mean is not None: imp += _get_pbox(min_max_mean,minimum,maximum,mean,debug=debug) - - if minimum is not None and maximum is not None and mode is not None: imp += _get_pbox(min_max_mode,minimum,maximum,mode,debug=debug) - - if minimum is not None and maximum is not None and median is not None: imp += _get_pbox(min_max_median,minimum,maximum,median,debug=debug) - - if minimum is not None and mean is not None and std is not None: imp += minimum+_get_pbox(pos_mean_std,mean-minimum,std,debug=debug) - - if maximum is not None and mean is not None and std is not None: imp += _get_pbox(pos_mean_std,maximum-mean,std,debug=debug) - maximum - - if minimum is not None and maximum is not None and mean is not None and std is not None: imp += _get_pbox(min_max_mean_std,minimum,maximum,mean,std,debug=debug) - - if minimum is not None and maximum is not None and mean is not None and var is not None: imp += _get_pbox(min_max_mean_var,minimum,maximum,mean,var,debug=debug) - - if mean is not None and std is not None: imp += _get_pbox(mean_std,mean,std,debug=debug) - - if mean is not None and var is not None: imp += _get_pbox(mean_var,mean,var,debug=debug) - - if mean is not None and cv is not None: imp += _get_pbox(mean_std,mean,cv*mean,debug=debug) + + if minimum is not None and maximum is not None: + imp += _get_pbox(min_max, minimum, maximum, debug=debug) + + if minimum is not None and mean is not None: + imp += _get_pbox(min_mean, minimum, mean, debug=debug) + + if maximum is not None and mean is not None: + imp += _get_pbox(max_mean, maximum, mean, debug=debug) + + if minimum is not None and maximum is not None and mean is not None: + imp += _get_pbox(min_max_mean, minimum, maximum, mean, debug=debug) + + if minimum is not None and maximum is not None and mode is not None: + imp += _get_pbox(min_max_mode, minimum, maximum, mode, debug=debug) + + if minimum is not None and maximum is not None and median is not None: + imp += _get_pbox(min_max_median, minimum, maximum, median, debug=debug) + + if minimum is not None and mean is not None and std is not None: + imp += minimum + _get_pbox(pos_mean_std, mean - minimum, std, debug=debug) + + if maximum is not None and mean is not None and std is not None: + imp += _get_pbox(pos_mean_std, maximum - mean, std, debug=debug) - maximum + + if ( + minimum is not None + and maximum is not None + and mean is not None + and std is not None + ): + imp += _get_pbox(min_max_mean_std, minimum, maximum, mean, std, debug=debug) + + if ( + minimum is not None + and maximum is not None + and mean is not None + and var is not None + ): + imp += _get_pbox(min_max_mean_var, minimum, maximum, mean, var, debug=debug) + + if mean is not None and std is not None: + imp += _get_pbox(mean_std, mean, std, debug=debug) + + if mean is not None and var is not None: + imp += _get_pbox(mean_var, mean, var, debug=debug) + + if mean is not None and cv is not None: + imp += _get_pbox(mean_std, mean, cv * mean, debug=debug) if len(imp) == 0: raise Exception("No valid p-boxes found") return imposition(imp) - + + def box( - a: Union[Interval,float,int], - b: Union[Interval,float,int] = None, - steps = Pbox.STEPS, - shape = 'box' - ) -> Pbox: - ''' + a: Union[Interval, float, int], + b: Union[Interval, float, int] = None, + steps=Pbox.STEPS, + shape="box", +) -> Pbox: + """ Returns a box shaped Pbox. This is equivalent to an Interval expressed as a Pbox. - + **Parameters**: ``a`` : Left side of box ``b``: Right side of box - - + + **Returns**: ``Pbox`` - - ''' + + """ if b == None: b = a - i = Interval(a,b) + i = Interval(a, b) return Pbox( - left = np.repeat(i.left,steps), - right = np.repeat(i.right,steps), - mean_left = i.left, - mean_right= i.right, - var_left = 0, - var_right=((i.right-i.left)**2)/4, - steps = steps, - shape=shape + left=np.repeat(i.left, steps), + right=np.repeat(i.right, steps), + mean_left=i.left, + mean_right=i.right, + var_left=0, + var_right=((i.right - i.left) ** 2) / 4, + steps=steps, + shape=shape, ) + min_max = box -def min_max_mean( - minimum: Union[Interval,float,int], - maximum: Union[Interval,float,int], - mean: Union[Interval,float,int], - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + +def min_max_mean( + minimum: Union[Interval, float, int], + maximum: Union[Interval, float, int], + mean: Union[Interval, float, int], + steps: int = Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the minimum, maximum and mean of the variable - + **Parameters**: ``minimum`` : minimum value of the variable @@ -186,180 +216,181 @@ def min_max_mean( ``mean`` : mean value of the variable - + **Returns**: ``Pbox`` - ''' - mid = (maximum-mean)/(maximum-minimum) - ii = [i/steps for i in range(steps)] - left = [minimum if i <= mid else ((mean-maximum)/i + maximum) for i in ii] - jj = [j/steps for j in range(1,steps+1)] + """ + mid = (maximum - mean) / (maximum - minimum) + ii = [i / steps for i in range(steps)] + left = [minimum if i <= mid else ((mean - maximum) / i + maximum) for i in ii] + jj = [j / steps for j in range(1, steps + 1)] right = [maximum if mid <= j else (mean - minimum * j) / (1 - j) for j in jj] print(len(left)) - return Pbox( - left = np.array(left), - right = np.array(right), - steps = steps) + return Pbox(left=np.array(left), right=np.array(right), steps=steps) + def min_mean( - minimum: Union[Interval,float,int], - mean: Union[Interval,float,int], - steps = Pbox.STEPS - ) -> Pbox: - ''' + minimum: Union[Interval, float, int], + mean: Union[Interval, float, int], + steps=Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the minimum and mean of the variable - + **Parameters**: ``minimum`` : minimum value of the variable ``mean`` : mean value of the variable - + **Returns**: ``Pbox`` - ''' - jjj = np.array([j/steps for j in range(1,steps-1)] + [1-1/steps]) + """ + jjj = np.array([j / steps for j in range(1, steps - 1)] + [1 - 1 / steps]) - right = [((mean-minimum)/(1-j) + minimum) for j in jjj] + right = [((mean - minimum) / (1 - j) + minimum) for j in jjj] return Pbox( - left = np.repeat(minimum,steps), - right = right, - mean_left = mean, - mean_right = mean, - steps=steps + left=np.repeat(minimum, steps), + right=right, + mean_left=mean, + mean_right=mean, + steps=steps, ) + def max_mean( - maximum: Union[Interval,float,int], - mean: Union[Interval,float,int], - steps = Pbox.STEPS - ) -> Pbox: - ''' + maximum: Union[Interval, float, int], + mean: Union[Interval, float, int], + steps=Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the minimum and mean of the variable - + **Parameters**: ``minimum`` : minimum value of the variable ``mean`` : mean value of the variable - + **Returns**: ``Pbox`` - ''' - return min_mean(-maximum,-mean).__neg__() - + """ + return min_mean(-maximum, -mean).__neg__() def mean_std( - mean: Union[Interval,float,int], - std: Union[Interval,float,int], - steps = Pbox.STEPS - ) -> Pbox: - ''' + mean: Union[Interval, float, int], + std: Union[Interval, float, int], + steps=Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the mean and standard deviation of the variable - + **Parameters**: ``mean`` : mean of the variable ``std`` : standard deviation of the variable - + **Returns**: ``Pbox`` - ''' - iii = [1/steps] + [i/steps for i in range(1,steps-1)] - jjj = [j/steps for j in range(1,steps-1)] + [1-1/steps] - - left = [mean - std * np.sqrt(1/i - 1) for i in iii] + """ + iii = [1 / steps] + [i / steps for i in range(1, steps - 1)] + jjj = [j / steps for j in range(1, steps - 1)] + [1 - 1 / steps] + + left = [mean - std * np.sqrt(1 / i - 1) for i in iii] right = [mean + std * np.sqrt(j / (1 - j)) for j in jjj] - + return Pbox( - left = left, - right = right, - steps = steps, - mean_left = mean, - mean_right =mean, - var_left = std**2, - var_right = std**2 + left=left, + right=right, + steps=steps, + mean_left=mean, + mean_right=mean, + var_left=std**2, + var_right=std**2, ) + def mean_var( - mean: Union[Interval,float,int], - var: Union[Interval,float,int], - steps = Pbox.STEPS - ) -> Pbox: - ''' + mean: Union[Interval, float, int], + var: Union[Interval, float, int], + steps=Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the mean and variance of the variable - + Equivalent to `mean_std(mean,np.sqrt(var))` - + **Parameters**: ``mean`` : mean of the variable ``var`` : variance of the variable - + **Returns**: ``Pbox`` - ''' - return mean_std(mean,np.sqrt(var),steps) + """ + return mean_std(mean, np.sqrt(var), steps) + def pos_mean_std( - mean: Union[Interval,float,int], - std: Union[Interval,float,int], - steps = Pbox.STEPS - ) -> Pbox: - ''' + mean: Union[Interval, float, int], + std: Union[Interval, float, int], + steps=Pbox.STEPS, +) -> Pbox: + """ Generates a positive distribution-free p-box based upon the mean and standard deviation of the variable - + **Parameters**: ``mean`` : mean of the variable ``std`` : standard deviation of the variable - + **Returns**: ``Pbox`` - ''' - iii = [1/steps] + [i/steps for i in range(1,steps-1)] - jjj = [j/steps for j in range(1,steps-1)] + [1-1/steps] - - left = [max((0,mean - std * np.sqrt(1/i - 1))) for i in iii] - right = [min((mean/(1-j), mean + std * np.sqrt(j / (1 - j)))) for j in jjj] - + """ + iii = [1 / steps] + [i / steps for i in range(1, steps - 1)] + jjj = [j / steps for j in range(1, steps - 1)] + [1 - 1 / steps] + + left = [max((0, mean - std * np.sqrt(1 / i - 1))) for i in iii] + right = [min((mean / (1 - j), mean + std * np.sqrt(j / (1 - j)))) for j in jjj] + return Pbox( - left = left, - right = right, - steps = steps, - mean_left = mean, - mean_right =mean, - var_left = std**2, - var_right = std**2 + left=left, + right=right, + steps=steps, + mean_left=mean, + mean_right=mean, + var_left=std**2, + var_right=std**2, ) + def min_max_mode( - minimum: Union[Interval,float,int], - maximum: Union[Interval,float,int], - mode: Union[Interval,float,int], - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + minimum: Union[Interval, float, int], + maximum: Union[Interval, float, int], + mode: Union[Interval, float, int], + steps: int = Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the minimum, maximum, and mode of the variable - + **Parameters**: ``minimum`` : minimum value of the variable @@ -368,37 +399,37 @@ def min_max_mode( ``mode`` : mode value of the variable - + **Returns**: ``Pbox`` - ''' + """ if minimum == maximum: return box(minimum, maximum) - ii = np.array([i/steps for i in range(steps)]) - jj = np.array([j/steps for j in range(1,steps+1)]) - + ii = np.array([i / steps for i in range(steps)]) + jj = np.array([j / steps for j in range(1, steps + 1)]) + return Pbox( - left = ii * (mode - minimum) + minimum, - right = jj * (maximum - mode) + mode, - mean_left = (minimum+mode)/2, - mean_right = (mode+maximum)/2, - var_left = 0, - var_right = (maximum-minimum)*(maximum-minimum)/12 + left=ii * (mode - minimum) + minimum, + right=jj * (maximum - mode) + mode, + mean_left=(minimum + mode) / 2, + mean_right=(mode + maximum) / 2, + var_left=0, + var_right=(maximum - minimum) * (maximum - minimum) / 12, ) def min_max_median( - minimum: Union[Interval,float,int], - maximum: Union[Interval,float,int], - median: Union[Interval,float,int], - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + minimum: Union[Interval, float, int], + maximum: Union[Interval, float, int], + median: Union[Interval, float, int], + steps: int = Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the minimum, maximum and median of the variable - + **Parameters**: ``minimum`` : minimum value of the variable @@ -407,37 +438,37 @@ def min_max_median( ``median`` : median value of the variable - + **Returns**: ``Pbox`` - ''' + """ if minimum == maximum: return box(minimum, maximum) - ii = np.array([i/steps for i in range(steps)]) - jj = np.array([j/steps for j in range(1,steps+1)]) - - + ii = np.array([i / steps for i in range(steps)]) + jj = np.array([j / steps for j in range(1, steps + 1)]) + return Pbox( - left = np.array([p if p>0.5 else minimum for p in ii]), - right = np.array([p if p<=0.5 else minimum for p in jj]), - mean_left = (minimum+median)/2, - mean_right = (median+maximum)/2, - var_left = 0, - var_right = (maximum-minimum)*(maximum-minimum)/4 + left=np.array([p if p > 0.5 else minimum for p in ii]), + right=np.array([p if p <= 0.5 else minimum for p in jj]), + mean_left=(minimum + median) / 2, + mean_right=(median + maximum) / 2, + var_left=0, + var_right=(maximum - minimum) * (maximum - minimum) / 4, ) + def min_max_median_is_mode( - minimum: Union[Interval,float,int], - maximum: Union[Interval,float,int], - m: Union[Interval,float,int], - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + minimum: Union[Interval, float, int], + maximum: Union[Interval, float, int], + m: Union[Interval, float, int], + steps: int = Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the minimum, maximum and median/mode of the variable when median = mode. - + **Parameters**: ``minimum`` : minimum value of the variable @@ -446,241 +477,264 @@ def min_max_median_is_mode( ``m`` : m = median = mode value of the variable - + **Returns**: ``Pbox`` - ''' - ii = np.array([i/steps for i in range(steps)]) - jjj = [j/steps for j in range(1,steps-1)] + [1-1/steps] - - u = [p*2*(m-minimum) + minimum if p <= 0.5 else m for p in ii] - - d = [(p-0.5)*2*(maximum - m) + m if p > 0.5 else m for p in jjj] + """ + ii = np.array([i / steps for i in range(steps)]) + jjj = [j / steps for j in range(1, steps - 1)] + [1 - 1 / steps] + + u = [p * 2 * (m - minimum) + minimum if p <= 0.5 else m for p in ii] + + d = [(p - 0.5) * 2 * (maximum - m) + m if p > 0.5 else m for p in jjj] - return Pbox( - left = u, - right = d, - mean_left=(minimum + 3 + m)/4, - mean_right=(3*m +maximum)/4, + left=u, + right=d, + mean_left=(minimum + 3 + m) / 4, + mean_right=(3 * m + maximum) / 4, var_left=0, - var_right=(maximum - minimum)*(maximum - minimum)/4, + var_right=(maximum - minimum) * (maximum - minimum) / 4, ) + def symmetric_mean_std( - mean: Union[Interval,float,int], - std: Union[Interval,float,int], - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + mean: Union[Interval, float, int], + std: Union[Interval, float, int], + steps: int = Pbox.STEPS, +) -> Pbox: + """ Generates a symmetrix distribution-free p-box based upon the mean and standard deviation of the variable - + **Parameters**: ``mean`` : mean value of the variable - ``std`` : standard deviation of the variable - + ``std`` : standard deviation of the variable + **Returns** - ``Pbox`` - - ''' - iii = [1/steps] + [i/steps for i in range(1,steps-1)] - jjj = [j/steps for j in range(1,steps-1)] + [1-1/steps] - + ``Pbox`` + + """ + iii = [1 / steps] + [i / steps for i in range(1, steps - 1)] + jjj = [j / steps for j in range(1, steps - 1)] + [1 - 1 / steps] + u = [mean - std / np.sqrt(2 * p) if p <= 0.5 else mean for p in iii] d = [mean + std / np.sqrt(2 * (1 - p)) if p > 0.5 else mean for p in jjj] - + return Pbox( - left = u, - right = d, + left=u, + right=d, mean_left=mean, mean_right=mean, var_left=std**2, - var_right=std**2 + var_right=std**2, ) - def min_max_mean_std( - minimum: Union[Interval,float,int], - maximum: Union[Interval,float,int], - mean: Union[Interval,float,int], - std: Union[Interval,float,int], - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + minimum: Union[Interval, float, int], + maximum: Union[Interval, float, int], + mean: Union[Interval, float, int], + std: Union[Interval, float, int], + steps: int = Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the minimum, maximum, mean and standard deviation of the variable - + **Parameters** ``minimum`` : minimum value of the variable ``maximum`` : maximum value of the variable ``mean`` : mean value of the variable - ``std`` :standard deviation of the variable - + ``std`` :standard deviation of the variable + **Returns** ``Pbox`` .. seealso:: - + :func:`min_max_mean_var` - - ''' + + """ if minimum == maximum: return box(minimum, maximum) - def _left(x): - if type(x) in [int,float]: + + def _left(x): + if type(x) in [int, float]: return x if x.__class__.__name__ == "Interval": return x.left if x.__class__.__name__ == "Pbox": return min(x.left) - - def _right(x): - if type(x) in [int,float]: + + def _right(x): + if type(x) in [int, float]: return x if x.__class__.__name__ == "Interval": return x.right if x.__class__.__name__ == "Pbox": - return max(x.right) - - def _imp(a,b) : - return Interval(max(_left(a),_left(b)),min(_right(a),_right(b))) - def _env(a,b) : - return Interval(min(_left(a),_left(b)),max(_right(a),_right(b))) - + return max(x.right) + + def _imp(a, b): + return Interval(max(_left(a), _left(b)), min(_right(a), _right(b))) + + def _env(a, b): + return Interval(min(_left(a), _left(b)), max(_right(a), _right(b))) + def _constrain(a, b, msg): - if ((_right(a) < _left(b)) or (_right(b) < _left(a))) : + if (_right(a) < _left(b)) or (_right(b) < _left(a)): print("Math Problem: impossible constraint", msg) - return _imp(a,b) - - zero = 0.0 + return _imp(a, b) + + zero = 0.0 one = 1.0 - ran = maximum - minimum; - m = _constrain(mean, Interval(minimum,maximum), "(mean)"); - s = _constrain(std, _env(Interval(0.0),(abs(ran*ran/4.0 - (maximum-mean-ran/2.0)**2))**0.5)," (dispersion)") - ml = (m.left-minimum)/ran - sl = s.left/ran - mr = (m.right-minimum)/ran - sr = s.right/ran + ran = maximum - minimum + m = _constrain(mean, Interval(minimum, maximum), "(mean)") + s = _constrain( + std, + _env( + Interval(0.0), + (abs(ran * ran / 4.0 - (maximum - mean - ran / 2.0) ** 2)) ** 0.5, + ), + " (dispersion)", + ) + ml = (m.left - minimum) / ran + sl = s.left / ran + mr = (m.right - minimum) / ran + sr = s.right / ran z = box(minimum, maximum) - n = len(z.left) + n = len(z.left) L = [0.0] * n R = [1.0] * n - for i in range(n) : + for i in range(n): p = i / n - if (p <= zero) : + if p <= zero: x2 = zero - else : x2 = ml - sr * (one / p - one)**0.5 - if (ml + p <= one) : + else: + x2 = ml - sr * (one / p - one) ** 0.5 + if ml + p <= one: x3 = zero - else : - x5 = p*p + sl*sl - p - if (x5 >= zero) : - x4 = one - p + x5**0.5 - if (x4 < ml) : x4 = ml - else : x4 = ml - x3 = (p + sl*sl + x4*x4 - one) / (x4 + p - one) - if ((p <= zero) or (p <= (one - ml))) : x6 = zero - else : x6 = (ml - one) / p + one - L[i] = max(max(max(x2,x3),x6),zero) * ran + minimum; - - p = (i+1)/n - if (p >= one) : x2 = one - else : x2 = mr + sr * (one/(one/p - one))**0.5 - if (mr + p >= one) : x3 = one - else : - x5 = p*p + sl*sl - p - if (x5 >= zero) : - x4 = one - p - x5**0.5 - if (x4 > mr) : x4 = mr - else : x4 = mr - x3 = (p + sl*sl + x4*x4 - one) / (x4 + p - one) - one - - if (((one - mr) <= p) or (one <= p)) : x6 = one - else : x6 = mr / (one - p) - R[i] = min(min(min(x2,x3),x6),one) * ran + minimum - + else: + x5 = p * p + sl * sl - p + if x5 >= zero: + x4 = one - p + x5**0.5 + if x4 < ml: + x4 = ml + else: + x4 = ml + x3 = (p + sl * sl + x4 * x4 - one) / (x4 + p - one) + if (p <= zero) or (p <= (one - ml)): + x6 = zero + else: + x6 = (ml - one) / p + one + L[i] = max(max(max(x2, x3), x6), zero) * ran + minimum + + p = (i + 1) / n + if p >= one: + x2 = one + else: + x2 = mr + sr * (one / (one / p - one)) ** 0.5 + if mr + p >= one: + x3 = one + else: + x5 = p * p + sl * sl - p + if x5 >= zero: + x4 = one - p - x5**0.5 + if x4 > mr: + x4 = mr + else: + x4 = mr + x3 = (p + sl * sl + x4 * x4 - one) / (x4 + p - one) - one + + if ((one - mr) <= p) or (one <= p): + x6 = one + else: + x6 = mr / (one - p) + R[i] = min(min(min(x2, x3), x6), one) * ran + minimum + v = s**2 return Pbox( - left = np.array(L), - right = np.array(R), - mean_left = _left(m), - mean_right = _right(m), - var_left = _left(v), - var_right = _right(v), - steps = steps) + left=np.array(L), + right=np.array(R), + mean_left=_left(m), + mean_right=_right(m), + var_left=_left(v), + var_right=_right(v), + steps=steps, + ) + def min_max_mean_var( - minimum: Union[Interval,float,int], - maximum: Union[Interval,float,int], - mean: Union[Interval,float,int], - var: Union[Interval,float,int], - steps: int = Pbox.STEPS - ) -> Pbox: - ''' + minimum: Union[Interval, float, int], + maximum: Union[Interval, float, int], + mean: Union[Interval, float, int], + var: Union[Interval, float, int], + steps: int = Pbox.STEPS, +) -> Pbox: + """ Generates a distribution-free p-box based upon the minimum, maximum, mean and standard deviation of the variable - + **Parameters** - + ``minimum`` : minimum value of the variable ``maximum`` : maximum value of the variable ``mean`` : mean value of the variable - ``var`` :variance of the variable - + ``var`` :variance of the variable + **Returns** ``Pbox`` .. admonition:: Implementation - + Equivalent to ``min_max_mean_std(minimum,maximum,mean,np.sqrt(var))`` - - .. seealso:: - + + .. seealso:: + :func:`min_max_mean_std` - ''' - return min_max_mean_std(minimum,maximum,mean,np.sqrt(var)) + """ + return min_max_mean_std(minimum, maximum, mean, np.sqrt(var)) + def from_percentiles(percentiles: dict, steps: int = Pbox.STEPS) -> Pbox: - ''' + """ Generates a distribution-free p-box based upon percentiles of the variable - + **Parameters** ``percentiles`` : dictionary of percentiles and their values (e.g. {0: 0, 0.1: 1, 0.5: 2, 0.9: Interval(3,4), 1:5}) - + ``steps`` : number of steps to use in the p-box - + .. important:: - + The percentiles dictionary is of the form {percentile: value}. Where value can either be a number or an Interval. If value is a number, the percentile is assumed to be a point percentile. If value is an Interval, the percentile is assumed to be an interval percentile. - + .. warning:: - + If no keys for 0 and 1 are given, ``-np.inf`` and ``np.inf`` are used respectively. This will result in a p-box that is not bounded and raise a warning. - + If the percentiles are not increasing, the percentiles will be intersected. This may not be desired behaviour. - + .. error:: - + ``ValueError``: If any of the percentiles are not between 0 and 1. - + **Returns** - + ``Pbox`` - + **Example**: - + .. code-block:: python - + pba.from_percentiles( {0: 0, 0.25: 0.5, @@ -692,8 +746,8 @@ def from_percentiles(percentiles: dict, steps: int = Pbox.STEPS) -> Pbox: .. image:: https://github.com/Institute-for-Risk-and-Uncertainty/pba-for-python/blob/master/docs/images/from_percentiles.png?raw=true :scale: 35 % :align: center - - ''' + + """ # check if 0 and 1 are in the dictionary if 0 not in percentiles.keys(): percentiles[0] = -np.inf @@ -701,78 +755,84 @@ def from_percentiles(percentiles: dict, steps: int = Pbox.STEPS) -> Pbox: if 1 not in percentiles.keys(): percentiles[1] = np.inf warn("No value given for 1 percentile. Using np.inf") - + # sort the dictionary by percentile percentiles = dict(sorted(percentiles.items())) - + # transform values to intervals - for k,v in percentiles.items(): - if not isinstance(v,Interval): + for k, v in percentiles.items(): + if not isinstance(v, Interval): percentiles[k] = Interval(v) - - if any([p<0 or p>1 for p in percentiles.keys()]): + + if any([p < 0 or p > 1 for p in percentiles.keys()]): raise ValueError("Percentiles must be between 0 and 1") - - X = np.linspace(0,1,steps) - + + X = np.linspace(0, 1, steps) + if any([p not in percentiles.keys() for p in X]): - warn("Not all percentiles are possible with given steps. This may result in a p-box that does not cover the entire range.") - + warn( + "Not all percentiles are possible with given steps. This may result in a p-box that does not cover the entire range." + ) + left = [] right = [] for i in X: - smallest_key = min(key for key in percentiles.keys() if key >= i) - largest_key = max(key for key in percentiles.keys() if key <= i) + smallest_key = min(key for key in percentiles.keys() if key >= i) + largest_key = max(key for key in percentiles.keys() if key <= i) left.append(percentiles[largest_key].left) right.append(percentiles[smallest_key].right) - + try: - return Pbox(left,right,steps=steps, interpolation='outer') + return Pbox(left, right, steps=steps, interpolation="outer") except NotIncreasingError: warn("Percentiles are not increasing. Will take intersection of percentiles.") - + left = [] right = [] p = list(percentiles.keys()) - for i,j,k in zip(p,p[1:],p[2:]): + for i, j, k in zip(p, p[1:], p[2:]): if sometimes(percentiles[j] < percentiles[i]): - percentiles[j] = Interval(percentiles[i].right,percentiles[j].right) + percentiles[j] = Interval(percentiles[i].right, percentiles[j].right) if sometimes(percentiles[j] > percentiles[k]): - percentiles[j] = Interval(percentiles[j].left,percentiles[k].left) + percentiles[j] = Interval(percentiles[j].left, percentiles[k].left) - for i in np.linspace(0,1,steps): + for i in np.linspace(0, 1, steps): left_key = max(key for key in percentiles.keys() if key <= i) - right_key = min(key for key in percentiles.keys() if key >= i) + right_key = min(key for key in percentiles.keys() if key >= i) left.append(percentiles[left_key].left) right.append(percentiles[right_key].right) - - return Pbox(left,right,steps=steps, interpolation='outer') + + return Pbox(left, right, steps=steps, interpolation="outer") except: raise Exception("Unable to generate p-box") + ### ML-ME -''' +""" Maximum Likelihood methods __________________________ Maximum likelihood estimation (MLE) is a method of estimating the parameters of a probability distribution by maximizing a likelihood function, so that under the assumed statistical model the observed data is most probable. The point in the parameter space that maximizes the likelihood function is called the maximum likelihood estimate. The logic of maximum likelihood is both intuitive and flexible, and as such the method has become a dominant means of statistical inference. -''' -def MLnorm(data): - return norm(np.mean(data),np.std(data)) +""" + + +def MLnorm(data): + return norm(np.mean(data), np.std(data)) + def ME_min_max_mean_std( - minimum: Union[Interval,float,int], - maximum: Union[Interval,float,int], - mean: Union[Interval,float,int], - stddev: Union[Interval,float,int], - steps: int = Pbox.STEPS - ) -> Pbox: - - μ = ((mean- minimum) / (maximum - minimum)) - - σ = (stddev/(maximum - minimum) ) - - a = ((1-μ)/(σ**2) - 1/μ)*μ**2 - b = a*(1/μ - 1) - - return beta(a,b,steps=steps)* (maximum - minimum) + minimum \ No newline at end of file + minimum: Union[Interval, float, int], + maximum: Union[Interval, float, int], + mean: Union[Interval, float, int], + stddev: Union[Interval, float, int], + steps: int = Pbox.STEPS, +) -> Pbox: + + μ = (mean - minimum) / (maximum - minimum) + + σ = stddev / (maximum - minimum) + + a = ((1 - μ) / (σ**2) - 1 / μ) * μ**2 + b = a * (1 / μ - 1) + + return beta(a, b, steps=steps) * (maximum - minimum) + minimum diff --git a/pba/pbox_constructors/parametric.py b/pba/pbox_constructors/parametric.py index 290a15b..476118c 100644 --- a/pba/pbox_constructors/parametric.py +++ b/pba/pbox_constructors/parametric.py @@ -1,4 +1,4 @@ -''' +""" .. _pbox_constructors.sp: Automatically generated @@ -7,7 +7,7 @@ P-boxes based on scipy.stats distributions. These p-boxes are generated using the scipy.stats distributions. The p-boxes are generated by taking the minimum and maximum values of the inverse CDF of the distribution for a given set of parameters. -''' +""" from ..interval import * from ..pbox import * @@ -18,7 +18,7 @@ import sys from copy import deepcopy -#* The commented out distributions are overwritten in dists +# * The commented out distributions are overwritten in dists dist = [ "alpha", "anglit", @@ -131,58 +131,56 @@ "randint", "skellam", "zipf", - "yulesimon" - ] + "yulesimon", +] class parametric: def __init__(self, function_name): self.function_name = function_name - - def make_pbox(self, *args, steps = 200, support = Interval(1e-4,1-1e-4), **kwargs) -> Pbox: - f''' + + def make_pbox( + self, *args, steps=200, support=Interval(1e-4, 1 - 1e-4), **kwargs + ) -> Pbox: + f""" Generate a P-box from a scipy.stats.{self.function_name} distribution. - ''' + """ args = list(args) for i, a in enumerate(args): if not isinstance(a, Interval): args[i] = Interval(a) # define support - x = np.linspace(support.left,support.right,steps) + x = np.linspace(support.left, support.right, steps) - #get bound arguments + # get bound arguments new_args = list(product(*args)) - bounds = np.empty((len(new_args),steps)) + bounds = np.empty((len(new_args), steps)) means = [] variances = [] - for i,a in enumerate(new_args): + for i, a in enumerate(new_args): - bounds[i,:] = stats.__dict__[self.function_name].ppf(x,*a,**kwargs) - m, v = stats.__dict__[self.function_name].stats(*a,**kwargs,moments = 'mv') + bounds[i, :] = stats.__dict__[self.function_name].ppf(x, *a, **kwargs) + m, v = stats.__dict__[self.function_name].stats(*a, **kwargs, moments="mv") means.append(m) variances.append(v) - Left = np.array([min([b[i] for b in bounds]) for i in range(steps)]) Right = np.array([max([b[i] for b in bounds]) for i in range(steps)]) - mean = Interval(min(means),max(means)) - var = Interval(min(variances),max(variances)) - + mean = Interval(min(means), max(means)) + var = Interval(min(variances), max(variances)) + return Pbox( - left=Left, - right=Right, - mean = mean, - var = var, - shape= self.function_name + left=Left, right=Right, mean=mean, var=var, shape=self.function_name ) - + + for d in dist: locals()[d] = parametric(d).make_pbox -__all__ = deepcopy(dist) \ No newline at end of file +__all__ = deepcopy(dist) diff --git a/pyproject.toml b/pyproject.toml index 082bfe6..62d0315 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,7 @@ authors = [ {name = "Nick Gray", email = "ngg@liv.ac.uk"} ] description = "Probability Bounds Analysis in Python" -readme = "README.rst" +readme = "README.md" license = {file = "LICENSE"} [project.urls] Documentation = "https://pba-for-python.readthedocs.io/en/latest/"