feat: add students t-distribution

Author: sritchie

src/gen/distribution/commons_math.clj

@@ -11,6 +11,7 @@
             EnumeratedIntegerDistribution
             GammaDistribution
             NormalDistribution
+            TDistribution
             UniformIntegerDistribution
             UniformRealDistribution)
            (org.apache.commons.math3.random
@@ -36,6 +37,19 @@
   (sample [^AbstractIntegerDistribution obj]
     (.sample obj)))
 
+;; Small wrapper around a kixi T-distribution to allow for location and scale
+;; parameters.
+
+(defrecord LocationScaleT [^TDistribution t-dist location scale]
+  d/LogPDF
+  (logpdf [_ v]
+    (- (.logDensity t-dist (/ (- v location) scale))
+       (Math/log scale)))
+
+  d/Sample
+  (sample [_]
+    (+ location (* scale (.sample t-dist)))))
+
 ;; ## Primitive probability distributions
 
 (defn ^:no-doc rng ^RandomGenerator []
@@ -62,6 +76,13 @@
 (defn gamma-distribution [^double shape ^double scale]
   (GammaDistribution. (rng) shape scale))
 
+(defn student-t-distribution
+  ([^double nu]
+   (TDistribution. (rng) nu))
+  ([nu location scale]
+   (->LocationScaleT
+    (student-t-distribution nu) location scale)))
+
 (defn normal-distribution
   ([] (normal-distribution 0.0 1.0))
   ([^double mean ^double sd]
@@ -92,6 +113,9 @@
 (def gamma
   (d/->GenerativeFn gamma-distribution))
 
+(def student-t
+  (d/->GenerativeFn student-t-distribution))
+
 (def normal
   (d/->GenerativeFn normal-distribution))
 

src/gen/distribution/kixi.cljc

@@ -5,7 +5,7 @@
   #?(:clj
      (:import (kixi.stats.distribution Bernoulli Cauchy
                                        Exponential Beta
-                                       Gamma Normal Uniform))))
+                                       Gamma Normal Uniform T))))
 
 ;; ## Kixi.stats protocol implementations
 ;;
@@ -80,6 +80,26 @@
                  (.-sd this)
                  v)))
 
+(extend-type #?(:clj T :cljs k/T)
+  d/Sample
+  (sample [this] (k/draw this))
+
+  d/LogPDF
+  (logpdf [this v]
+    (ll/student-t (.-dof this) 0 1 v)))
+
+;; Small wrapper around a kixi T-distribution to allow for location and scale
+;; parameters.
+
+(defrecord LocationScaleT [^T t-dist location scale]
+  d/Sample
+  (sample [_]
+    (+ location (* scale (k/draw t-dist))))
+
+  d/LogPDF
+  (logpdf [_ v]
+    (ll/student-t (.-dof t-dist) location scale v)))
+
 ;; ## Primitive probability distributions
 
 (defn bernoulli-distribution
@@ -110,6 +130,14 @@
 (defn gamma-distribution [shape scale]
   (k/gamma {:shape shape :scale scale}))
 
+(defn student-t-distribution
+  ([nu]
+   (k/t {:v nu}))
+  ([nu location scale]
+   (->LocationScaleT (student-t-distribution nu)
+                     location
+                     scale)))
+
 ;; ## Primitive generative functions
 
 (def bernoulli
@@ -132,3 +160,6 @@
 
 (def gamma
   (d/->GenerativeFn gamma-distribution))
+
+(def student-t
+  (d/->GenerativeFn student-t-distribution))

src/gen/distribution/math/log_likelihood.cljc

@@ -167,7 +167,6 @@
                (Math/log variance)
                (/ v-mu-sq variance)))))
 
