gmm_1d model does not compile

[haohaiziround]

The GMM model does not compile. Here is the full model.

(defn sample-components [_ zs pi]
  (let [z (sample (categorical pi))]
    (conj zs z)))

(defn observe-data [n _ ys zs mus]
  (let [y (get ys n)
        z (get zs n)
        mu (get mus z)]
    (observe (normal mu 2) y)
    nil))

(let [ys      (vector -2.0  -2.5  -1.7  -1.9  -2.2
                      1.5  2.2  3  1.2  2.8)
      pi [0.5 0.5]
      zs  (loop 10 (vector) sample-components pi)
      mus (vector (sample (normal 0 2))
                  (sample (normal 0 2)))]
  (loop 10 nil observe-data ys zs mus)
  (vector mus zs))

The model is about clustering N data (ys) into K groups. Here N=10 and K=2. We have RVs of mus with length K and zs with length N. We first sample K prior cluster centers mus. Then for each data ys[n] where n = 1, ... , N, we sample its assignment zs[n] = k where k is an integer between 1 and K, and observe with corresponding mus[k].

[haohaiziround]

The loop syntax is a little bit tricky. It is defined as (loop n_iter val f param_1, ..., param_n). And the function f is defined as (defn f [n val param_1, ..., param_n] ...), where n is the n-th iteration and val is the return value for each and also final iteration.

[Tobias-Kohn]

I concur: loop is a ghastly beast, primarily meant as a substitute for recursion. You could try and use map instead. We should also add for to the language.
I think, we should tackle such comple models in the new implementation PyFOPPL-2. Due to its type system, it should be better equipped to handle such cases.

[bayesianbrad]

I concur. As of Monday morning, we shall switch to pyfoppl 2 and resolve all issues simultaneously.

…

On 20 Jan 2018 23:39, "Tobias Kohn" ***@***.***> wrote: I concur: loop is a ghastly beast, primarily meant as a substitute for recursion. You could try and use map instead. We should also add for to the language. I think, we should tackle such comple models in the new implementation PyFOPPL-2. Due to its type system, it should be better equipped to handle such cases. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#13 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/Ab2RWr-CNhTULG8-gpTbBASMs318RsUNks5tMnk5gaJpZM4RlkQC> .

[bayesianbrad]

