Add poisson primitive support

src/inferenceql/inference/gpm.cljc

@@ -12,6 +12,7 @@
             [inferenceql.inference.gpm.primitive-gpms.bernoulli :as bernoulli]
             [inferenceql.inference.gpm.primitive-gpms.categorical :as categorical]
             [inferenceql.inference.gpm.primitive-gpms.gaussian :as gaussian]
+            [inferenceql.inference.gpm.primitive-gpms.poisson :as poisson]
             [inferenceql.inference.gpm.proto :as gpm-proto]
             [inferenceql.inference.gpm.view :as view]))
 
@@ -152,6 +153,7 @@
    'inferenceql.inference.gpm.primitive_gpms.bernoulli.Bernoulli bernoulli/map->Bernoulli
    'inferenceql.inference.gpm.primitive_gpms.categorical.Categorical categorical/map->Categorical
    'inferenceql.inference.gpm.primitive_gpms.gaussian.Gaussian gaussian/map->Gaussian
+   'inferenceql.inference.gpm.primitive_gpms.poisson.Poisson poisson/map->Poisson
    'inferenceql.inference.gpm.view.View view/map->View})
 
 (defn as-gpm

src/inferenceql/inference/gpm/primitive_gpms.cljc

@@ -1,15 +1,17 @@
 (ns inferenceql.inference.gpm.primitive-gpms
   (:require [inferenceql.inference.gpm.primitive-gpms.bernoulli :as bernoulli]
             [inferenceql.inference.gpm.primitive-gpms.categorical :as categorical]
-            [inferenceql.inference.gpm.primitive-gpms.gaussian :as gaussian]))
+            [inferenceql.inference.gpm.primitive-gpms.gaussian :as gaussian]
+            [inferenceql.inference.gpm.primitive-gpms.poisson :as poisson]))
 
 (defn primitive?
   "Checks whether the given GPM is a primitive GPM."
   [stattype]
   (and (record? stattype)
        (or (bernoulli/bernoulli? stattype)
            (categorical/categorical? stattype)
-           (gaussian/gaussian? stattype))))
+           (gaussian/gaussian? stattype)
+           (poisson/poisson? stattype))))
 
 (defn hyper-grid
   [stattype data & {:keys [n-grid] :or {n-grid 30}}]
@@ -19,6 +21,7 @@
       :bernoulli (bernoulli/hyper-grid data n-grid)
       :categorical (categorical/hyper-grid data n-grid)
       :gaussian (gaussian/hyper-grid data n-grid)
+      :poisson (poisson/hyper-grid data n-grid)
       (throw (ex-info (str "pGPM doesn't exist: " stattype)
                       {:stattype stattype :data data})))))
 
@@ -39,5 +42,6 @@
     :bernoulli (bernoulli/spec->bernoulli var-name :suff-stats suff-stats :hyperparameters hyperparameters)
     :categorical (categorical/spec->categorical var-name :suff-stats suff-stats :hyperparameters hyperparameters :options options)
     :gaussian (gaussian/spec->gaussian var-name :suff-stats suff-stats :hyperparameters hyperparameters)
+    :poisson (poisson/spec->poisson var-name :suff-stats suff-stats :hyperparameters hyperparameters)
     (throw (ex-info (str "pGPM doesn't exist for var-name: " primitive " for " var-name)
                     {:primitive primitive}))))

src/inferenceql/inference/gpm/primitive_gpms/poisson.cljc

