Simulating a probability problem: 3 independent dice

Question

I decided to simulate a probability question from a textbook:

Three fair dice are rolled independently, what is the probability that one dice shows 6 and the other two show two non-equal numbers (and neither is equal 6)

The dice is assumed fair, so "theoretical" answer will be $\frac{\binom{5}{2}}{\binom{6}{3}}=0.5$; I decided to simulate this in Julia, here is the function I wrote for this:

function simulation(n_experiments::Int = 100)
    dice = [DiscreteUniform(1, 6) for i in 1:3] # 3 independent fair dice
    successes = 0
    for trial in 1:n_experiments
        experiment = rand.(dice) # roll
        one_six = sum(r == 6 for r in experiment) == 1 # check if there is only one "6"
        others = filter(num -> num != 6, experiment)
        no_duplicates = length(Set(others)) == 2 # check if other two are distinct
        if no_duplicates
            if one_six
                successes += 1 # count "success"
            end
        end
    end
    return successes / n_experiments
end

I expected this to return something around 0.5, but in fact it returns ~0.27 after about 10^5 iterations (n_experiments). Can anyone please help me find the flaw in the code?

Answer 1

I think your code is right, but your math wrong:

julia> function simulation(n_experiments=100)
           chains = rand.(fill(DiscreteUniform(1, 6), 3), n_experiments)
           return (1/n_experiments) * mapreduce(+, chains...) do d1, d2, d3
               (d1 == 6 && d2 != 6 && d3 != 6 && d2 != d3) || (d1 != 6 && d2 == 6 && d3 != 6 && d1 != d3) || (d1 != 6 && d2 != 6 && d3 == 6 && d1 != d2)
           end
       end
simulation (generic function with 1 method)

julia> simulation(10_000)
0.279

julia> count(Iterators.product(1:6, 1:6, 1:6)) do (d1, d2, d3)
           (d1 == 6 && d2 != 6 && d3 != 6 && d2 != d3) || (d1 != 6 && d2 == 6 && d3 != 6 && d1 != d3) || (d1 != 6 && d2 != 6 && d3 == 6 && d1 != d2)
       end 
60

julia> 60 / 6^3
0.2777777777777778

I'd love to see a concise way to encode the condition. I avoided your variant due to allocations, and mine is only verbose make sure it's correct -- but way too long...

---

Addendum

Here's the optimized variant with a generated condition. Not relevant for the question, but because I find it cool. Credit for the condition formula goes to Bogumił.

julia> @generated function condition(xs::NTuple{N,<:Any}) where {N}
           offdiagonals = ((i, j) for i = 1:N for j = (i+1):N)
           names = Symbol.(:xs, 1:N)
           assignments = [:($(names[i]) = xs[$i]) for i = 1:N]
           comparisons = [:($(names[i]) != $(names[j])) for (i, j) in offdiagonals]
           condition = Expr(:&&, :(maximum(xs) == 6), comparisons...)
           return quote
               $(assignments...)
               $condition
           end
       end
condition (generic function with 1 method)

julia> @code_warntype condition((1,2,3))
MethodInstance for condition(::Tuple{Int64, Int64, Int64})
  from condition(xs::Tuple{Vararg{var"#s17", N}} where var"#s17") where N in Main at REPL[71]:1
Static Parameters
  N = 3
Arguments
  #self#::Core.Const(condition)
  xs::Tuple{Int64, Int64, Int64}
Locals
  xs3::Int64
  xs2::Int64
  xs1::Int64
Body::Bool
1 ─       (xs1 = Base.getindex(xs, 1))
│         (xs2 = Base.getindex(xs, 2))
│         (xs3 = Base.getindex(xs, 3))
│   %4  = Main.maximum(xs)::Int64
│   %5  = (%4 == 6)::Bool
└──       goto #7 if not %5
2 ─ %7  = (xs1 != xs2)::Bool
└──       goto #6 if not %7
3 ─ %9  = (xs1 != xs3)::Bool
└──       goto #5 if not %9
4 ─ %11 = (xs2 != xs3)::Bool
└──       return %11
5 ─       return false
6 ─       return false
7 ─       return false


julia> function simulation(n_experiments=100)
           (1/n_experiments) * sum(1:n_experiments) do _
               dice = (rand(1:6), rand(1:6), rand(1:6))
               condition(dice)
           end
       end
simulation (generic function with 2 methods)

julia> @time simulation(10_000_000)
  0.157303 seconds
0.2777498

This has no allocations at all. sum is exactly as fast as a manual loop!