overlap statistic debugged

parent 67578fe1
...@@ -771,8 +771,8 @@ def invz(x): ...@@ -771,8 +771,8 @@ def invz(x):
class NormalDist: class NormalDist:
def __init__(self,m,s): def __init__(self,m,s):
self._m = m self._m = float(m)
self._ss = s*s self._ss = float(s*s)
def __repr__(self): def __repr__(self):
return 'NormDist(m=%f,s=%f)'%(self._m,math.sqrt(self._ss)) return 'NormDist(m=%f,s=%f)'%(self._m,math.sqrt(self._ss))
...@@ -875,9 +875,9 @@ class NormalDist: ...@@ -875,9 +875,9 @@ class NormalDist:
def cdf(self,x): def cdf(self,x):
if x < self._m: if x < self._m:
ret = 0.5 - math.erf((self._m - x) / math.sqrt(self._ss)) ret = 0.5 - math.erf((self._m - x) / math.sqrt(self._ss))/2.0
else: else:
ret = 0.5 + math.erf((x - self._m) / math.sqrt(self._ss)) ret = 0.5 + math.erf((x - self._m) / math.sqrt(self._ss))/2.0
return ret; return ret;
...@@ -888,12 +888,27 @@ class NormalDist: ...@@ -888,12 +888,27 @@ class NormalDist:
ss1 = self._ss ss1 = self._ss
ss2 = o._ss ss2 = o._ss
m1 = self._m m1 = self._m
m2 = self._m m2 = o._m
if (abs(ss1 - ss2) / math.sqrt(ss1 + ss2)) < 1e-6:
b = 2*(ss1*m2 - ss2*m1)
c = (ss2*m1*m1 - ss1*m2*m2)
if (abs(b)/math.sqrt(ss1+ss2)) < 1e-6 :
return 1
else:
x0 = -c / b
e1 = -(x0 - m1)/ss1*self.pdf(x0)
e2 = -(x0 - m2)/ss2*o.pdf(x0)
if e1 > e2:
return self.cdf(x0) + 1 - o.cdf(x0)
else:
return o.cdf(x0) + 1 - self.cdf(x0)
K = math.log(ss2 / ss1) K = math.log(ss2 / ss1)
R = ss2 - ss1 R = ss2 - ss1
b = (ss1*m2 - ss2*m1)/R b = (ss1*m2 - ss2*m1)/R
c = (ss2*m1*m1 - ss1*m1*m1 - K*ss1*ss2)/R c = (ss2*m1*m1 - ss1*m2*m2 - K*ss1*ss2)/R
D = math.sqrt(b*b - c) D = math.sqrt(b*b - c)
...@@ -906,8 +921,8 @@ class NormalDist: ...@@ -906,8 +921,8 @@ class NormalDist:
x2 = x1 x2 = x1
x1 = t x1 = t
e11 = -(x2 - m1)/ss1*self.pdf(x2) e11 = -(x1 - m1)/ss1*self.pdf(x1)
e12 = -(x2 - m2)/ss2*o.pdf(x2) e12 = -(x1 - m2)/ss2*o.pdf(x1)
e21 = -(x2 - m1)/ss1*self.pdf(x2) e21 = -(x2 - m1)/ss1*self.pdf(x2)
e22 = -(x2 - m2)/ss2*o.pdf(x2) e22 = -(x2 - m2)/ss2*o.pdf(x2)
......
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