Added student's t-distribution

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
+}

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)
+      }
+    }
+  }
+}