diff --git a/mcerp/__init__.py b/mcerp/__init__.py index 50f194c..d27c4a4 100644 --- a/mcerp/__init__.py +++ b/mcerp/__init__.py @@ -19,7 +19,7 @@ npts = 10000 -CONSTANT_TYPES = (float, int, long) +CONSTANT_TYPES = (float, int, long, np.ndarray, np.float64) class NotUpcast(Exception): 'Raised when an object cannot be converted to a number with uncertainty' @@ -27,7 +27,7 @@ class NotUpcast(Exception): def to_uncertain_func(x): """ Transforms x into an UncertainFunction-compatible object, - unless it is already an UncertainFunction (in which case x is returned + unless it is already an UncertainFunction (in which case x is returned unchanged). Raises an exception unless 'x' belongs to some specific classes of @@ -41,21 +41,22 @@ def to_uncertain_func(x): elif isinstance(x, CONSTANT_TYPES): # No variable => no derivative to define: return UncertainFunction([x]*npts) - + + print 'Picking up new version' raise NotUpcast("%s cannot be converted to a number with" " uncertainty" % type(x)) - - + + class UncertainFunction(object): """ - UncertainFunction objects represent the uncertainty of a result of + UncertainFunction objects represent the uncertainty of a result of calculations with uncertain variables. Nearly all basic mathematical operations are supported. - + This class is mostly intended for internal use. - - + + """ def __init__(self, mcpts): self._mcpts = np.atleast_1d(mcpts).flatten() @@ -68,7 +69,7 @@ def mean(self): """ mn = np.mean(self._mcpts) return mn - + @property def var(self): """ @@ -77,42 +78,42 @@ def var(self): mn = self.mean vr = np.mean((self._mcpts - mn)**2) return vr - + @property def std(self): """ - Standard deviation value as a result of an uncertainty calculation, + Standard deviation value as a result of an uncertainty calculation, defined as:: - + ________ std = \/variance - + """ return self.var**0.5 - + @property def skew(self): """ Skewness coefficient value as a result of an uncertainty calculation, defined as:: - + _____ m3 \/beta1 = ------ std**3 - + where m3 is the third central moment and std is the standard deviation """ mn = self.mean sd = self.std sk = 0.0 if abs(sd)<=1e-8 else np.mean((self._mcpts - mn)**3)/sd**3 return sk - + @property def kurt(self): """ Kurtosis coefficient value as a result of an uncertainty calculation, defined as:: - + m4 beta2 = ------ std**4 @@ -123,7 +124,7 @@ def kurt(self): sd = self.std kt = 0.0 if abs(sd)<=1e-8 else np.mean((self._mcpts - mn)**4)/sd**4 return kt - + @property def stats(self): """ @@ -135,23 +136,23 @@ def stats(self): sk = self.skew kt = self.kurt return [mn, vr, sk, kt] - + def percentile(self, val): """ Get the distribution value at a given percentile or set of percentiles. This follows the NIST method for calculating percentiles. - + Parameters ---------- val : scalar or array Either a single value or an array of values between 0 and 1. - + Returns ------- out : scalar or array The actual distribution value that appears at the requested percentile value or values - + """ try: # test to see if an input is given as an array @@ -169,11 +170,11 @@ def percentile(self, val): if isinstance(val, np.ndarray): out = np.array(out) return out - + def _to_general_representation(self,str_func): mn, vr, sk, kt = self.stats return ('uv({:}, {:}, {:}, {:})'.format( - str_func(mn), str_func(vr), str_func(sk), str_func(kt)) + str_func(mn), str_func(vr), str_func(sk), str_func(kt)) if any([vr, sk, kt]) else str_func(mn)) def __str__(self): @@ -190,17 +191,17 @@ def describe(self, name=None): - Variance - Standardized Skewness Coefficient - Standardized Kurtosis Coefficient - + For a standard Normal distribution, these are [0, 1, 0, 3]. - + If the object has an associated tag, this is presented. If the optional ``name`` kwarg is utilized, this is presented as with the moments. Otherwise, no unique name is presented. - + Example ======= :: - + >>> x = N(0, 1, 'x') >>> x.describe() # print tag since assigned MCERP Uncertain Value (x): @@ -209,7 +210,7 @@ def describe(self, name=None): >>> x.describe('foobar') # 'name' kwarg takes precedence MCERP Uncertain Value (foobar): ... - + >>> y = x**2 >>> y.describe('y') # print name since assigned MCERP Uncertain Value (y): @@ -232,22 +233,22 @@ def describe(self, name=None): s += ' > Skewness Coefficient... {: }\n'.format(sk) s += ' > Kurtosis Coefficient... {: }\n'.format(kt) print s - + def plot(self, hist=False, show=False, **kwargs): """ Plot the distribution of the UncertainFunction. By default, the distribution is shown with a kernel density estimate (kde). - + Optional -------- hist : bool If true, a density histogram is displayed (histtype='stepfilled') show : bool - If ``True``, the figure will be displayed after plotting the + If ``True``, the figure will be displayed after plotting the distribution. If ``False``, an explicit call to ``plt.show()`` is required to display the figure. kwargs : any valid matplotlib.pyplot.plot or .hist kwarg - + """ vals = self._mcpts low = min(vals) @@ -257,14 +258,14 @@ def plot(self, hist=False, show=False, **kwargs): xp = np.linspace(low,high,100) if hist: - h = plt.hist(vals, bins=np.round(np.sqrt(len(vals))), + h = plt.hist(vals, bins=np.round(np.sqrt(len(vals))), histtype='stepfilled', normed=True, **kwargs) plt.ylim(0, 1.1*h[0].max()) else: plt.plot(xp, p.evaluate(xp), **kwargs) plt.xlim(low - (high - low)*0.1, high + (high - low)*0.1) - + if show: self.show() @@ -272,77 +273,161 @@ def show(self): plt.show() def __add__(self, val): + if isinstance(val, np.ndarray): + uf = to_uncertain_func(self) + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts + uf._mcpts )) + return np.array(f(val)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[0]._mcpts + uf[1]._mcpts return UncertainFunction(mcpts) def __radd__(self, val): + if isinstance(self, np.ndarray): + uf = to_uncertain_func(val) + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts + uf._mcpts )) + return np.array(f(self)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[0]._mcpts + uf[1]._mcpts return UncertainFunction(mcpts) - + def __mul__(self, val): + if isinstance(val, np.ndarray): + uf = to_uncertain_func(self) + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts * uf._mcpts )) + return np.array(f(val)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[0]._mcpts * uf[1]._mcpts return UncertainFunction(mcpts) + def __rmul__(self, val): + if isinstance(self, np.ndarray): + uf = to_uncertain_func(val) + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts * uf._mcpts )) + return np.array(f(self)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[0]._mcpts * uf[1]._mcpts return UncertainFunction(mcpts) - + def __sub__(self, val): + if isinstance(val, np.ndarray): + uf = to_uncertain_func(self) + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts - uf._mcpts )) + return np.array(f(val)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[0]._mcpts - uf[1]._mcpts return UncertainFunction(mcpts) def __rsub__(self, val): + if isinstance(self, np.ndarray): + uf = to_uncertain_func(val) + f = np.vectorize( lambda x: UncertainFunction( uf._mcpts - to_uncertain_func(x)._mcpts )) + return np.array(f(self)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[1]._mcpts - uf[0]._mcpts return UncertainFunction(mcpts) - + def __div__(self, val): + if isinstance(val, np.ndarray): + uf = to_uncertain_func(self) + f = np.vectorize( lambda x: UncertainFunction( uf._mcpts / to_uncertain_func(x)._mcpts )) + return np.array(f(val)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[0]._mcpts/uf[1]._mcpts return UncertainFunction(mcpts) def __rdiv__(self, val): + if isinstance(self, np.ndarray): + uf = to_uncertain_func(val) + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts / uf._mcpts )) + return np.array(f(self)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[1]._mcpts/uf[0]._mcpts return UncertainFunction(mcpts) - + def __truediv__(self, val): + if isinstance(val, np.ndarray): + uf = to_uncertain_func(self) + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts / uf._mcpts )) + return np.array(f(val)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[0]._mcpts/uf[1]._mcpts return UncertainFunction(mcpts) def __rtruediv__(self, val): + if isinstance(self, np.ndarray): + uf = to_uncertain_func(val) + f = np.vectorize( lambda x: UncertainFunction( uf._mcpts / to_uncertain_func(x)._mcpts )) + return np.array(f(self)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[1]._mcpts/uf[0]._mcpts return UncertainFunction(mcpts) - + def __pow__(self, val): + if isinstance(val, np.ndarray): + uf = to_uncertain_func(self) + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts ** uf._mcpts )) + return np.array(f(val)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[0]._mcpts**uf[1]._mcpts return UncertainFunction(mcpts) + def pow(self, val): + return __pow__(self, val) + def __rpow__(self, val): + if isinstance(self, np.ndarray): + uf = to_uncertain_func(val) + f = np.vectorize( lambda x: UncertainFunction( uf._mcpts ** to_uncertain_func(x)._mcpts )) + return np.array(f(self)) + uf = map(to_uncertain_func, [self, val]) mcpts = uf[1]._mcpts**uf[0]._mcpts return UncertainFunction(mcpts) - + def __neg__(self): + if isinstance(self, np.ndarray): + f = np.vectorize( lambda x: UncertainFunction( -to_uncertain_func(x)._mcpts )) + return np.array(f(self)) + mcpts = -self._mcpts return UncertainFunction(mcpts) - + + def exp(self): + if isinstance(self, np.ndarray): + f = np.vectorize( lambda x: UncertainFunction( np.exp(to_uncertain_func(x)._mcpts) )) + return np.array(f(self)) + + mcpts = np.exp(self._mcpts) + return UncertainFunction(mcpts) + def __pos__(self): + if isinstance(self, np.ndarray): + f = np.vectorize( lambda x: UncertainFunction( to_uncertain_func(x)._mcpts )) + return np.array(f(self)) + mcpts = self._mcpts return UncertainFunction(mcpts) - + def __abs__(self): + if isinstance(self, np.ndarray): + f = np.vectorize( lambda x: UncertainFunction( np.abs(to_uncertain_func(x)._mcpts) )) + return np.array(f(self)) + mcpts = np.abs(self._mcpts) return UncertainFunction(mcpts) - + def __eq__(self, val): """ If we are comparing two distributions, check the resulting moments. If @@ -350,21 +435,21 @@ def __eq__(self, val): we can know that it is actually the same distribution we are comparing ``self`` to, otherwise, at least one statistical moment will be non- zero. - + If we are comparing ``self`` to a scalar, just do a normal comparison so that if the underlying distribution looks like a PMF, a meaningful probability of self==val is returned. This can still work quite safely for PDF distributions since the likelihood of comparing self to an actual sampled value is negligible when mcerp.npts is large. - + Examples: - + >>> h = H(50, 5, 10) # Hypergeometric distribution (PMF) >>> h==4 # what's the probability of getting 4 of the 5? 0.004 >>> sum([h==i for i in (0, 1, 2, 3, 4, 5)]) # sum of all discrete probabilities 1.0 - + >>> n = N(0, 1) # Normal distribution (PDF) >>> n==0 # what's the probability of being exactly 0.0? 0.0 @@ -380,13 +465,13 @@ def __eq__(self, val): return not (diff.mean or diff.var or diff.skew or diff.kurt) else: return len(self._mcpts[self._mcpts==val])/float(npts) - + def __ne__(self, val): if isinstance(val, UncertainFunction): return not self==val else: return 1 - (self==val) - + def __lt__(self, val): """ If we are comparing two distributions, perform statistical tests, @@ -402,13 +487,13 @@ def __lt__(self, val): return True if sgn==-1 else False else: return len(self._mcpts[self._mcptsval @@ -438,8 +523,8 @@ def __nonzero__(self): class UncertainVariable(UncertainFunction): """ - UncertainVariable objects track the effects of uncertainty, characterized - in terms of the first four standard moments of statistical distributions + UncertainVariable objects track the effects of uncertainty, characterized + in terms of the first four standard moments of statistical distributions (mean, variance, skewness and kurtosis coefficients). Monte Carlo simulation, in conjunction with Latin-hypercube based sampling performs the calculations. @@ -447,34 +532,34 @@ class UncertainVariable(UncertainFunction): ---------- rv : scipy.stats.distribution A distribution to characterize the uncertainty - + tag : str, optional A string identifier when information about this variable is printed to the screen - + Notes ----- - + The ``scipy.stats`` module contains many distributions which we can use to perform any necessary uncertainty calculation. It is important to follow the initialization syntax for creating any kind of distribution object: - - - *Location* and *Scale* values must use the kwargs ``loc`` and + + - *Location* and *Scale* values must use the kwargs ``loc`` and ``scale`` - - *Shape* values are passed in as non-keyword arguments before the + - *Shape* values are passed in as non-keyword arguments before the location and scale, (see below for syntax examples).. - - The mathematical operations that can be performed on uncertain objects will - work for any distribution supplied, but may be misleading if the supplied - moments or distribution is not accurately defined. Here are some guidelines - for creating UncertainVariable objects using some of the most common + + The mathematical operations that can be performed on uncertain objects will + work for any distribution supplied, but may be misleading if the supplied + moments or distribution is not accurately defined. Here are some guidelines + for creating UncertainVariable objects using some of the most common statistical distributions: - + +---------------------------+-------------+-------------------+-----+---------+ | Distribution | scipy.stats | args | loc | scale | | | class name | (shape params) | | | +===========================+=============+===================+=====+=========+ - | Normal(mu, sigma) | norm | | mu | sigma | + | Normal(mu, sigma) | norm | | mu | sigma | +---------------------------+-------------+-------------------+-----+---------+ | Uniform(a, b) | uniform | | a | b-a | +---------------------------+-------------+-------------------+-----+---------+ @@ -506,19 +591,19 @@ class UncertainVariable(UncertainFunction): +---------------------------+-------------+-------------------+-----+---------+ | Poisson(lamda) | poisson | lamda | | | +---------------------------+-------------+-------------------+-----+---------+ - + Thus, each distribution above would have the same call signature:: - + >>> import scipy.stats as ss >>> ss.your_dist_here(args, loc=loc, scale=scale) - + ANY SCIPY.STATS.DISTRIBUTION SHOULD WORK! IF ONE DOESN'T, PLEASE LET ME KNOW! - - Convenient constructors have been created to make assigning these + + Convenient constructors have been created to make assigning these distributions easier. They follow the parameter notation found in the respective Wikipedia articles: - + +---------------------------+---------------------------------------------------------------+ | MCERP Distibution | Wikipedia page | +===========================+===============================================================+ @@ -563,11 +648,11 @@ class UncertainVariable(UncertainFunction): Examples -------- A three-part assembly - + >>> x1 = N(24, 1) >>> x2 = N(37, 4) >>> x3 = Exp(2) # Exp(mu=0.5) works too - + >>> Z = (x1*x2**2)/(15*(1.5 + x3)) >>> Z uv(1161.46231679, 116646.762981, 0.345533974771, 3.00791101068) @@ -575,89 +660,89 @@ class UncertainVariable(UncertainFunction): The result shows the mean, variance, and standardized skewness and kurtosis of the output variable Z, which will vary from use to use due to the random nature of Monte Carlo simulation and latin-hypercube sampling techniques. - - Basic math operations may be applied to distributions, where all + + Basic math operations may be applied to distributions, where all statistical calculations are performed using latin-hypercube enhanced Monte - Carlo simulation. Nearly all of the built-in trigonometric-, logarithm-, - etc. functions of the ``math`` module have uncertainty-compatible - counterparts that should be used when possible since they support both - scalar values and uncertain objects. These can be used after importing the + Carlo simulation. Nearly all of the built-in trigonometric-, logarithm-, + etc. functions of the ``math`` module have uncertainty-compatible + counterparts that should be used when possible since they support both + scalar values and uncertain objects. These can be used after importing the ``umath`` module:: - + >>> from mcerp.umath import * # sin(), sqrt(), etc. >>> sqrt(x1) uv(4.89791765647, 0.0104291897681, -0.0614940614672, 3.00264937735) - + At any time, the standardized statistics can be retrieved using:: - + >>> x1.mean >>> x1.var # x1.std (standard deviation) is also available >>> x1.skew >>> x1.kurt - + or all four together with:: - + >>> x1.stats - + By default, the Monte Carlo simulation uses 10000 samples, but this can be changed at any time with:: - + >>> mcerp.npts = number_of_samples - + Any value from 1,000 to 1,000,000 is recommended (more samples means more - accurate, but also means more time required to perform the calculations). + accurate, but also means more time required to perform the calculations). Although it can be changed, since variables retain their samples from one - calculation to the next, this parameter should be changed before any - calculations are performed to ensure parameter compatibility (this may + calculation to the next, this parameter should be changed before any + calculations are performed to ensure parameter compatibility (this may change to be more dynamic in the future, but for now this is how it is). - + Also, to see the underlying distribution of the variable, and if matplotlib is installed, simply call its plot method:: - + >>> x1.plot() - + Optional kwargs can be any valid kwarg used by matplotlib.pyplot.plot - + See Also -------- - N, U, Exp, Gamma, Beta, LogN, X2, F, Tri, PERT, T, Weib, Bern, B, G, H, + N, U, Exp, Gamma, Beta, LogN, X2, F, Tri, PERT, T, Weib, Bern, B, G, H, Pois - + """ - + def __init__(self, rv, tag=None): - + assert hasattr(rv, 'dist'), 'Input must be a distribution from ' + \ 'the scipy.stats module.' self.rv = rv - + # generate the latin-hypercube points self._mcpts = lhd(dist=self.rv, size=npts).flatten() self.tag = tag - + def plot(self, hist=False, show=False, **kwargs): """ - Plot the distribution of the UncertainVariable. Continuous + Plot the distribution of the UncertainVariable. Continuous distributions are plotted with a line plot and discrete distributions are plotted with discrete circles. - + Optional -------- hist : bool If true, a histogram is displayed show : bool - If ``True``, the figure will be displayed after plotting the + If ``True``, the figure will be displayed after plotting the distribution. If ``False``, an explicit call to ``plt.show()`` is required to display the figure. kwargs : any valid matplotlib.pyplot.plot kwarg - + """ - + if hist: vals = self._mcpts low = vals.min() high = vals.max() - h = plt.hist(vals, bins=np.round(np.sqrt(len(vals))), + h = plt.hist(vals, bins=np.round(np.sqrt(len(vals))), histtype='stepfilled', normed=True, **kwargs) plt.ylim(0, 1.1*h[0].max()) else: @@ -673,11 +758,11 @@ def plot(self, hist=False, show=False, **kwargs): vals = np.linspace(low, high, 500) plt.plot(vals, self.rv.pdf(vals), **kwargs) plt.xlim(low - (high - low)*0.1, high + (high - low)*0.1) - + if show: self.show() - - + + uv = UncertainVariable # a nicer form for the user # DON'T MOVE THIS IMPORT!!! The prior definitions must be in place before @@ -688,7 +773,7 @@ def plot(self, hist=False, show=False, **kwargs): ############################################################################### # Define some convenience constructors for common statistical distributions. -# Hopefully these are a little easier/more intuitive to use than the +# Hopefully these are a little easier/more intuitive to use than the # scipy.stats.distributions. ############################################################################### @@ -699,14 +784,14 @@ def plot(self, hist=False, show=False, **kwargs): def Beta(alpha, beta, low=0, high=1, tag=None): """ A Beta random variate - + Parameters ---------- alpha : scalar The first shape parameter beta : scalar The second shape parameter - + Optional -------- low : scalar @@ -723,14 +808,14 @@ def Beta(alpha, beta, low=0, high=1, tag=None): def BetaPrime(alpha, beta, tag=None): """ A BetaPrime random variate - + Parameters ---------- alpha : scalar The first shape parameter beta : scalar The second shape parameter - + """ assert alpha>0 and beta>0, 'BetaPrime "alpha" and "beta" parameters must be greater than zero' x = Beta(alpha, beta, tag) @@ -741,7 +826,7 @@ def BetaPrime(alpha, beta, tag=None): def Bradford(q, low=0, high=1, tag=None): """ A Bradford random variate - + Parameters ---------- q : scalar @@ -760,14 +845,14 @@ def Bradford(q, low=0, high=1, tag=None): def Burr(c, k, tag=None): """ A Burr random variate - + Parameters ---------- c : scalar The first shape parameter k : scalar The second shape parameter - + """ assert c>0 and k>0, 'Burr "c" and "k" parameters must be greater than zero' return uv(ss.burr(c, k), tag=tag) @@ -777,7 +862,7 @@ def Burr(c, k, tag=None): def ChiSquared(k, tag=None): """ A Chi-Squared random variate - + Parameters ---------- k : int @@ -792,12 +877,12 @@ def ChiSquared(k, tag=None): def Erf(h, tag=None): """ - An Error Function random variate. - - This distribution is derived from a normal distribution by setting - m = 0 and s = 1/(h*sqrt(2)), and thus is used in similar situations + An Error Function random variate. + + This distribution is derived from a normal distribution by setting + m = 0 and s = 1/(h*sqrt(2)), and thus is used in similar situations as the normal distribution. - + Parameters ---------- h : scalar @@ -811,12 +896,12 @@ def Erf(h, tag=None): def Erlang(k, lamda, tag=None): """ An Erlang random variate. - - This distribution is the same as a Gamma(k, theta) distribution, but + + This distribution is the same as a Gamma(k, theta) distribution, but with the restriction that k must be a positive integer. This is provided for greater compatibility with other simulation tools, but provides no advantage over the Gamma distribution in its applications. - + Parameters ---------- k : int @@ -833,7 +918,7 @@ def Erlang(k, lamda, tag=None): def Exponential(lamda, tag=None): """ An Exponential random variate - + Parameters ---------- lamda : scalar @@ -849,7 +934,7 @@ def Exponential(lamda, tag=None): def ExtValueMax(mu, sigma, tag=None): """ An Extreme Value Maximum random variate. - + Parameters ---------- mu : scalar @@ -868,7 +953,7 @@ def ExtValueMax(mu, sigma, tag=None): def ExtValueMin(mu, sigma, tag=None): """ An Extreme Value Minimum random variate. - + Parameters ---------- mu : scalar @@ -887,7 +972,7 @@ def ExtValueMin(mu, sigma, tag=None): def Fisher(d1, d2, tag=None): """ An F (fisher) random variate - + Parameters ---------- d1 : int @@ -906,7 +991,7 @@ def Fisher(d1, d2, tag=None): def Gamma(k, theta, tag=None): """ A Gamma random variate - + Parameters ---------- k : scalar @@ -922,7 +1007,7 @@ def Gamma(k, theta, tag=None): def LogNormal(mu, sigma, tag=None): """ A Log-Normal random variate - + Parameters ---------- mu : scalar @@ -940,7 +1025,7 @@ def LogNormal(mu, sigma, tag=None): def Normal(mu, sigma, tag=None): """ A Normal (or Gaussian) random variate - + Parameters ---------- mu : scalar @@ -958,7 +1043,7 @@ def Normal(mu, sigma, tag=None): def Pareto(q, a, tag=None): """ A Pareto random variate (first kind) - + Parameters ---------- q : scalar @@ -976,7 +1061,7 @@ def Pareto2(q, b, tag=None): """ A Pareto random variate (second kind). This form always starts at the origin. - + Parameters ---------- q : scalar @@ -992,7 +1077,7 @@ def Pareto2(q, b, tag=None): def PERT(low, peak, high, g=4.0, tag=None): """ A PERT random variate - + Parameters ---------- low : scalar @@ -1001,12 +1086,12 @@ def PERT(low, peak, high, g=4.0, tag=None): The location of the distribution's peak (low <= peak <= high) high : scalar Upper bound of the distribution support - + Optional -------- g : scalar Controls the uncertainty of the distribution around the peak. Smaller - values make the distribution flatter and more uncertain around the + values make the distribution flatter and more uncertain around the peak while larger values make it focused and less uncertain around the peak. (Default: 4) """ @@ -1019,7 +1104,7 @@ def PERT(low, peak, high, g=4.0, tag=None): else: a1 = ((mu - a)*(2*b - a - c))/((b - mu)*(c - a)) a2 = a1*(c - mu)/(mu - a) - + return Beta(a1, a2, a, c, tag) ############################################################################### @@ -1027,7 +1112,7 @@ def PERT(low, peak, high, g=4.0, tag=None): def StudentT(v, tag=None): """ A Student-T random variate - + Parameters ---------- v : int @@ -1043,7 +1128,7 @@ def StudentT(v, tag=None): def Triangular(low, peak, high, tag=None): """ A triangular random variate - + Parameters ---------- low : scalar @@ -1055,7 +1140,7 @@ def Triangular(low, peak, high, tag=None): """ assert low<=peak<=high, 'Triangular "peak" must lie between "low" and "high"' low, peak, high = [float(x) for x in [low, peak, high]] - return uv(ss.triang((1.0*peak - low)/(high - low), loc=low, + return uv(ss.triang((1.0*peak - low)/(high - low), loc=low, scale=(high - low)), tag=tag) Tri = Triangular # for more concise use @@ -1065,7 +1150,7 @@ def Triangular(low, peak, high, tag=None): def Uniform(low, high, tag=None): """ A Uniform random variate - + Parameters ---------- low : scalar @@ -1083,7 +1168,7 @@ def Uniform(low, high, tag=None): def Weibull(lamda, k, tag=None): """ A Weibull random variate - + Parameters ---------- lamda : scalar @@ -1103,7 +1188,7 @@ def Weibull(lamda, k, tag=None): def Bernoulli(p, tag=None): """ A Bernoulli random variate - + Parameters ---------- p : scalar @@ -1119,7 +1204,7 @@ def Bernoulli(p, tag=None): def Binomial(n, p, tag=None): """ A Binomial random variate - + Parameters ---------- n : int @@ -1138,7 +1223,7 @@ def Binomial(n, p, tag=None): def Geometric(p, tag=None): """ A Geometric random variate - + Parameters ---------- p : scalar @@ -1154,7 +1239,7 @@ def Geometric(p, tag=None): def Hypergeometric(N, n, K, tag=None): """ A Hypergeometric random variate - + Parameters ---------- N : int @@ -1163,7 +1248,7 @@ def Hypergeometric(N, n, K, tag=None): The number of individuals of interest in the population K : int The number of individuals that will be chosen from the population - + Example ------- (Taken from the wikipedia page) Assume we have an urn with two types of @@ -1171,19 +1256,19 @@ def Hypergeometric(N, n, K, tag=None): close your eyes and draw 10 marbles without replacement. What is the probability that exactly 4 of the 10 are white? :: - + >>> black = 45 >>> white = 5 >>> draw = 10 - + # Now we create the distribution >>> h = H(black + white, white, draw) - + # To check the probability, in this case, we can use the underlying # scipy.stats object >>> h.rv.pmf(4) # What is the probability that white count = 4? 0.0039645830580151975 - + """ assert int(N)==N and N>0, 'Hypergeometric total population size "N" must be an integer greater than zero.' assert int(n)==n and 0>> x = N(1, 0.1) >>> y = N(10, 0.1) >>> z = x + 2*y @@ -1249,7 +1334,7 @@ def covariance_matrix(nums_with_uncert): coef = np.mean((expr1._mcpts - mean1)*(expr2._mcpts - mean2)) coefs_expr1.append(coef) cov_matrix.append(coefs_expr1) - + # We symmetrize the matrix: for (i, covariance_coefs) in enumerate(cov_matrix): covariance_coefs.extend(cov_matrix[j][i] @@ -1263,20 +1348,20 @@ def correlation_matrix(nums_with_uncert): """ Calculate the correlation matrix of uncertain variables, oriented by the order of the inputs - + Parameters ---------- nums_with_uncert : array-like A list of variables that have an associated uncertainty - + Returns ------- corr_matrix : 2d-array-like A nested list containing covariance values - + Example ------- - + >>> x = N(1, 0.1) >>> y = N(10, 0.1) >>> z = x + 2*y @@ -1289,11 +1374,11 @@ def correlation_matrix(nums_with_uncert): ufuncs = map(to_uncertain_func, nums_with_uncert) data = np.vstack([ufunc._mcpts for ufunc in ufuncs]) return np.corrcoef(data.T, rowvar=0) - + if __name__=='__main__': - + import mcerp.umath as umath - + print '*'*80 print 'TEST FUNCTIONS USING DERIVED MOMENTS FROM SCIPY DISTRIBUTIONS' print '*'*80 @@ -1303,7 +1388,7 @@ def correlation_matrix(nums_with_uncert): x3 = Exp(2) # Exp(mu=0.5) is the same Z = (x1*x2**2)/(15*(1.5 + x3)) Z.describe() - + print '*'*80 print 'Example of volumetric gas flow through orifice meter' H = N(64, 0.5) diff --git a/mcerp/unumpy.py b/mcerp/unumpy.py new file mode 100644 index 0000000..7a5544e --- /dev/null +++ b/mcerp/unumpy.py @@ -0,0 +1,155 @@ +""" +================================================================================ +mcerp: Real-time latin-hypercube-sampling-based Monte Carlo Error Propagation +================================================================================ + +Generalizes mathematical operators that work on numeric objects (from the math +module or numpy) compatible with objects with uncertainty distributions + +Author: Abraham Lee +Copyright: 2013 +""" +from mcerp import UncertainFunction, to_uncertain_func +import numpy as np +import umath as um + +__author__ = 'Tim Galvin' + +#!!!!!!!!!!!!!!!!!!!!!!!!!!! +# Conflicts with the abs keyword +#!!!!!!!!!!!!!!!!!!!!!!!!!!! + +# abs = np.vectorize(um.abs) + + + + +""" +Inverse cosine +""" +acos = np.vectorize(um.acos) + +""" +Inverse hyperbolic cosine +""" +acosh = np.vectorize(um.acosh) + +""" +Inverse sine +""" +asin = np.vectorize(um.asin) + +""" +Inverse hyperbolic sine +""" +asinh = np.vectorize(um.asinh) + +""" +Inverse tangent +""" +atan = np.vectorize(um.atan) + +""" +Inverse hyperbolic tangent +""" +atanh = np.vectorize(um.atanh) + +""" +Ceiling function (round towards positive infinity) +""" +ceil = np.vectorize(um.ceil) + +""" +Cosine +""" +cos = np.vectorize(um.cos) + +""" +Hyperbolic cosine +""" +cosh = np.vectorize(um.cosh) + +""" +Convert radians to degrees +""" +degrees = np.vectorize(um.degrees) + +""" +Exponential function +""" +exp = np.vectorize(um.exp) + +""" +Calculate exp(x) - 1 +""" +expm1 = np.vectorize(um.expm1) + +""" +Absolute value function +""" +fabs = np.vectorize(um.fabs) + +""" +Floor function (round towards negative infinity) +""" +floor = np.vectorize(um.floor) + +""" +Calculate the hypotenuse given two "legs" of a right triangle +""" +hypot = np.vectorize(um.hypot) + +""" +Natural logarithm (same as "log(x)") +""" +ln = np.vectorize(um.ln) + +""" +Natural logarithm +""" +log = np.vectorize(um.log) + +""" +Base-10 logarithm +""" +log10 = np.vectorize(um.log10) + +""" +Natural logarithm of (1 + x) +""" +log1p = np.vectorize(um.log1p) + +""" +Convert degrees to radians +""" +radians = np.vectorize(um.radians) + +""" +Sine +""" +sin = np.vectorize(um.sin) + +""" +Hyperbolic sine +""" +sinh = np.vectorize(um.sinh) + +""" +Square-root function +""" +sqrt = np.vectorize(um.sqrt) + +""" +Tangent +""" +tan = np.vectorize(um.tan) + +""" +Hyperbolic tangent +""" +tanh = np.vectorize(um.tanh) + +""" +Truncate the values to the integer value without rounding +""" +trunc = np.vectorize(um.trunc)