If you want an accurate answer, then you have to loop over all possibilities of outcomes for Peter’s rolls and compare each one with all the possibilities of outcomes for Colin’s rolls. This way, you will get the exact required probability. The algorithm of using random numbers converges very slowly, that is why you don’t get a correct answer in a reasonable time.
This can be done as follows (you can shrink all loops into just one using Iterator but I’ll keep it simple and verbose here):
using Printf
function dice()
count = 0
total = 0
for p1 = 1:4, p2 = 1:4, p3 = 1:4, p4 = 1:4, p5 = 1:4, p6 = 1:4, p7 = 1:4, p8 = 1:4, p9 = 1:4
for c1 = 1:6, c2 = 1:6, c3 = 1:6, c4 = 1:6, c5 = 1:6, c6 = 1:6
count += p1+p2+p3+p4+p5+p6+p7+p8+p9 > c1+c2+c3+c4+c5+c6
total += 1
end
end
@printf("%.7f\n", count/total) # 0.5731441
end
julia> dice()
0.5731441