模具总和的蒙特卡洛模拟

Question

我是编程新手，所以我为我的知识不足提前道歉。

我想求掷m个骰子得到和k的概率。我不是在寻找直接的答案，我只是想问问我是否在正确的轨道上，我可以改进什么。

我从计算 m 个骰子数组总和的函数开始：

function dicesum(m)
j = rand((1:6), m)
sum(j)
end

现在我正在尝试特定值以查看是否可以找到模式（但不太走运）。我试过 m = 2 （两个死）。我想做的是编写一个函数来检查两个骰子的总和是否为 k，如果是，则计算概率。我的尝试很天真，但我希望有人能指出我正确的方向：

m = 2
x, y = rand(1:6), rand(1:6)
z = x+y
if z == dicesum(m)
Probability = ??/6^m

我想以某种方式在 dicesum(2) 中找到 'elements' 的数量，以便计算概率。例如，考虑 dicesum(2) = 8 的情况。有两个骰子，可能的结果是 (2,6),(6,2), (5,3), (3,5), (4,4 ), (4,4)。概率为 (2/36)*3.

我知道一般情况要复杂得多，但我只想知道如何解决这个问题。在此先感谢您的帮助。

Answer 1

如果我没理解错的话，你想用模拟来估计掷m个骰子得到总和为k的概率。我建议创建一个将 k 和 m 作为参数并多次重复模拟的函数。以下内容可能会帮助您入门：

function Simulate(m,k,Nsim=10^4)
    #Initialize the counter
    cnt=0 
    #Repeat the experiment Nsim times
    for sim in 1:Nsim
        #Simulate roll of m dice
        s = sum(rand(1:6,m))
        #Increment counter if sum matches k
        if s == k 
            cnt += 1
        end
    end
    #Return the estimated probability
    return cnt/Nsim
end

prob = Simulate(3,4)

估计值约为 .0131。

Answer 2

您还可以如下所示以矢量化方式执行模拟。它在内存分配方面效率较低，因为它创建了一个长度为 Nsim 的向量 s，而循环代码使用单个整数 cnt 进行计数。有时不必要的内存分配会导致性能问题。在这种情况下，事实证明矢量化代码的速度大约是原来的两倍。通常，循环会更快一些。更熟悉 Julia 内部结构的人可能会提供解释。

function Simulate1(m,k,Nsim=10^4)
    #Simulate roll of m dice Nsim times
    s = sum(rand(1:6,Nsim,m),2)
    #Relative frequency of matches
    prob = mean(s .== k)
    return prob
end

模具总和的蒙特卡洛模拟

Monte-Carlo Simulation for the sum of die

julia

montecarlo