-
 (defn uniform
   "Returns the log-likelihood of the continuous [uniform
   distribution](https://en.wikipedia.org/wiki/Continuous_uniform_distribution)
@@ -177,3 +176,26 @@
   (if (<= a v b)
     (- (Math/log (- b a)))
     ##-Inf))
+
+(defn student-t
+  "Returns the log-likelihood of the [non-standardized Student's
+  t-distribution](https://en.wikipedia.org/wiki/Student's_t-distribution#Location-scale_transformation)
+  parametrized by `location`, `scale` and degrees-of-freedom `nu` at the value
+  `v`.
+
+  This distribution is also known as the location-scale t-distribution.
+
+  The implementation follows the algorithm described on the
+  distribution's [Wikipedia
+  page](https://en.wikipedia.org/wiki/Student's_t-distribution#Location-scale_transformation)."
+  ([nu v] (student-t nu 0 1 v))
+  ([nu location scale v]
+   (let [inc-nu      (inc nu)
+         half-inc-nu (* 0.5 inc-nu)
+         normalized  (/ (- v location) scale)
+         norm**2     (* normalized normalized)]
+     (- (log-gamma-fn half-inc-nu)
+        (log-gamma-fn (* 0.5 nu))
+        (Math/log scale)
+        (* 0.5 (+ log-pi (Math/log nu)))
+        (* half-inc-nu (Math/log (inc (/ norm**2 nu))))))))

test/gen/distribution/commons_math_test.clj

@@ -18,3 +18,6 @@
 
 (deftest gamma-tests
   (dt/gamma-tests commons/gamma-distribution))
+
+(deftest student-t-tests
+  (dt/student-t-tests commons/student-t-distribution))

test/gen/distribution/kixi_test.cljc

@@ -24,3 +24,6 @@
 
 (deftest gamma-tests
   (dt/gamma-tests kixi/gamma-distribution))
+
+(deftest student-t-tests
+  (dt/student-t-tests kixi/student-t-distribution))

test/gen/distribution/math/log_likelihood_test.cljc

@@ -4,16 +4,9 @@
             [gen.distribution.math.log-likelihood :as ll]
             [gen.distribution :as distribution]
             [gen.distribution-test :as dt]
-            [gen.test-check-util :refer [gen-double]]
+            [gen.test-check-util :refer [gen-double within]]
             [same.core :refer [ish? with-comparator]]))
 
-(defn within
-  "Returns a function that tests whether two values are within `eps` of each
-  other."
-  [^double eps]
-  (fn [^double x ^double y]
-    (< (Math/abs (- x y)) eps)))
-
 (defn factorial
   "Factorial implementation for testing."
   [n]
@@ -79,3 +72,6 @@
 
 (deftest uniform-tests
   (dt/uniform-tests (->logpdf ll/uniform)))
+
+(deftest student-t-tests
+  (dt/student-t-tests (->logpdf ll/student-t)))

test/gen/distribution_test.cljc

@@ -7,14 +7,27 @@
             [gen.dynamic.choice-map :as choice-map]
             [gen.generative-function :as gf]
             [gen.trace :as trace]
-            [gen.test-check-util :refer [gen-double]]
-            [same.core :refer [ish? zeroish?]]))
+            [gen.test-check-util :refer [gen-double within]]
+            [same.core :refer [ish? zeroish? with-comparator]]))
 
 (defn gamma-tests [->gamma]
   (testing "spot checks"
     (is (= -6.391804444241573 (dist/logpdf (->gamma 0.001 1) 0.4)))
     (is (= -393.0922447210179 (dist/logpdf (->gamma 1 0.001) 0.4)))))
 
+(defn student-t-tests [->student-t]
+  (testing "spot checks"
+    (with-comparator (within 1e-12)
+      (is (ish? -1.7347417805005154 (dist/logpdf (->student-t 2 2.1 2) 2)))
+      (is (ish? -2.795309741614719 (dist/logpdf (->student-t 1 0.8 4) 3)))))
+
+  (checking "Student's T matches generalized logpdf"
+            [v  (gen-double -10 10)
+             nu (gen/fmap inc gen/nat)]
+            (is (= (dist/logpdf (->student-t nu 0 1) v)
+                   (dist/logpdf (->student-t nu) v))
+                "these two paths should produce the same results")))
+
 (defn beta-tests [->beta]
   (testing "spot checks"
     (is (= -5.992380837839856 (dist/logpdf (->beta 0.001 1) 0.4)))
@@ -161,9 +174,10 @@
                       (dist/logpdf (->normal 0.0 sigma) (- v)))
                 "Normal is symmetric about the mean")
 
-            (is (ish? (dist/logpdf (->normal mu sigma) v)
-                      (dist/logpdf (->normal (+ mu shift) sigma) (+ v shift)))
-                "shifting by the mean is a symmetry"))
+            (with-comparator (within 1e-12)
+              (is (ish? (dist/logpdf (->normal mu sigma) v)
+                        (dist/logpdf (->normal (+ mu shift) sigma) (+ v shift)))
+                  "shifting by the mean is a symmetry")))
 
   (testing "spot checks"
     (is (= -1.0439385332046727 (dist/logpdf (->normal 0 1) 0.5)))

test/gen/test_check_util.cljc

@@ -2,6 +2,13 @@
   "Utilities for using Gen distributions and types with test.check."
   (:require [clojure.test.check.generators :as gen]))
 
+(defn within
+  "Returns a function that tests whether two values are within `eps` of each
+  other."
+  [^double eps]
+  (fn [^double x ^double y]
+    (< (Math/abs (- x y)) eps)))
+
 (defn gen-double
   "Returns a generator that produces numerical doubles between `min` and
   `max` (inclusive)."