NormDist added in statistics

parent 713629c3
......@@ -96,6 +96,7 @@ SOURCES = \
language/python/module/_sha3_512.scm \
language/python/module/_shake_128.scm \
language/python/module/_shake_256.scm \
language/python/module/hashlib.py \
language/python/module/csv.py \
language/python/module/datetime.py \
language/python/module/textwrap.py \
......
......@@ -81,13 +81,14 @@ A single exception is defined: StatisticsError is a subclass of ValueError.
__all__ = [ 'StatisticsError','geometric_mean','multimode',
'pstdev', 'pvariance', 'stdev', 'variance',
'median', 'median_low', 'median_high', 'median_grouped',
'mean', 'fmean', 'mode', 'harmonic_mean',
'mean', 'fmean', 'mode', 'harmonic_mean', 'NormalDist'
]
import collections
import decimal
import math
import numbers
import random
from fractions import Fraction
from decimal import Decimal
......@@ -732,7 +733,11 @@ def stdev(data, xbar=None):
except AttributeError:
return math.sqrt(var)
f = 1/math.sqrt(2.0*math.pi)
def z(x):
return f * math.exp(-x*x/2.0)
def pstdev(data, mu=None):
"""Return the square root of the population variance.
......@@ -747,3 +752,178 @@ def pstdev(data, mu=None):
return var.sqrt()
except AttributeError:
return math.sqrt(var)
def invz(x):
if x < 0:
sgn = -1
else:
sgn = 1
x = (1-x)*(1+x)
lnx = math.log(x)
tt1 = 2/(math.pi*0.147) + +.5*lnx
tt2 = 1 / (0.147) * lnx
return sgn*math.sqrt(tt1*tt1 - tt2)
class NormalDist:
def __init__(self,m,s):
self._m = m
self._ss = s*s
def __repr__(self):
return 'NormDist(m=%f,s=%f)'%(self._m,math.sqrt(self._ss))
def __add__(self,o):
if isinstance(o,NormalDist):
m = self._m + o._m
ss = self._ss + o._ss
return NormalDist(m,math.sqrt(ss))
else:
m = self._m + o
ss = self._ss
return NormalDist(m,math.sqrt(ss))
def __radd__(self,o):
if isinstance(o,NormalDist):
m = self._m + o._m
ss = self._ss + o._ss
return NormalDist(m,math.sqrt(ss))
else:
m = self._m + o
ss = self._ss
return NormalDist(m,math.sqrt(ss))
def __mul__(self,o):
m = self._m * o
ss = self._ss * o * o
return NormalDist(m,math.sqrt(ss))
def __rmul__(self,o):
m = self._m * o
ss = self._ss * o * o
return NormalDist(m,math.sqrt(ss))
def __sub__(self,o):
if isinstance(o,NormalDist):
m = self._m - o._m
ss = self._ss + o._ss
return NormalDist(m,math.sqrt(ss))
else:
m = self._m - o
ss = self._ss
return NormalDist(m,math.sqrt(ss))
def __neg__(self):
m = -self._m
ss = self._ss
return NormalDist(m,math.sqrt(ss))
@classmethod
def from_samples(cls,data):
n = len(data)
s = 0
ss = 0
for x in data:
s += x
ss += x*x
m = s / n
ss = (ss - s*s/n)/(n-1)
return cls(m,math.sqrt(ss))
def _get_mean(self):
return self._m
def _set_mean(self):
raise AttributeError('NormDist mean attribute read only')
def _get_std(self):
return math.sqrt(self._ss)
def _set_std(self):
raise AttributeError('NormDist std attribute read only')
def _get_var(self):
return self._ss
def _set_var(self):
raise AttributeError('NormDist var attribute read only')
mean = property(_get_mean,_set_mean)
stddev = property(_get_std,_set_std)
variance = property(_get_var,_set_var)
def samples(self,n,*,seed=None):
if not seed==None:
random.seed(seed)
ret = []
for i in range(n):
ret.append(random.normalvariate(self._m, math.sqrt(self._ss)))
return ret
def pdf(self,x):
s = math.sqrt(self._ss)
return z((x-self._m)/s) / s
def cdf(self,x):
if x < self._m:
ret = 0.5 - math.erf((self._m - x) / math.sqrt(self._ss))
else:
ret = 0.5 + math.erf((x - self._m) / math.sqrt(self._ss))
return ret;
def inv_cdf(self,x):
return math.sqrt(self._ss) * (invz(x) + self._m)
def overlap(self,o):
ss1 = self._ss
ss2 = o._ss
m1 = self._m
m2 = self._m
K = math.log(ss2 / ss1)
R = ss2 - ss1
b = (ss1*m2 - ss2*m1)/R
c = (ss2*m1*m1 - ss1*m1*m1 - K*ss1*ss2)/R
D = math.sqrt(b*b - c)
x0 = -b
x1 = x0 - D
x2 = x0 + D
if x2 < x1:
t = x2
x2 = x1
x1 = t
e11 = -(x2 - m1)/ss1*self.pdf(x2)
e12 = -(x2 - m2)/ss2*o.pdf(x2)
e21 = -(x2 - m1)/ss1*self.pdf(x2)
e22 = -(x2 - m2)/ss2*o.pdf(x2)
ret = 0
if e11 < e12:
ret += o.cdf(x1)
ret += self.cdf(x2) - self.cdf(x1)
else:
ret += self.cdf(x1)
ret += o.cdf(x2) - o.cdf(x1)
if e21 < e22:
ret += 1 - self.cdf(x2)
else:
ret += 1 - o.cdf(x2)
return ret
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment