Added beta distribution

beta.go

@@ -0,0 +1,213 @@
+package distributions
+
+import (
+  "math"
+  "math/rand"
+)
+
+//TheBeta  Distribution is a continuous probability distribution
+// with parameters α > 0, β >= 0.
+//
+// See: https://en.wikipedia.org/wiki/Beta_distribution
+type Beta struct {
+  Alpha   float64
+  Beta    float64
+}
+
+func (dist Beta) validate() error {
+  if dist.Alpha <= 0 {
+    return InvalidParamsError{ "Alpha must be greater than zero." }
+  }
+  if dist.Beta <= 0 {
+    return InvalidParamsError{ "Beta must be greater than zero." }
+  }
+  return nil
+}
+
+func (dist Beta) Mean() (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  result := dist.Alpha / (dist.Alpha + dist.Beta)
+  return result, nil
+}
+
+func (dist Beta) Variance() (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  comb := dist.Alpha + dist.Beta
+  result := (dist.Alpha * dist.Beta) / (comb * comb * (comb + 1))
+  return result, nil
+}
+
+func (dist Beta) Skewness() (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  numer := 2 * (dist.Beta - dist.Alpha) * math.Sqrt(dist.Alpha + dist.Beta + 1)
+  result := numer / ((dist.Alpha + dist.Beta + 2) * math.Sqrt(dist.Alpha * dist.Beta))
+  return result, nil
+}
+
+func (dist Beta) Kurtosis() (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  first := (dist.Alpha - dist.Beta) * (dist.Alpha - dist.Beta) * (dist.Alpha + dist.Beta + 1)
+  second := dist.Alpha * dist.Beta * (dist.Alpha + dist.Beta + 2)
+  denom := dist.Alpha * dist.Beta * (dist.Alpha + dist.Beta + 2) * (dist.Alpha + dist.Beta + 3)
+  result := 6 * (first - second) / denom
+  return result, nil
+}
+
+func (dist Beta) StdDev() (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  variance, _ := dist.Variance()
+  result := math.Sqrt(variance)
+  return result, nil
+}
+
+func (dist Beta) RelStdDev() (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  mean, _ := dist.Mean()
+  stdDev, _ := dist.StdDev()
+  result := stdDev / mean
+  return result, nil
+}
+
+func (dist Beta) Pdf(x float64) (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  beta := BetaFn(dist.Alpha, dist.Beta)
+  result := math.Pow(x, dist.Alpha - 1) * math.Pow(1 - x, dist.Beta - 1) / beta
+  return result, nil
+}
+
+func (dist Beta) Cdf(x float64) (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  result := RegBetaInc(dist.Alpha, dist.Beta, x)
+  return result, nil
+}
+
+// Ref: https://compbio.soe.ucsc.edu/gen_sequence/gen_beta.c
+// Ref: https://github.com/e-dard/godist/blob/master/beta.go
+func (dist Beta) random() (float64, error) {
+  if err := dist.validate(); err != nil {
+    return math.NaN(), err
+  }
+  min := math.Min(dist.Alpha, dist.Beta)
+  max := math.Max(dist.Alpha, dist.Beta)
+  // Use Joehnk's algorithm.
+  if max < 0.5 {
+    u1 := rand.Float64()
+    u2 := rand.Float64()
+    result := math.Pow(u1, 1 / dist.Alpha) / (math.Pow(u1, 1 / dist.Alpha) + math.Pow(u2, 1 / dist.Beta))
+    return result, nil
+  // Use Cheng's BC Algorithm
+  } else if min <= 1.0 {
+    var u1, u2, v, w, y, z float64
+    alpha := min + max
+    beta := 1.0 / min
+    delta := 1.0 + max - min
+  	k1 := delta * (0.0138889 + (0.0416667 * min)) / ((max * beta) - 0.777778)
+  	k2 := 0.25 + ((0.5 + 0.25) * min / delta)
+    setvw := func() {
+      v = beta * math.Log(u1 / (1.0 - u1))
+      if v <= 709.78 {
+        w = dist.Alpha * math.Exp(v)
+        if math.IsInf(w,0) {
+          w = math.MaxFloat64
+        }
+      } else {
+        w = math.MaxFloat64
+      }
+    }
+    for {
+      u1 = rand.Float64()
+      u2 = rand.Float64()
+      if u1 < 0.5 {
+        y = u1 * u2
+        z = u1 * y
+        if (0.25 * u2) + z - y >= k1 {
+          continue
+        }
+      } else {
+        z = u1 * u1 * u2
+        if z <= 0.25 {
+          setvw()
+          break
+        }
+        if z >= k2 {
+          continue
+        }
+      }
+      setvw()
+      if alpha * (math.Log(alpha / (min + w)) + v) - 1.3862944 >= math.Log(z) {
+        break
+      }
+    }
+    if dist.Alpha == min {
+      return min / (min + w), nil
+    }
+    return w / (min + w), nil
+  // Use Cheng's BB Algorithm
+  } else {
+    alpha := min + max
+    beta := math.Sqrt((alpha - 2.0) / ((2.0 * min * max) - alpha))
+    gamma := min + (1.0 / beta)
+    var r, s, t, v, w, z float64
+    for {
+      u1 := rand.Float64()
+      u2 := rand.Float64()
+      v = beta * math.Log(u1 / (1.0 - u1))
+      if v <= 709.78 {
+        w = dist.Alpha * math.Exp(v)
+        if math.IsInf(w,0) {
+          w = math.MaxFloat64
+        }
+      } else {
+        w = math.MaxFloat64
+      }
+      z = u1 * u1 * u2
+      r = (gamma * v) - 1.3862944
+      s = min + r - w
+      if s + 2.609438 >= 5.0 * z {
+        break
+      }
+      t = math.Log(z)
+      if r + (alpha * math.Log(alpha / (max + w))) < t {
+        break
+      }
+    }
+    if dist.Alpha != min {
+      return max / (max + w), nil
+    }
+    return w / (max + w), nil
+  }
+}
+
+func (dist Beta) 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
+}

