Turing constrain on observations of deterministic function

Specific Domains Probabilistic Programming
1

I have the following model

function f(x)
    x^2
end


@model function simplemodel()
    x ~ Normal(5, 1)
    y := f(x)
    return y
end

sample(simplemodel(), NUTS(), 50) |> histogram

As can be seen it involves the deterministic function f. Now lets say that I have some observation for y = f(x), such as y=3.

What I am trying to do is to constrain the model, like this simplemodel() | (y = 3.0,).

However the results are not as I would have expected. In fact this observation does not seem to change anything on the model. I would have expected similar results as when setting simplemodel() | (x = sqrt(3.0),), after which there is only one value for y. Instead nothing happens.

Follow up question: What if my observation of y is not a single value, but a value with an uncertainty, like y = 3 \pm 0.1 = \mathcal N(3,0.1)?

2

I don’t know the details of the DynamicPPL/Turing implementation, but it would appear that the probabilistic model can’t reason backwards from deterministically derived quantities. (Which is not generally something you want to do anyway.)

If you replace y := f(x) with y ~ Normal(f(x), 0.1), it will behave how you want it to, as it’s now a standard Bayesian problem with a prior and likelihood. Making the observation error very small (e.g. 1e-6) will give you a good approximation of your first, deterministic, situation.

1 Like
3

I also tried that, I think its a good idea. For that application there is also the Dirac Distribution, if I want to say that y is precicely known.

However the sampling does not seem to work at all for small standard deviations. It seems to only obtain a single inital sample and then the markov chain does not move any way from it:

using Turing, Plots, StatsPlots

function f(x)
    global nevals+=1
    x^2
end


@model function simplemodel()
    x ~ Normal(5, 1)
    y ~ Normal(f(x), 0.001)
    # y ~ Dirac(f(x))
    # return y
end

nevals=0
chain = sample(simplemodel(), MH(), 1000)
histogram(chain)
plot(chain)
@show nevals
4

Yeah, that’s not surprising to me. I wouldn’t expect the Dirac distribution to work for sampling at all, and when the standard deviations of the two parameters are very different, the sampler (NUTS, and especially MH) may struggle to get good proposal step lengths for both of them.

If you can share a bit more information about your actual problem/objective, it may be possible to give better advice on how to approach it…