Moved Gaussian to stats package and improved tests

Author: fasmat

stats/cdf.go

@@ -2,10 +2,30 @@ package stats
 
 import "math"
 
-// NormalCDF calculates the Cumulative Distribution Function for the
-// standardized normal distribution
+// NormGaussAt returns the probability of x in a standardized normal
+// distribution
+func NormGaussAt(x float64) float64 {
+	return GaussAt(x, 0, 1)
+}
+
+// GaussAt returns the probability of x for the normal distribution with mean mu
+// and standard deviation sigma
+func GaussAt(x, mu, sigma float64) float64 {
+	multiplier := 1.0 / (sigma * math.Sqrt(2*math.Pi))
+	expPart := math.Exp((-1.0 * math.Pow(x-mu, 2.0)) / (2 * sigma * sigma))
+	return multiplier * expPart
+}
+
+// NormalCDF calculates the Cumulative Distribution function for the
+// standardized normal distribution.
 func NormalCDF(x float64) float64 {
-	return (1.0 / 2.0) * (1 + math.Erf(x/math.Sqrt2))
+	return CDF(x, 0, 1)
+}
+
+// CDF calculates the cumulative distribution function for any standard
+// distribution.
+func CDF(x, mu, sigma float64) float64 {
+	return 0.5 * (1 + math.Erf((x-mu)/(sigma*math.Sqrt2)))
 }
 
 // InverseCDF calculates the Inverse Cumulative Distribution Function for the
@@ -15,14 +35,18 @@ func NormalCDF(x float64) float64 {
 func InverseCDF(p float64) float64 {
 	var r, x, pp, dp float64
 
-	dp = p - 0.5
+	if p > 1 || p < 0 {
+		return math.NaN()
+	}
+
 	switch p {
 	case 1.0:
 		return math.Inf(1)
 	case 0.0:
 		return math.Inf(-1)
 	}
 
+	dp = p - 0.5
 	if math.Abs(dp) <= 0.425 {
 		x = small(dp)
 		return x

stats/cdf_test.go

@@ -5,22 +5,83 @@ import (
 	"testing"
 )
 
+const (
+	ErrTolerance = 1e-6
+)
+
+func TestNormGaussAt(t *testing.T) {
+	in := []float64{
+		0,
+		1,
+		-1,
+		3.7,
+		-3.7,
+	}
+
+	// These values were generated using dnorm on http://www.r-fiddle.org/
+	out := []float64{
+		0.3989422804,
+		0.24197072452,
+		0.24197072452,
+		0.00042478027055,
+		0.00042478027055,
+	}
+
+	for i, v := range in {
+		result := NormGaussAt(v)
+		if math.Abs(result-out[i]) > ErrTolerance {
+			t.Error("Expected", out[i], "got", result, "for", v)
+		}
+	}
+}
+
+func TestGaussAt(t *testing.T) {
+	in := []float64{
+		20,
+		25,
+		15,
+		38.5,
+		1.5,
+	}
+
+	// These values were generated using dnorm on http://www.r-fiddle.org/
+	out := []float64{
+		0.07978845608,
+		0.048394144904,
+		0.048394144904,
+		8.495605411e-5,
+		8.495605411e-5,
+	}
+
+	for i, v := range in {
+		result := GaussAt(v, 20, 5)
+		if math.Abs(result-out[i]) > ErrTolerance {
+			t.Error("Expected", out[i], "got", result, "for", v)
+		}
+	}
+}
+
 func TestNormalCDF(t *testing.T) {
 	in := []float64{
 		math.Inf(1),
 		0,
 		math.Inf(-1),
-		InverseCDF(0.1),
-		InverseCDF(0.9),
+		1,
+		-1,
+		2,
+		-2,
 		InverseCDF(0.86),
 	}
 
+	// These results were calculated using pnorm(x) on www.r-fiddle.org
 	out := []float64{
 		1,
 		0.5,
 		0,
-		0.1,
-		0.9,
+		0.841344746069,
+		0.158655253931,
+		0.977249868052,
+		0.022750131948,
 		0.86,
 	}
 
@@ -31,7 +92,7 @@ func TestNormalCDF(t *testing.T) {
 	for i := range in {
 		o := NormalCDF(in[i])
 
-		if math.Abs(o-out[i]) > errTolerance {
+		if math.Abs(o-out[i]) > ErrTolerance {
 			t.Errorf("NormalCDF(%f) == %f, want %f", in[i], o, out[i])
 		}
 	}
@@ -42,17 +103,24 @@ func TestInverseCDF(t *testing.T) {
 		1,
 		0.5,
 		0,
-		NormalCDF(1),
-		NormalCDF(-1),
+		0.2,
+		0.3,
+		0.9,
+		-1,
+		2,
 		NormalCDF(4),
 	}
 
+	// These results were calculated using qnorm(x) on www.r-fiddle.org
 	out := []float64{
 		math.Inf(1),
 		0,
 		math.Inf(-1),
-		1,
-		-1,
+		-0.84162123357,
+		-0.52440051271,
+		1.2815515655,
+		math.NaN(),
+		math.NaN(),
 		4,
 	}
 
@@ -63,7 +131,7 @@ func TestInverseCDF(t *testing.T) {
 	for i := range in {
 		o := InverseCDF(in[i])
 
-		if math.Abs(o-out[i]) > errTolerance {
+		if math.Abs(o-out[i]) > ErrTolerance {
 			t.Errorf("InverseCDF(%f) == %f, want %f", in[i], o, out[i])
 		}
 	}

gauss.go stats/gauss.go

@@ -1,4 +1,4 @@
-package trueskill
+package stats
 
 import (
 	"fmt"
@@ -61,7 +61,7 @@ func (g *Gaussian) GetConservativeEstimate() float64 {
 
 // String can be used to print a Gaussian
 func (g *Gaussian) String() string {
-	text := fmt.Sprintf("μ=%1.4f, σ=%1.4f", g.GetMu(), g.GetSigma())
+	text := fmt.Sprintf("μ=%2.4f, σ=%2.4f", g.GetMu(), g.GetSigma())
 	return text
 }