beta_test.go

@@ -0,0 +1,126 @@
+package distributions
+
+import (
+  "testing"
+)
+
+type betaTest 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/1180573225
+func Test_Beta(t *testing.T) {
+  examples := []betaTest{
+    betaTest{
+      dist:       Beta{1, 2},
+      mean:       1.0 / 3.0,
+      variance:   1.0 / 18.0,
+      stdDev:     0.2357022603955158414669,
+      relStdDev:  0.7071067811865475244008,
+      skewness:   0.5656854249492380195207,
+      kurtosis:   -3.0 / 5.0,
+      pdf: []inOut{
+        inOut{ in: 0.4,   out: 1.2 },
+        inOut{ in: 0.6,   out: 0.8 },
+        inOut{ in: 0.14,  out: 1.72 },
+
+      },
+      cdf: []inOut{
+        inOut{ in: 0.4,   out: 0.64 },
+        inOut{ in: 0.6,   out: 0.84 },
+        inOut{ in: 0.14,  out: 0.2604 },
+      },
+    },
+    betaTest{
+      dist:       Beta{5, 4},
+      mean:       5.0 / 9.0,
+      variance:   2.0 / 81.0,
+      stdDev:     0.1571348402636772276446,
+      relStdDev:  0.258198889747161125678,
+      skewness:   -0.1285648693066450044365,
+      kurtosis:   -21.0 / 44.0,
+      pdf: []inOut{
+        inOut{ in: 0.4,   out: 1.548288 },
+        inOut{ in: 0.6,   out: 2.322432 },
+        inOut{ in: 0.14,  out: 0.0684172364288 },
+
+      },
+      cdf: []inOut{
+        inOut{ in: 0.4,   out: 0.1736704 },
+        inOut{ in: 0.6,   out: 0.5940864 },
+        inOut{ in: 0.14,  out: 0.002079010303104 },
+      },
+    },
+  }
+
+  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)
+    }
+    sampleVar := varianceFloats(samples, sampleMean)
+    if !floatsIntegerEqual(example.variance, sampleVar) {
+      t.Fatalf("\nSample variance:\n  Expected: %f\n  Got: %f\n", example.variance, sampleVar)
+    }
+  }
+}