diff --git a/studentst.go b/studentst.go new file mode 100755 index 0000000..d472d0d --- /dev/null +++ b/studentst.go @@ -0,0 +1,192 @@ +package distributions + +import ( + "math" +) + +//The Student's t-Distribution is a continuous probability distribution +// with parameters df > 0. +// +// See: https://en.wikipedia.org/wiki/Student's_t-distribution +type StudentsT struct { + Degrees float64 +} + +func (dist StudentsT) validate() error { + if dist.Degrees <= 0 { + return InvalidParamsError{ "Degrees must be greater than zero." } + } + return nil +} + +func (dist StudentsT) Mean() (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + if (dist.Degrees <= 1) { + return math.NaN(), nil + } + return 0.0, nil +} + +func (dist StudentsT) Variance() (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + if (dist.Degrees < 1) { + return math.NaN(), nil + } + if (dist.Degrees <= 2) { + return math.Inf(1), nil + } + result := dist.Degrees / (dist.Degrees - 2) + return result, nil +} + +func (dist StudentsT) Skewness() (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + if (dist.Degrees <= 3) { + return math.NaN(), nil + } + return 0.0, nil +} + +func (dist StudentsT) Kurtosis() (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + if (dist.Degrees <= 4) { + return math.NaN(), nil + } + result := 3 * (dist.Degrees - 2) / (dist.Degrees - 4) + return result, nil +} + +func (dist StudentsT) StdDev() (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + if (dist.Degrees < 1) { + return math.NaN(), nil + } + if (dist.Degrees <= 2) { + return math.Inf(1), nil + } + result := math.Sqrt(dist.Degrees / (dist.Degrees - 2)) + return result, nil +} + +func (dist StudentsT) RelStdDev() (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + return math.NaN(), nil +} + +func (dist StudentsT) Pdf(x float64) (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + lg1, _ := math.Lgamma(dist.Degrees / 2) + lg2, _ := math.Lgamma((dist.Degrees + 1) / 2) + result := math.Exp(lg2 - lg1) * math.Pow(1 + (x * x / dist.Degrees), -(dist.Degrees + 1) / 2) / math.Sqrt(math.Pi * dist.Degrees) + return result, nil +} + +// Ref: https://github.com/chbrown/nlp/blob/master/src/main/java/cc/mallet/util/StatFunctions.java +func (dist StudentsT) Cdf(x float64) (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + g1 := 1.0 / math.Pi + idf := dist.Degrees + a := x / math.Sqrt(idf) + b := idf / (idf + x * x) + im2 := dist.Degrees - 2.0 + ioe := math.Mod(idf, 2.0) + s := 1.0 + c := 1.0 + idf = 1.0 + ks := 2.0 + ioe + fk := ks + if im2 >= 2.0 { + for k := ks; k <= im2; k += 2.0 { + c = c * b * (fk - 1.0) / fk + s += c + if s != idf { + idf = s + fk += 2.0 + } + } + } + if ioe != 1 { + result := 0.5 + (0.5 * a * math.Sqrt(b) * s) + return result, nil + } + if dist.Degrees == 1 { + s = 0 + } + result := 0.5 + (((a * b * s) + math.Atan(a)) * g1) + return result, nil +} + +// Ref: https://github.com/ampl/gsl/blob/master/randist/tdist.c +func (dist StudentsT) random() (float64, error) { + if err := dist.validate(); err != nil { + return math.NaN(), err + } + if (dist.Degrees <= 2) { + normal := Normal{ Mu: 0, Sigma: 1} + chiSquared := ChiSquared{ Degrees: dist.Degrees } + y1, _, err := normal.random() + if err != nil { + return math.NaN(), err + } + y2, err := chiSquared.random() + if err != nil { + return math.NaN(), err + } + result := y1 / math.Sqrt(y2 / dist.Degrees) + return result, nil + } else { + var y1, y2, z float64 + var err error + normal := Normal{ Mu: 0, Sigma: 1} + exponential := Exponential{ Lambda: 1 / ((dist.Degrees / 2) - 1) } + ok := true + for ok { + y1, _, err = normal.random() + if err != nil { + return math.NaN(), err + } + y1, err = exponential.random() + if err != nil { + return math.NaN(), err + } + z = y1 * y2 / (dist.Degrees - 2) + ok = 1 - z < 0 || math.Exp(-y2 - z) > 1 - z + } + result := y1 / math.Sqrt((1 - (2 / dist.Degrees)) * (1 - z)) + return result, nil + } +} + +func (dist StudentsT) Sample(n int) ([]float64, error) { + if err := dist.validate(); err != nil { + return []float64{}, err + } + if n <= 0 { + return []float64{}, nil + } + result := make([]float64, n) + for i := 0; i < n; i++ { + value, err := dist.random() + if err != nil { + return []float64{}, nil + } + result[i] = value + } + return result, nil +} diff --git a/studentst_test.go b/studentst_test.go new file mode 100755 index 0000000..d2b315c --- /dev/null +++ b/studentst_test.go @@ -0,0 +1,128 @@ +package distributions + +import ( + "math" + "testing" +) + +type studentstTest struct { + dist Distribution + mean float64 + variance float64 + stdDev float64 + relStdDev float64 + skewness float64 + kurtosis float64 + pdf []inOut + cdf []inOut +} + +// Test at http://keisan.casio.com/exec/system/1180573203 +func Test_StudentsT(t *testing.T) { + examples := []studentstTest{ + studentstTest{ + dist: StudentsT{10}, + mean: 0.0, + variance: 1.25, + stdDev: 1.118033988749895, + relStdDev: math.NaN(), + skewness: 0.0, + kurtosis: 4.0, + pdf: []inOut{ + inOut{ in: 9.0, out: 0.0000020670116801089978 }, + inOut{ in: 2.5, out: 0.0269387276282444589776 }, + inOut{ in: 4.0, out: 0.00203103391104121595875 }, + }, + cdf: []inOut{ + inOut{ in: 9.0, out: 0.9999979309754751589939 }, + inOut{ in: 2.5, out: 0.9842765778816955978753 }, + inOut{ in: 4.0, out: 0.9987408336876316538681 }, + }, + }, + // This is a Exponential distribution ;P + studentstTest{ + dist: StudentsT{2}, + mean: 0.0, + variance: math.Inf(1), + stdDev: math.Inf(1), + relStdDev: math.NaN(), + skewness: math.NaN(), + kurtosis: math.NaN(), + pdf: []inOut{ + inOut{ in: 9.0, out: 0.00132246096373120901175 }, + inOut{ in: 2.5, out: 0.0422006438680479607703 }, + inOut{ in: 4.0, out: 0.0130945700219731023037 }, + }, + cdf: []inOut{ + inOut{ in: 9.0, out: 0.9939391699536065658886 }, + inOut{ in: 2.5, out: 0.935194139889244595443 }, + inOut{ in: 4.0, out: 0.9714045207910316829339 }, + }, + }, + } + + for _, example := range examples { + mean, err := example.dist.Mean() + if err != nil || !floatsPicoEqual(mean, example.mean) { + if !checkInf(mean, example.mean) && !checkNaN(mean, example.mean) { + t.Fatalf("\nMean:\n Expected: %f\n Got: %f\n", example.mean, mean) + } + } + variance, err := example.dist.Variance() + if err != nil || !floatsPicoEqual(variance, example.variance) { + if !checkInf(variance, example.variance) && !checkNaN(variance, example.variance) { + t.Fatalf("\nVariance:\n Expected: %f\n Got: %f\n", example.variance, variance) + } + } + stdDev, err := example.dist.StdDev() + if err != nil || !floatsPicoEqual(stdDev, example.stdDev) { + if !checkInf(stdDev, example.stdDev) && !checkNaN(stdDev, example.stdDev) { + t.Fatalf("\nStdDev:\n Expected: %f\n Got: %f\n", example.stdDev, stdDev) + } + } + relStdDev, err := example.dist.RelStdDev() + if err != nil || !floatsPicoEqual(relStdDev, example.relStdDev) { + if !checkInf(relStdDev, example.relStdDev) && !checkNaN(relStdDev, example.relStdDev) { + t.Fatalf("\nRelStdDev:\n Expected: %f\n Got: %f\n", example.relStdDev, relStdDev) + } + } + skewness, err := example.dist.Skewness() + if err != nil || !floatsPicoEqual(skewness, example.skewness) { + if !checkInf(skewness, example.skewness) && !checkNaN(skewness, example.skewness) { + t.Fatalf("\nSkewness:\n Expected: %f\n Got: %f\n", example.skewness, skewness) + } + } + kurtosis, err := example.dist.Kurtosis() + if err != nil || !floatsPicoEqual(kurtosis, example.kurtosis) { + if !checkInf(kurtosis, example.kurtosis) && !checkNaN(kurtosis, example.kurtosis) { + t.Fatalf("\nKurtosis:\n Expected: %f\n Got: %f\n", example.kurtosis, kurtosis) + } + } + for _, pdf := range example.pdf { + out, err := example.dist.Pdf(pdf.in) + if err != nil || !floatsPicoEqual(out, pdf.out) { + t.Fatalf("\nPdf of %f:\n Expected: %f\n Got: %f\n", pdf.in, pdf.out, out) + } + } + for _, cdf := range example.cdf { + out, err := example.dist.Cdf(cdf.in) + if err != nil || !floatsPicoEqual(out, cdf.out) { + t.Fatalf("\nCdf of %f:\n Expected: %f\n Got: %f\n", cdf.in, cdf.out, out) + } + } + samples, err := example.dist.Sample(1000000) + if err != nil { + t.Fatalf("\nCould not generate 1,000,000 samples.") + } + sampleMean := averageFloats(samples) + if !floatsIntegerEqual(example.mean, sampleMean) { + t.Fatalf("\nSample average:\n Expected: %f\n Got: %f\n", example.mean, sampleMean) + } + if !math.IsInf(example.variance,0) { + sampleVar := varianceFloats(samples, sampleMean) + if !floatsEqual(example.variance, sampleVar, 1.5) { + t.Fatalf("\nSample variance:\n Expected: %f\n Got: %f\n", example.variance, sampleVar) + } + } + } +}