implement Poisson based on PR #259 (#304)

* stat/distuv: implement Poisson

AUTHORS

@@ -38,6 +38,7 @@ Joseph Watson <jtwatson@linux-consulting.us>
 Josh Wilson <josh.craig.wilson@gmail.com>
 Julien Roland <juroland@gmail.com>
 Kent English <kent.english@gmail.com>
+Kevin C. Zimmerman <kevinczimmerman@gmail.com>
 Konstantin Shaposhnikov <k.shaposhnikov@gmail.com>
 Leonid Kneller <recondite.matter@gmail.com>
 Lyron Winderbaum <lyron.winderbaum@student.adelaide.edu.au>

CONTRIBUTORS

@@ -45,6 +45,7 @@ Joseph Watson <jtwatson@linux-consulting.us>
 Josh Wilson <josh.craig.wilson@gmail.com>
 Julien Roland <juroland@gmail.com>
 Kent English <kent.english@gmail.com>
+Kevin C. Zimmerman <kevinczimmerman@gmail.com>
 Konstantin Shaposhnikov <k.shaposhnikov@gmail.com>
 Leonid Kneller <recondite.matter@gmail.com>
 Lyron Winderbaum <lyron.winderbaum@student.adelaide.edu.au>

stat/distuv/poisson.go

@@ -0,0 +1,109 @@
+// Copyright ©2017 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package distuv
+
+import (
+	"math"
+	"math/rand"
+
+	"gonum.org/v1/gonum/mathext"
+)
+
+// Poisson implements the Poisson distribution, a discrete probability distribution
+// that expresses the probability of a given number of events occurring in a fixed
+// interval.
+// The poisson distribution has density function:
+//  f(k) = λ^k / k! e^(-λ)
+// For more information, see https://en.wikipedia.org/wiki/Poisson_distribution.
+type Poisson struct {
+	// Lambda is the average number of events in an interval.
+	// Lambda must be greater than 0.
+	Lambda float64
+
+	Source *rand.Rand
+}
+
+// CDF computes the value of the cumulative distribution function at x.
+func (p Poisson) CDF(x float64) float64 {
+	if x < 0 {
+		return 0
+	}
+	return mathext.GammaIncComp(math.Floor(x+1), p.Lambda)
+}
+
+// ExKurtosis returns the excess kurtosis of the distribution.
+func (p Poisson) ExKurtosis() float64 {
+	return 1 / p.Lambda
+}
+
+// LogProb computes the natural logarithm of the value of the probability
+// density function at x.
+func (p Poisson) LogProb(x float64) float64 {
+	if x < 0 || math.Floor(x) != x {
+		return math.Inf(-1)
+	}
+	lg, _ := math.Lgamma(math.Floor(x) + 1)
+	return x*math.Log(p.Lambda) - p.Lambda - lg
+}
+
+// Mean returns the mean of the probability distribution.
+func (p Poisson) Mean() float64 {
+	return p.Lambda
+}
+
+// NumParameters returns the number of parameters in the distribution.
+func (Poisson) NumParameters() int {
+	return 1
+}
+
+// Prob computes the value of the probability density function at x.
+func (p Poisson) Prob(x float64) float64 {
+	return math.Exp(p.LogProb(x))
+}
+
+// Rand returns a random sample drawn from the distribution.
+func (p Poisson) Rand() float64 {
+	// poisson generator based upon the multiplication of
+	// uniform random variates.
+	// see:
+	//  Non-Uniform Random Variate Generation,
+	//  Luc Devroye (p504)
+	//  http://www.eirene.de/Devroye.pdf
+	x := 0.0
+	prod := 1.0
+	exp := math.Exp(-p.Lambda)
+	rnd := rand.Float64
+	if p.Source != nil {
+		rnd = p.Source.Float64
+	}
+
+	for {
+		prod *= rnd()
+		if prod <= exp {
+			return x
+		}
+		x++
+	}
+}
+
+// Skewness returns the skewness of the distribution.
+func (p Poisson) Skewness() float64 {
+	return 1 / math.Sqrt(p.Lambda)
+}
+
+// StdDev returns the standard deviation of the probability distribution.
+func (p Poisson) StdDev() float64 {
+	return math.Sqrt(p.Variance())
+}
+
+// Survival returns the survival function (complementary CDF) at x.
+func (p Poisson) Survival(x float64) float64 {
+	return 1 - p.CDF(x)
+}
+
+// Variance returns the variance of the probability distribution.
+func (p Poisson) Variance() float64 {
+	return p.Lambda
+}

