Monte-Carlo Simulation for the sum of die

Question

I am very new to programming so I apologise in advance for my lack of knowledge.

I want to find the probability of obtaining the sum k when throwing m die. I am not looking for a direct answer, I just want to ask if I am on the right track and what I can improve.

I begin with a function that calculates the sum of an array of m die:

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

Now I am trying specific values to see if I can find a pattern (but without much luck). I have tried m = 2 (two die). What I am trying to do is to write a function which checks whether the sum of the two die is k and if it is, it calculates the probability. My attempt is very naive but I am hoping someone can point me in the right direction:

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

I want to somehow find the number of 'elements' in dicesum(2) in order to calculate the probability. For example, consider the case when dicesum(2) = 8. With two die, the possible outcomes are (2,6),(6,2), (5,3), (3,5), (4,4), (4,4). The probability being (2/36)*3.

I understand that the general case is far more complicated but I just want an idea of how to being this problem. Thanks in advance for any help.

Answer 1

If I understand correctly, you want to use simulation to approximate the probability of obtaining a sum of k when roll m dice. What I recommend is creating a function that will take k and m as arguments and repeat the simulation a large number of times. The following might help you get started:

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)

The estimate is approximately .0131.