'Simulating a probability problem: 3 independent dice

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?



Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source