diff --git a/AUTHORS b/AUTHORS index 4df34329..27434b9b 100644 --- a/AUTHORS +++ b/AUTHORS @@ -38,6 +38,7 @@ Joseph Watson Josh Wilson Julien Roland Kent English +Kevin C. Zimmerman Konstantin Shaposhnikov Leonid Kneller Lyron Winderbaum diff --git a/CONTRIBUTORS b/CONTRIBUTORS index 90cd0764..949cac81 100644 --- a/CONTRIBUTORS +++ b/CONTRIBUTORS @@ -45,6 +45,7 @@ Joseph Watson Josh Wilson Julien Roland Kent English +Kevin C. Zimmerman Konstantin Shaposhnikov Leonid Kneller Lyron Winderbaum diff --git a/stat/distuv/poisson.go b/stat/distuv/poisson.go new file mode 100644 index 00000000..2ae859e2 --- /dev/null +++ b/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 +} diff --git a/stat/distuv/poisson_test.go b/stat/distuv/poisson_test.go new file mode 100644 index 00000000..3fa918e9 --- /dev/null +++ b/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) +}