Turing: Observing sum of two random dice

Specific Domains Probabilistic Programming
1

I’m (still) a Turing beginner and tries to build some simple models (porting my WebPPL models). I now wonder how observing a sum of two (DiscreteUniform) random variables works, i.e. to force the sum to a specific value.

Here’s a simple model which sums the value of two dice:

@model two_dice(ss) = begin
    d1 ~ DiscreteUniform(1, 6) # die 1
    d2 ~ DiscreteUniform(1, 6) # die 2
    s ~ DiscreteUniform(2, 12) # sum of d1+d2

    s = d1 + d2

    # Force s to be the parameter ss 
    s == ss || begin Turing.@addlogprob! -Inf; return end
end

# Now we observe that the sum is 2
model = two_dice(2)

num_chains = 1 
chains = sample(model, MH(), MCMCThreads(), 10000, num_chains)
display(chains)

This model works, kind of: The display of chains is:

Chains MCMC chain (10000×4×1 Array{Float64,3}):

Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = d1, d2, s
internals         = lp

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat 
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64 

          d1    1.0032    0.0799     0.0008    0.0032   574.4010    1.0015
          d2    1.0064    0.1599     0.0016    0.0064   574.4010    1.0015
           s    6.9814    3.1706     0.0317    0.2331   145.6114    1.0078

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

          d1    1.0000    1.0000    1.0000    1.0000    1.0000
          d2    1.0000    1.0000    1.0000    1.0000    1.0000
           s    2.0000    4.0000    7.0000   10.0000   12.0000

I.e. both d1 and d2 are (correctly) forced to be 1, but s is all over the place (values in 2…12). I had assumed that s should be forced to be 2.

Note that this only works with MH(); what I can see none of the other samplers forces d1 and d2 to 1.

Here are my questions:

  • Is s = d1 + d2 the correct way to making the “deterministic” sum of the two dice d1 and d2? Or are there better way to do this?
  • Is the code
    s == ss || begin Turing.@addlogprob! -Inf; return end
    the best (or even good) way of observing (forcing) the sum s?
  • Does it make sense that only MH() can handle these kind of problem: neither IS, PG or SMC forces d1 and d2.

A variant of the above model is the following where s is observed via the model parameter sdirectly (i.e. without any use of @addlogprob! via the paramter ss):

@model two_dice(s) = begin
    d1 ~ DiscreteUniform(1, 6)
    d2 ~ DiscreteUniform(1, 6)
    s ~ DiscreteUniform(2, 12)

    # Observe the parameter directly 
    s = d1 + d2
end

model = two_dice(2)
# ...

This was actually the first approach I tried, but this don’t seems to connect s with the input value (here 2) at all: both d1 and d2 are all even distributed in 1…6 (s is not available since it’s an observed variable).

Related: I have probably been a little spoiled by some other PPLs I’ve tested (such as WebPPL, BLOG, and PSI) which all have some kind of explicit observe function that seems to do the same thing as the code in the first model i.e. adding -Inf to the logprob. I tried to make an @observe macro (or observe function) for this but didn’t succeed. It would be nice if such a macro could be supported.

I might very well have missed some obvious thing here, but I could no find discussions about this in the excellent Turing documentation.

2

I’m (still) a Turing beginner and tries to build some simple models (porting my WebPPL models).

Welcome!

s ~ DiscreteUniform(2, 12)

Note that you are sampling a value for s using this line which is the value reported in the chain but then you define s to be d1 + d2 later anyways. So the sampled value is thrown away. The way you are doing sampling is correct but it won’t work with samplers that require differentiating the log probability because it’s either 0 or -Inf in your case. If other samplers that don’t require derivatives don’t work, please open an issue about that.

In Turing we don’t yet have a way to store intermediate results like your s = d1 + d2. There were some recent discussions and there are plans but it hasn’t been done yet. So for now, you need to recompute any intermediate values from the chain’s sample.

1 Like
3

@mohamed82008 Thanks for your answer.

It make sense now that you explain it, i.e. that s = d1 + d2 replace the value of s which was generated in s ~ DiscreteUniform(2,12).

So, until this is (perhaps) fixed, I will use this constructs which is OK for me:

@model two_dice(ss) = begin
   d1 ~ DiscreteUniform(1, 6)
   d2 ~ DiscreteUniform(1, 6)

   d1+d2 == ss || begin Turing.@addlogprob! -Inf; return end
end

model = two_dice(4)
# ...

Regarding that only MH() that support this now: what I understand the two particle based samplers IS and PG should also work. Is that correct?

Also, do you have an idea how to write the macro (@observe) to do the Turing.@addlogprob! in a neater way?

4

We don’t have an @observe macro. But perhaps we should. Please open an issue to discuss this further.

1 Like
5

@mohamed82008 Thanks. I’ll open an issue (but first I’ll try to implement it myself. :-))

1 Like