DEFMODEL line expansion

gistfile1.txt

(defun line
       (observations
        &key (n-particles 100) (steps 100)
        (jitter-scales '(0.01d0 0.35d0 1.0d0)))
  "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 (init-parameter-vector n-particles -10.0d0 10.0d0))
        (#:c (init-parameter-vector n-particles -10.0d0 10.0d0))
        (#:σ (init-parameter-vector n-particles 1.1102230246251568d-16 5.0d0)))
    (labels ((log-likelihood (m c σ x y)
               (declare (type r m c σ x y))
               (declare (inline gaussian-log-likelihood))
               (gaussian-log-likelihood (+ c (* m x)) σ y))
             (particle-log-likelihood (i x y)
               (declare (type r 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 jitter 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 = (+ m
                                             (* (gaussian-random) jitter
                                                (/ 2.2222222222222223d0
                                                   (expt n-particles 1/3))))
                              for #:c.p = (+ c
                                             (* (gaussian-random) jitter
                                                (/ 2.2222222222222223d0
                                                   (expt n-particles 1/3))))
                              for #:σ.p = (+ σ
                                             (* (gaussian-random) jitter
                                                (/ 0.5555555555555556d0
                                                   (expt n-particles 1/3))))
                              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)))))
      (values
       (smc/likelihood-tempering n-particles observations :log-likelihood
                                 #'particle-log-likelihood :respawn! #'respawn!
                                 :jitter! #'jitter! :temp-step (/ 1.0d0 steps))
       (dict 'm #:m 'c #:c 'σ #:σ)))))