stat/distuv/poisson_test.go

@@ -0,0 +1,121 @@
+// Copyright ©2017 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package distuv
+
+import (
+	"math/rand"
+	"sort"
+	"testing"
+
+	"gonum.org/v1/gonum/floats"
+)
+
+func TestPoissonProb(t *testing.T) {
+	const tol = 1e-10
+	for i, tt := range []struct {
+		k      float64
+		lambda float64
+		want   float64
+	}{
+		{0, 1, 3.678794411714423e-01},
+		{1, 1, 3.678794411714423e-01},
+		{2, 1, 1.839397205857211e-01},
+		{3, 1, 6.131324019524039e-02},
+		{4, 1, 1.532831004881010e-02},
+		{5, 1, 3.065662009762020e-03},
+		{6, 1, 5.109436682936698e-04},
+		{7, 1, 7.299195261338139e-05},
+		{8, 1, 9.123994076672672e-06},
+		{9, 1, 1.013777119630298e-06},
+
+		{0, 2.5, 8.208499862389880e-02},
+		{1, 2.5, 2.052124965597470e-01},
+		{2, 2.5, 2.565156206996838e-01},
+		{3, 2.5, 2.137630172497365e-01},
+		{4, 2.5, 1.336018857810853e-01},
+		{5, 2.5, 6.680094289054267e-02},
+		{6, 2.5, 2.783372620439277e-02},
+		{7, 2.5, 9.940616501568845e-03},
+		{8, 2.5, 3.106442656740263e-03},
+		{9, 2.5, 8.629007379834082e-04},
+
+		{0.5, 2.5, 0},
+		{1.5, 2.5, 0},
+		{2.5, 2.5, 0},
+		{3.5, 2.5, 0},
+		{4.5, 2.5, 0},
+		{5.5, 2.5, 0},
+		{6.5, 2.5, 0},
+		{7.5, 2.5, 0},
+		{8.5, 2.5, 0},
+		{9.5, 2.5, 0},
+	} {
+		p := Poisson{Lambda: tt.lambda}
+		got := p.Prob(tt.k)
+		if !floats.EqualWithinAbs(got, tt.want, tol) {
+			t.Errorf("test-%d: got=%e. want=%e\n", i, got, tt.want)
+		}
+	}
+}
+
+func TestPoissonCDF(t *testing.T) {
+	const tol = 1e-10
+	for i, tt := range []struct {
+		k      float64
+		lambda float64
+		want   float64
+	}{
+		{0, 1, 0.367879441171442},
+		{1, 1, 0.735758882342885},
+		{2, 1, 0.919698602928606},
+		{3, 1, 0.981011843123846},
+		{4, 1, 0.996340153172656},
+		{5, 1, 0.999405815182418},
+		{6, 1, 0.999916758850712},
+		{7, 1, 0.999989750803325},
+		{8, 1, 0.999998874797402},
+		{9, 1, 0.999999888574522},
+
+		{0, 2.5, 0.082084998623899},
+		{1, 2.5, 0.287297495183646},
+		{2, 2.5, 0.543813115883329},
+		{3, 2.5, 0.757576133133066},
+		{4, 2.5, 0.891178018914151},
+		{5, 2.5, 0.957978961804694},
+		{6, 2.5, 0.985812688009087},
+		{7, 2.5, 0.995753304510655},
+		{8, 2.5, 0.998859747167396},
+		{9, 2.5, 0.999722647905379},
+	} {
+		p := Poisson{Lambda: tt.lambda}
+		got := p.CDF(tt.k)
+		if !floats.EqualWithinAbs(got, tt.want, tol) {
+			t.Errorf("test-%d: got=%e. want=%e\n", i, got, tt.want)
+		}
+	}
+}
+
+func TestPoisson(t *testing.T) {
+	src := rand.New(rand.NewSource(1))
+	for i, b := range []Poisson{
+		{3, src},
+		{1.5, src},
+		{0.9, src},
+	} {
+		testPoisson(t, b, i)
+	}
+}
+
+func testPoisson(t *testing.T, p Poisson, i int) {
+	tol := 1e-2
+	const n = 1e6
+	x := make([]float64, n)
+	generateSamples(x, p)
+	sort.Float64s(x)
+
+	checkMean(t, i, x, p, tol)
+	checkVarAndStd(t, i, x, p, tol)
+	checkExKurtosis(t, i, x, p, 7e-2)
+}