Composing probabilistic models

Specific Domains Probabilistic programming
1

Hi, I’m quite new to Julia and I have stumbled upon Turing. It looks great and the documentation is nice but either I have missed something or simply it’s not possible.

I want to define models in a modular way and combine them. For example imagine I want to estimate the speed of a car moving at constant speed and I have two sensors: gps for measuring positions and radar for directly measuring speed. The following simple script shows that only copy-pasting in a single function works

using Distributions
using Turing

"""
x: position measure 
t: time measure 
speed: parameter
"""
@model function gps(t, x, speed)
    for i=1:length(x)
        x[i] ~ Normal(speed * t[i], 1.0)
    end
end

"""
v: speed measure
speed: parameter
"""
@model function radar(v, speed)
    for i=1:length(v)
        v[i] ~ Normal(speed, 5.0)
    end
end

"Compose gps + radar"
@model function fullmodel1(t, x, v)
    speed ~ Normal(0.0, 100.0)
    gps(t, x, speed)
    radar(v, speed)
end

"Copy-paste GPS + radar"
@model function fullmodel2(t, x, v)
    speed ~ Normal(0.0, 100.0)
    for i=1:length(x)
        x[i] ~ Normal(speed * t[i], 1.0)
    end
    for i=1:length(v)
        v[i] ~ Normal(speed, 5.0)
    end
end

"Compare fullmodel1 vs fullmodel2"
function testmodel(N=100)    
    speed = 5.0  # true value of speed
    # generate measurements
    t = LinRange(0, 120, N)
    x = speed .* t .+ randn(N)
    v = speed .+ randn(N)
    # estimate speed
    c1 = sample(fullmodel1(t, x, v), NUTS(0.6), 1000)
    c2 = sample(fullmodel2(t, x, v), NUTS(0.6), 1000)
    c1, c2
end

c1, c2 = testmodel()
println("Results for fullmodel1")
display(c1)
println("Results for fullmodel2")
display(c2)

The script will run but only fullmodel2 gives correct results. Eliminating speed as a parameter from gps or radar gives identical results.

Since I don’t know what magic is @model really doing behind the scenes I don’t know how to proceed and any help is welcomed.

I know that for such a simple example splitting into several functions is overkill but for more complex models I think is extremely important. For example, thinking again in terms of localization problems, I may have different models for accelerometers, gyroscopes, GPS just depending on the model or type of sensor (laser or MEMS gyros, differential GPS, relative GPS…) although all of them give back the same types of measures and so the structure of the problem doesn’t change, just the internal details of each part.

Thanks!

2 Likes
2

Soss.jl targets exactly this kind of use. You could write your example in Soss as

gps = @model speed,n,t begin
    x ~ For(1:n) do i
        Normal(μ = speed * t[i])
    end
    return x
end

radar = @model speed,n begin
    v ~ For(1:n) do i
        Normal(speed, 5.0)
    end
    return v
end

fullmodel = @model n,t begin
    speed ~ Normal(0.0, 100.0)
    x ~ gps(;speed, n, t)
    v ~ radar(;speed, n)
end

You can sample it like this:

n = 10
t = 1:10

rand(fullmodel(;n, t))

We’re in the middle of some big changes:

This is all in the dev branch. It looks like I have a couple of things to get back into place before I can sample from the posterior for this example. But it’s a great example, and I think we’ll add it as a unit test.

6 Likes
3

Thank you for your answer. Yes, this is exactly what I’m looking for. I have tested the dev branch and I can sample the model although as you say it’s not possible yet to sample from the posterior, so for the time being is not a solution for me :slight_smile: . I will keep an eye on it.

1 Like
4

We’ve had this kind of thing working before, just need to reconnect some wires. Shouldn’t be too long.

This kind of compositional programming model is not too common in PPL, so I’d like to have some examples to show the value. Please let me know if you have other ideas for nice examples like this.

1 Like
5

I had a silly mistake. There’s something that comes up a lot in most PPLs that you can only observe values that are sampled. But in this line:

radar = @model speed,n begin
    v ~ For(1:n) do i
        Normal(speed, 5.0)
    end
    return v
end

that return v could have computed any function of v. In this case it’s clear what we meant to do, but it takes some work to set things up in a way that generalizes well.

For now, something like this works:

gps(n,speed,t) = For(1:n) do i
    Normal(μ = speed * t[i])
end

radar(n, speed) = For(1:n) do i
    Normal(speed, 5.0)
end

fullmodel2 = @model n,t begin
    speed ~ Normal(0.0, 100.0)
    x ~ gps(n,speed, t)
    v ~ radar(n, speed)
end

n = 10
t = 1:10

truth = rand(fullmodel2(;n, t))
v = truth.v

post = fullmodel2(;n, t) | (;v)

dynamicHMC(post)

You could also pass gps and radar as parameters, or have the For constructors in the model so gps and radar don’t need to worry about n.

6

I think a killer example would be composing bolting on a GLM onto a domain specific model (e.g. of behaviour).

7

Sounds interesting, but I don’t really understand yet. Could you give some more details?