模拟一个概率问题：3个独立的骰子

Question

我决定模拟一道教科书上的概率题：

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)

假设骰子是公平的，因此“理论”答案将是 $\frac{\binom{5}{2}}{\binom{6}{3}}=0.5$；我决定在 Julia 中对此进行模拟，这是我为此编写的函数：

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

我预计这是 return 大约 0.5，但实际上在大约 10^5 次迭代后 returns ~0.27 (n_experiments)。谁能帮我找出代码中的缺陷？

Answer 1

我认为你的代码是正确的，但你的数学错误：

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

我希望看到一种简洁的方式来对条件进行编码。由于分配，我避免了你的变体，而我的变体只是冗长确保它是正确的——但太长了……

附录

这是具有生成条件的优化变体。与问题无关，但因为我觉得它很酷。条件公式归功于 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

这根本没有分配。 sum 与手动循环一样快！

模拟一个概率问题：3个独立的骰子

Simulating a probability problem: 3 independent dice

julia

probability

montecarlo

附录