Ping, or, Bayesian random bits

[lukego]

Here is a random post of stuff that I have in my Emacs buffers today.

Cute type declarations:

(deftype R   (&rest args) "Real number."    `(double-float ,@args))
(deftype R[] ()           "Real vector."    '(simple-array R (*)))
(deftype P   ()           "Probability."    '(R 0d0 1d0))
(deftype L   ()           "Log-likelihood." '(R * #.(log 1)))

Example of DEFMODEL macro to define statistical models for Bayesian parameter inference and model selection:

(defmodel line (x y &param m c σ)
  "Linear relationship between X and Y with Gaussian noise of constant scale:
    y = m*x + c + N(0,σ)
   Infers parameters M (gradient), C (intercept), and σ (standard deviation.)"
  (gaussian-log-likelihood (+ c (* m x)) σ y))

Macroexpansion of the above into a bespoke Sequential Monte Carlo simulation:

(defun line (&key n-particles observations (jitter-scales '(0.01 0.1 0.5)))
  "Linear relationship between X and Y with Gaussian noise of constant scale:
    y = m*x + c + N(0,σ)
   Infers parameters M (gradient), C (intercept), and σ (standard deviation.)"
  (let ((#:m (make-array (list n-particles) :element-type 'r))
        (#:c (make-array (list n-particles) :element-type 'r))
        (#:σ (make-array (list n-particles) :element-type 'r)))
    (labels ((log-likelihood (m c σ x y)
               (gaussian-log-likelihood (+ c (* m x)) σ y))
             (particle-log-likelihood (i x y)
               (log-likelihood (aref #:m i) (aref #:c i) (aref #:σ i) x y))
             (respawn! (parents)
               (reorder! parents #:m #:c #:σ))
             (jitter! (metropolis-accept?)
               (loop for stddev in jitter-scales
                     do (loop for i below n-particles
                              for m = (aref #:m i)
                              for c = (aref #:c i)
                              for σ = (aref #:σ i)
                              for ll.old = (partial #'log-likelihood m c σ)
                              for #:m.p = (add-noise m stddev)
                              for #:c.p = (add-noise c stddev)
                              for #:σ.p = (add-noise σ stddev)
                              for ll.new = (partial #'log-likelihood #:m.p
                                                    #:c.p #:σ.p)
                              when (funcall metropolis-accept? ll.old ll.new)
                              do (setf (aref #:m i) #:m.p
                                       (aref #:c i) #:c.p
                                       (aref #:σ i) #:σ.p))))
             (add-noise (x stddev)
               (+ x (* stddev (gaussian-random)))))
      (smc/likelihood-tempering n-particles observations :log-likelihood
       #'particle-log-likelihood :respawn! #'respawn! :jitter! #'jitter!))))

I am having fun 👍