# Turing: Observing sum of two random dice

**URL:** <https://discourse.julialang.org/t/turing-observing-sum-of-two-random-dice/50486>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [November 20, 2020, 10:54am UTC](https://discourse.julialang.org/t/turing-observing-sum-of-two-random-dice/50486 "2020-11-20T10:54:31Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Hakan\_Kjellerstrand](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hakan_kjellerstrand/32/8447_2.png) [@Hakan\_Kjellerstrand](https://discourse.julialang.org/u/Hakan_Kjellerstrand)\
**Post date:** [November 20, 2020, 10:54am UTC](https://discourse.julialang.org/t/turing-observing-sum-of-two-random-dice/50486/1 "2020-11-20T10:54:31Z")

</div>

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:

```julia
@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:

```julia
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 `s`directly (i.e. without any use of @addlogprob! via the paramter `ss`):

```julia
@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.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [November 20, 2020, 12:50pm UTC](https://discourse.julialang.org/t/turing-observing-sum-of-two-random-dice/50486/2 "2020-11-20T12:50:07Z")

</div>

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

Welcome!

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

---

<div class="post-metadata">

**Author:** ![Hakan\_Kjellerstrand](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hakan_kjellerstrand/32/8447_2.png) [@Hakan\_Kjellerstrand](https://discourse.julialang.org/u/Hakan_Kjellerstrand)\
**Post date:** [November 20, 2020, 4:26pm UTC](https://discourse.julialang.org/t/turing-observing-sum-of-two-random-dice/50486/3 "2020-11-20T16:26:13Z")

</div>

@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:

```julia
@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?

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [November 21, 2020, 2:53am UTC](https://discourse.julialang.org/t/turing-observing-sum-of-two-random-dice/50486/4 "2020-11-21T02:53:33Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Hakan\_Kjellerstrand](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hakan_kjellerstrand/32/8447_2.png) [@Hakan\_Kjellerstrand](https://discourse.julialang.org/u/Hakan_Kjellerstrand)\
**Post date:** [November 21, 2020, 8:29am UTC](https://discourse.julialang.org/t/turing-observing-sum-of-two-random-dice/50486/5 "2020-11-21T08:29:03Z")

</div>

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