stat/distuv: more accurate calculation of Normal.CDF and LogNormal.CDF (

#580)

* A+C: add Taesu Pyo

* stat/distuv: more accurate calculation of Normal.CDF and Lognormal.CDF

Fixes #577

AUTHORS

@@ -64,6 +64,7 @@ source{d} <hello@sourced.tech>
 Shawn Smith <shawnpsmith@gmail.com>
 Spencer Lyon <spencerlyon2@gmail.com>
 Steve McCoy <mccoyst@gmail.com>
+Taesu Pyo <pyotaesu@gmail.com>
 Takeshi Yoneda <cz.rk.t0415y.g@gmail.com>
 The University of Adelaide
 The University of Minnesota

CONTRIBUTORS

@@ -70,6 +70,7 @@ Sebastien Binet <seb.binet@gmail.com>
 Shawn Smith <shawnpsmith@gmail.com>
 Spencer Lyon <spencerlyon2@gmail.com>
 Steve McCoy <mccoyst@gmail.com>
+Taesu Pyo <pyotaesu@gmail.com>
 Takeshi Yoneda <cz.rk.t0415y.g@gmail.com>
 Tobin Harding <me@tobin.cc>
 Vladimír Chalupecký <vladimir.chalupecky@gmail.com>

stat/distuv/lognormal.go

@@ -21,7 +21,7 @@ type LogNormal struct {
 
 // CDF computes the value of the cumulative density function at x.
 func (l LogNormal) CDF(x float64) float64 {
-	return 0.5 + 0.5*math.Erf((math.Log(x)-l.Mu)/(math.Sqrt2*l.Sigma))
+	return 0.5 * math.Erfc(-(math.Log(x)-l.Mu)/(math.Sqrt2*l.Sigma))
 }
 
 // Entropy returns the differential entropy of the distribution.

stat/distuv/lognormal_test.go

@@ -42,3 +42,16 @@ func TestLognormal(t *testing.T) {
 		checkProbQuantContinuous(t, i, x, f, tol)
 	}
 }
+
+// See https://github.com/gonum/gonum/issues/577 for details.
+func TestLognormalIssue577(t *testing.T) {
+	x := 1.0e-16
+	max := 1.0e-295
+	cdf := LogNormal{Mu: 0, Sigma: 1}.CDF(x)
+	if cdf <= 0 {
+		t.Errorf("LogNormal{0,1}.CDF(%e) should be positive. got: %e", x, cdf)
+	}
+	if cdf > max {
+		t.Errorf("LogNormal{0,1}.CDF(%e) is greater than %e. got: %e", x, max, cdf)
+	}
+}

stat/distuv/norm.go

@@ -29,7 +29,7 @@ type Normal struct {
 
 // CDF computes the value of the cumulative density function at x.
 func (n Normal) CDF(x float64) float64 {
-	return 0.5 * (1 + math.Erf((x-n.Mu)/(n.Sigma*math.Sqrt2)))
+	return 0.5 * math.Erfc(-(x-n.Mu)/(n.Sigma*math.Sqrt2))
 }
 
 // ConjugateUpdate updates the parameters of the distribution from the sufficient

stat/distuv/norm_test.go

@@ -169,3 +169,16 @@ func BenchmarkNormalQuantile(b *testing.B) {
 		}
 	}
 }
+
+// See https://github.com/gonum/gonum/issues/577 for details.
+func TestNormalIssue577(t *testing.T) {
+	x := -36.0
+	max := 1.e-282
+	cdf := Normal{Mu: 0, Sigma: 1}.CDF(x)
+	if cdf <= 0 {
+		t.Errorf("Normal{0,1}.CDF(%e) should be positive. got: %e", x, cdf)
+	}
+	if cdf > max {
+		t.Errorf("Normal{0,1}.CDF(%e) is greater than %e. got: %e", x, max, cdf)
+	}
+}