# Composing probabilistic models

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
"""
for i=1:length(v)
v[i] ~ Normal(speed, 5.0)
end
end

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

@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

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

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)
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

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 . I will keep an eye on it.

1 Like

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

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)
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`.

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

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