@@ -0,0 +1,77 @@
+(ns inferenceql.inference.gpm.primitive-gpms.poisson
+  (:require [clojure.math :as math]
+            [inferenceql.inference.gpm.proto :as gpm.proto]
+            [inferenceql.inference.primitives :as primitives]
+            [inferenceql.inference.distributions :as dist]
+            [inferenceql.inference.utils :as utils]))
+
+(defn posterior-hypers
+ [n sum-x  a  b]
+  [(+ a sum-x) (+ b n)])
+
+(defn calc-log-Z
+  [a b]
+  (- (dist/log-gamma a) (* a (math/log b))))
+
+(defrecord Poisson [var-name suff-stats hyperparameters]
+  gpm.proto/GPM
+  (logpdf [_ targets constraints]
+    (let [x (get targets var-name)
+          x' (get constraints var-name)
+          constrained? (not (nil? x'))]
+      (cond
+        (nil? x) 0
+        constrained? (if (= x x') 0 ##-Inf)
+        :else (let [n (:n suff-stats)
+                    sum-x  (:sum-x suff-stats)
+                    sum-log-fact (:sum-log-fact suff-stats)
+                    a  (:a hyperparameters)
+                    b  (:b hyperparameters)
+                    [an bn] (posterior-hypers n  sum-x  a  b)
+                    [am bm] (posterior-hypers (+ n 1)  (+ sum-x x)  a  b)
+                    Zn (calc-log-Z an bn)
+                    Zm (calc-log-Z am bm)]
+          (- Zm Zn (dist/log-gamma (+ x 1)))))))
+
+  (simulate [this _ _]
+    (throw (Exception. "Poisson simulate not implemented")))
+
+  gpm.proto/Incorporate
+  (incorporate [this values]
+    (let [x (get values var-name)]
+      (assoc this :suff-stats (-> suff-stats
+                                  (update :n inc)
+                                  (update :sum-x #(+ % x))
+                                  (update :sum-log-fact #(+ % (dist/log-gamma (+ x 1))))))))
+  (unincorporate [this values]
+    (let [x (get values var-name)]
+      (assoc this :suff-stats (-> suff-stats
+                                  (update :n dec)
+                                  (update :sum-x #(- % x))
+                                  (update :sum-log-fact #(- % (dist/log-gamma (+ x 1))))))))
+
+  gpm.proto/Variables
+  (variables [{:keys [var-name]}]
+    #{var-name}))
+
+(defn poisson?
+  "Checks if the given pGPM is Poisson."
+  [stattype]
+  (and (record? stattype)
+       (instance? Poisson stattype)))
+
+(defn hyper-grid
+  "Hyperparameter grid for the Poisson variable, used in column hyperparameter inference
+  for Column GPMs."
+  [data n-grid]
+  (let [grid (utils/log-linspace 1 (count data) n-grid)]
+    {:alpha grid}))
+
+(defn spec->poisson
+  "Casts a CrossCat category spec to a Poisson pGPM.
+  Requires a variable name, optionally takes by key
+  sufficient statistics, options, and hyperparameters."
+  [var-name & {:keys [hyperparameters suff-stats options]}]
+  (let [suff-stats' (if-not (nil? suff-stats) suff-stats {:n 0 :sum-x 0 :sum-log-fact 0})
+        hyperparameters' (if-not (nil? hyperparameters) hyperparameters {:a 1 :b 1})]
+    (->Poisson var-name suff-stats' hyperparameters')))

test/inferenceql/inference/gpm/primitive_gpms/poisson_test.cljc

@@ -0,0 +1,22 @@
+(ns inferenceql.inference.gpm.primitive-gpms.poisson-test
+  (:require [clojure.math :as math]
+            [clojure.test :as test :refer [deftest is]]
+            [inferenceql.inference.gpm :as gpm]
+            [inferenceql.inference.gpm.primitive-gpms.poisson :as poisson]
+            [inferenceql.inference.gpm.proto :as gpm.proto]
+            [inferenceql.inference.utils :as utils]))
+
+(def var-name "poisson")
+
+(def poisson-pgpm
+  (let [suff-stats {:n 0 :sum-x 0 :sum-log-fact 0}
+        hypers {:a 2 :b 2}]
+    (poisson/spec->poisson var-name :suff-stats suff-stats :hyperparameters hypers)))
+
+(deftest logpdf
+  (let [targets {"poisson" 0}
+        constraints {"poisson" 1}]
+    (is (= 1.0 (math/exp (gpm.proto/logpdf poisson-pgpm {} {}))))
+    (is (= 1.0 (math/exp (gpm.proto/logpdf poisson-pgpm targets targets))))
+    (is (= ##-Inf (gpm.proto/logpdf poisson-pgpm targets constraints))))
+    (is (utils/almost-equal? -1.2163953243244932 (gpm.proto/logpdf poisson-pgpm {"poisson" 1} {}) utils/relerr 1E-8)))