diff --git a/AUTHORS b/AUTHORS index 6a90f815..19c64ab5 100644 --- a/AUTHORS +++ b/AUTHORS @@ -9,6 +9,7 @@ # Please keep the list sorted. Brendan Tracey +Bill Gray Bill Noon Chih-Wei Chang Chris Tessum diff --git a/CONTRIBUTORS b/CONTRIBUTORS index 3b563e7a..1ffaad97 100644 --- a/CONTRIBUTORS +++ b/CONTRIBUTORS @@ -16,6 +16,7 @@ # Please keep the list sorted. Brendan Tracey +Bill Gray Bill Noon Chih-Wei Chang Chris Tessum diff --git a/stat/distuv/triangle.go b/stat/distuv/triangle.go new file mode 100644 index 00000000..157b9ea1 --- /dev/null +++ b/stat/distuv/triangle.go @@ -0,0 +1,192 @@ +// 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" +) + +// Triangle represents a triangle distribution (https://en.wikipedia.org/wiki/Triangular_distribution). +type Triangle struct { + a, b, c float64 + Source *rand.Rand +} + +// NewTriangle constructs a new triangle distribution with lower limit a, upper limit b, and mode c. +// Constraints are a < b and a ≤ c ≤ b. +// This distribution is uncommon in nature, but may be useful for simulation. +func NewTriangle(a, b, c float64) Triangle { + checkTriangleParameters(a, b, c) + return Triangle{a, b, c, nil} +} + +func checkTriangleParameters(a, b, c float64) { + if a >= b { + panic("triangle: constraint of a < b violated") + } + if a > c { + panic("triangle: constraint of a <= c violated") + } + if c > b { + panic("triangle: constraint of c <= b violated") + } +} + +// CDF computes the value of the cumulative density function at x. +func (t Triangle) CDF(x float64) float64 { + switch { + case x <= t.a: + return 0 + case x <= t.c: + d := x - t.a + return (d * d) / ((t.b - t.a) * (t.c - t.a)) + case x < t.b: + d := t.b - x + return 1 - (d*d)/((t.b-t.a)*(t.b-t.c)) + default: + return 1 + } +} + +// Entropy returns the entropy of the distribution. +func (t Triangle) Entropy() float64 { + return 0.5 + math.Log(t.b-t.a) - math.Ln2 +} + +// ExKurtosis returns the excess kurtosis of the distribution. +func (Triangle) ExKurtosis() float64 { + return -3.0 / 5.0 +} + +// Fit is not appropriate for Triangle, because the distribution is generally used when there is little data. + +// LogProb computes the natural logarithm of the value of the probability density function at x. +func (t Triangle) LogProb(x float64) float64 { + return math.Log(t.Prob(x)) +} + +// Mean returns the mean of the probability distribution. +func (t Triangle) Mean() float64 { + return (t.a + t.b + t.c) / 3 +} + +// Median returns the median of the probability distribution. +func (t Triangle) Median() float64 { + if t.c >= (t.a+t.b)/2 { + return t.a + math.Sqrt((t.b-t.a)*(t.c-t.a)/2) + } + return t.b - math.Sqrt((t.b-t.a)*(t.b-t.c)/2) +} + +// Mode returns the mode of the probability distribution. +func (t Triangle) Mode() float64 { + return t.c +} + +// NumParameters returns the number of parameters in the distribution. +func (Triangle) NumParameters() int { + return 3 +} + +// Prob computes the value of the probability density function at x. +func (t Triangle) Prob(x float64) float64 { + switch { + case x < t.a: + return 0 + case x < t.c: + return 2 * (x - t.a) / ((t.b - t.a) * (t.c - t.a)) + case x == t.c: + return 2 / (t.b - t.a) + case x <= t.b: + return 2 * (t.b - x) / ((t.b - t.a) * (t.b - t.c)) + default: + return 0 + } +} + +// Quantile returns the inverse of the cumulative probability distribution. +func (t Triangle) Quantile(p float64) float64 { + if p < 0 || p > 1 { + panic(badPercentile) + } + + f := (t.c - t.a) / (t.b - t.a) + + if p < f { + return t.a + math.Sqrt(p*(t.b-t.a)*(t.c-t.a)) + } + return t.b - math.Sqrt((1-p)*(t.b-t.a)*(t.b-t.c)) +} + +// Rand returns a random sample drawn from the distribution. +func (t Triangle) Rand() float64 { + var rnd float64 + if t.Source == nil { + rnd = rand.Float64() + } else { + rnd = t.Source.Float64() + } + + return t.Quantile(rnd) +} + +// Skewness returns the skewness of the distribution. +func (t Triangle) Skewness() float64 { + n := math.Sqrt2 * (t.a + t.b - 2*t.c) * (2*t.a - t.b - t.c) * (t.a - 2*t.b + t.c) + d := 5 * math.Pow(t.a*t.a+t.b*t.b+t.c*t.c-t.a*t.b-t.a*t.c-t.b*t.c, 3.0/2.0) + + return n / d +} + +// StdDev returns the standard deviation of the probability distribution. +func (t Triangle) StdDev() float64 { + return math.Sqrt(t.Variance()) +} + +// Survival returns the survival function (complementary CDF) at x. +func (t Triangle) Survival(x float64) float64 { + return 1 - t.CDF(x) +} + +// MarshalParameters implements the ParameterMarshaler interface +func (t Triangle) MarshalParameters(p []Parameter) { + if len(p) != t.NumParameters() { + panic("triangle: improper parameter length") + } + p[0].Name = "A" + p[0].Value = t.a + p[1].Name = "B" + p[1].Value = t.b + p[2].Name = "C" + p[2].Value = t.c +} + +// UnmarshalParameters implements the ParameterMarshaler interface +func (t *Triangle) UnmarshalParameters(p []Parameter) { + if len(p) != t.NumParameters() { + panic("triangle: incorrect number of parameters to set") + } + if p[0].Name != "A" { + panic("triangle: " + panicNameMismatch) + } + if p[1].Name != "B" { + panic("triangle: " + panicNameMismatch) + } + if p[2].Name != "C" { + panic("triangle: " + panicNameMismatch) + } + + checkTriangleParameters(p[0].Value, p[1].Value, p[2].Value) + + t.a = p[0].Value + t.b = p[1].Value + t.c = p[2].Value +} + +// Variance returns the variance of the probability distribution. +func (t Triangle) Variance() float64 { + return (t.a*t.a + t.b*t.b + t.c*t.c - t.a*t.b - t.a*t.c - t.b*t.c) / 18 +} diff --git a/stat/distuv/triangle_test.go b/stat/distuv/triangle_test.go new file mode 100644 index 00000000..a11152e1 --- /dev/null +++ b/stat/distuv/triangle_test.go @@ -0,0 +1,83 @@ +// 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" + "testing" +) + +func TestTriangleConstraint(t *testing.T) { + defer func() { + if r := recover(); r == nil { + t.Errorf("The constraints were violated, but not caught") + } + }() + + // test b < a + NewTriangle(3, 1, 2) + // test c > b + NewTriangle(1, 2, 3) +} + +func TestTriangle(t *testing.T) { + for i, test := range []struct { + a, b, c float64 + }{ + { + a: 0.0, + b: 1.0, + c: 0.5, + }, + { + a: 0.1, + b: 0.3, + c: 0.2, + }, + { + a: 1.0, + b: 2.0, + c: 1.5, + }, + { + a: 0.0, + b: 1.0, + c: 0.0, + }, + } { + dist := NewTriangle(test.a, test.b, test.c) + testFullDist(t, dist, i, true) + } +} + +func TestTriangleProb(t *testing.T) { + pts := []univariateProbPoint{ + { + loc: 0.5, + prob: 0, + cumProb: 0, + logProb: math.Inf(-1), + }, + { + loc: 1, + prob: 0, + cumProb: 0, + logProb: math.Inf(-1), + }, + { + loc: 2, + prob: 1.0, + cumProb: 0.5, + logProb: 0, + }, + { + loc: 3, + prob: 0, + cumProb: 1, + logProb: math.Inf(-1), + }, + } + testDistributionProbs(t, NewTriangle(1, 3, 2), "Standard 1,2,3 Triangle", pts) +}