Also there is a small typo in the model it should be the following (notice that (vector sample...) has changed to (vector) sample...

In Brookes' FOPPL the current compiler output is as follows:

| Vertices V: #{y30793 x30676 y30775 x30701 y30739 y30766 y30784 x30706 x30671 y30811 y30757 y30748 y30730 y30802 x30719 x30686 x30666 x30691 x30716 x30681 x30696 x30711}
-- | --

Arcs A: #{}
--

Conditional densities P:
--
y30793 -> (fn [] (normal nil 2))
x30676 -> (fn [] (categorical [0.5 0.5]))
y30775 -> (fn [] (normal nil 2))
x30701 -> (fn [] (categorical [0.5 0.5]))
y30739 -> (fn [] (normal nil 2))
y30766 -> (fn [] (normal nil 2))
y30784 -> (fn [] (normal nil 2))
x30706 -> (fn [] (categorical [0.5 0.5]))
x30671 -> (fn [] (categorical [0.5 0.5]))
y30811 -> (fn [] (normal nil 2))
y30757 -> (fn [] (normal nil 2))
y30748 -> (fn [] (normal nil 2))
y30730 -> (fn [] (normal nil 2))
y30802 -> (fn [] (normal nil 2))
x30719 -> (fn [] (normal 0 2))
x30686 -> (fn [] (categorical [0.5 0.5]))
x30666 -> (fn [] (categorical [0.5 0.5]))
x30691 -> (fn [] (categorical [0.5 0.5]))
x30716 -> (fn [] (normal 0 2))
x30681 -> (fn [] (categorical [0.5 0.5]))
x30696 -> (fn [] (categorical [0.5 0.5]))
x30711 -> (fn [] (categorical [0.5 0.5]))


Observed values O:
--
y30793 -> 3
y30775 -> 1.5
y30739 -> -2.5
y30766 -> -2.2
y30784 -> 2.2
y30811 -> 2.8
y30757 -> -1.9
y30748 -> -1.7
y30730 -> -2.0
y30802 -> 1.2

So there is actually a problem with the way the model is written as the arcs {} are empty.

The following model in Brookes' code, which is very similar, compiles fine in Brookes code, but not with our compiler:

(defn sample-likelihoods [_ likes]
      (let [precision (sample (gamma 1.0 1.0))
            mean (sample (normal 0.0 precision))
            sigma (/ (sqrt precision))]
        (conj likes
              (normal mean sigma))))

    (defn sample-components [_ zs prior]
      (let [z (sample prior)]
        (conj zs z)))

    (defn observe-data [n _ ys zs likes]
      (let [y (nth ys n)
            z (nth zs n)]
        (observe (nth likes z) y)
        nil))

    (let [ys (vector 1.1 2.1 2.0 1.9 0.0 -0.1 -0.05)
          z-prior (discrete
                    (sample (dirichlet (vector 1.0 1.0 1.0))))
          zs (loop 7 (vector) sample-components z-prior)
          likes (loop 3 (vector) sample-likelihoods)]
      (loop 7 nil observe-data ys zs likes)
      zs)

which compiles (in Brookes' compiler) to the following:

; Vertices V: #{x35233 x35236 x35224 x35239 y35285 x35221 x35249 y35299 x35250 x35261 x35227 y35271 y35292 x35242 x35262 x35230 x35255 x35256 y35306 y35278 y35313}
;
; Arcs A: #{[x35255 y35313] [x35255 x35256] [x35250 y35299] [x35255 y35299] [x35261 y35306] [x35249 x35250] [x35250 y35278] [x35255 y35306] [x35262 y35313] [x35242 y35313] [x35221 x35230] [x35227 y35278] [x35249 y35285] [x35261 y35313] [x35256 y35299] [x35261 y35271] [x35261 y35292] [x35256 y35306] [x35249 y35278] [x35255 y35271] [x35250 y35292] [x35262 y35299] [x35221 x35236] [x35236 y35299] [x35256 y35292] [x35255 y35285] [x35230 y35285] [x35250 y35306] [x35250 y35313] [x35261 y35278] [x35250 y35271] [x35249 y35313] [x35221 x35227] [x35262 y35292] [x35262 y35271] [x35221 x35233] [x35233 y35292] [x35249 y35306] [x35261 y35285] [x35224 y35271] [x35262 y35306] [x35239 y35306] [x35262 y35278] [x35256 y35271] [x35221 x35242] [x35256 y35278] [x35250 y35285] [x35255 y35292] [x35261 y35299] [x35249 y35299] [x35256 y35313] [x35221 x35224] [x35249 y35271] [x35221 x35239] [x35256 y35285] [x35249 y35292] [x35261 x35262] [x35255 y35278] [x35262 y35285]}
;
; Conditional densities P:
; x35233 -> (fn [x35221] (discrete x35221))
; x35236 -> (fn [x35221] (discrete x35221))
; x35224 -> (fn [x35221] (discrete x35221))
; x35239 -> (fn [x35221] (discrete x35221))
; y35285 -> (fn [x35230 x35250 x35255 x35261 x35262 x35249 x35256] (nth [(normal x35250 (/ (sqrt x35249))) (normal x35256 (/ (sqrt x35255))) (normal x35262 (/ (sqrt x35261)))] x35230))
; x35221 -> (fn [] (dirichlet [1.0 1.0 1.0]))
; x35249 -> (fn [] (gamma 1.0 1.0))
; y35299 -> (fn [x35250 x35236 x35255 x35261 x35262 x35249 x35256] (nth [(normal x35250 (/ (sqrt x35249))) (normal x35256 (/ (sqrt x35255))) (normal x35262 (/ (sqrt x35261)))] x35236))
; x35250 -> (fn [x35249] (normal 0.0 x35249))
; x35261 -> (fn [] (gamma 1.0 1.0))
; x35227 -> (fn [x35221] (discrete x35221))
; y35271 -> (fn [x35250 x35255 x35224 x35261 x35262 x35249 x35256] (nth [(normal x35250 (/ (sqrt x35249))) (normal x35256 (/ (sqrt x35255))) (normal x35262 (/ (sqrt x35261)))] x35224))
; y35292 -> (fn [x35250 x35255 x35261 x35262 x35233 x35249 x35256] (nth [(normal x35250 (/ (sqrt x35249))) (normal x35256 (/ (sqrt x35255))) (normal x35262 (/ (sqrt x35261)))] x35233))
; x35242 -> (fn [x35221] (discrete x35221))
; x35262 -> (fn [x35261] (normal 0.0 x35261))
; x35230 -> (fn [x35221] (discrete x35221))
; x35255 -> (fn [] (gamma 1.0 1.0))
; x35256 -> (fn [x35255] (normal 0.0 x35255))
; y35306 -> (fn [x35239 x35250 x35255 x35261 x35262 x35249 x35256] (nth [(normal x35250 (/ (sqrt x35249))) (normal x35256 (/ (sqrt x35255))) (normal x35262 (/ (sqrt x35261)))] x35239))
; y35278 -> (fn [x35227 x35250 x35255 x35261 x35262 x35249 x35256] (nth [(normal x35250 (/ (sqrt x35249))) (normal x35256 (/ (sqrt x35255))) (normal x35262 (/ (sqrt x35261)))] x35227))
; y35313 -> (fn [x35242 x35250 x35255 x35261 x35262 x35249 x35256] (nth [(normal x35250 (/ (sqrt x35249))) (normal x35256 (/ (sqrt x35255))) (normal x35262 (/ (sqrt x35261)))] x35242))
;
; Observed values O:
; y35271 -> 1.1
; y35278 -> 2.1
; y35285 -> 2.0
; y35292 -> 1.9
; y35299 -> 0.0
; y35306 -> -0.1
; y35313 -> -0.05

where > == -->

[Tobias-Kohn]

I have significantly enhanced the new PyFOPPL-2 implementation to handle Yuan's example -- or, to be precise, a syntactically simplified version of it. In particular, I added functions such as map, repeat and the for-loop to PyFOPPL-2, so that the model can be written as follows:

(defn sample-components [pi]
  (sample (categorical pi)))

(defn observe-data [y z mus]
  (let [mu (get mus z)]
    (observe (normal mu 2) y)))

(let [ys (vector -2.0  -2.5  -1.7  -1.9  -2.2
                  1.5   2.2   3.0   1.2   2.8)
      pi [0.5 0.5]
      zs  (map sample-components (repeat 10 pi))
      mus (vector (sample (normal 0 2))
                  (sample (normal 0 2)))]
  (for [[y z] (interleave ys zs)] (observe-data y z mus))
  (vector mus zs))

Given the complexity, I think the current PyFOPPL-implementation here will definitely not be able to handle such models.