I made this short code to calculate the chances of a success rolling dice, and it worked very well... but not in big numbers. Se the code, I'll explain better below.
def calc_dados(f_sucessos = 1, faces = 6, n_dados = 1):
p_max = ((f_sucessos/faces)**n_dados) #chance de todos
fator = 1
p_meio = 0
for i in range(n_dados-1):
p_meio += (((f_sucessos/faces)**(n_dados-fator) * ((faces-f_sucessos)/faces)**(n_dados-(n_dados-fator))) * n_dados)
fator += 1
p = p_max + p_meio
return p*100
So, ok, it works, why not see how my chances are better in function of adding dice? More the dice, better the chance. So I made this tiny table with pandas:
f_sucessos = 1 # how many faces are success
faces = 2 # faces of the dice
n_dados = 10 # n de dados lançados
suc_list = []
for i in range(0,n_dados): suc_list.append(f_sucessos)
fac_list = []
for i in range(0,n_dados): fac_list.append(faces)
cha_list = []
for i in range(0,n_dados): cha_list.append(calc_dados(f_sucessos, faces, i+1))
df = pd.DataFrame(
{
"n_dados" : range(1,n_dados+1),
"faces" : fac_list,
"sucessos" : suc_list,
"chance" : cha_list
}
)
df
The results were very strange... So I wrote an coin probability table and tested as the coin was an 2 faced dice. The right table is this: table of right brute force tested results
But if you use my code to create this table the result will be this: table of the results of my code
Please, anybody can help me to understood why in a certain moment the probabilities just fall when they should be higher? For example:The chance of at least 1 'head' in 4 coins should be 93,75%, but my code says it is 81,25%...
Hello people, it's finally solved!
First I would like to thank Raibek, who solved it using combinations, I didn't realize it was solved when he did it and below I'll tell you how and why.
If you are not following the history of this code, you just need to know that it is used to calculate the probability of getting at least ns successes when rolling d amount of dice. Solution codes are at the end of this answer.
I found out how to solve the problem by talking to a friend, Eber, who pointed me to an alternative to check the data, anydice.com. I quickly realized that my visual check, assembling tables in Excel/Calc was wrong, but why?
Well, here comes my friend who, reading the table of large numbers with 7d6, where the error was already very evident, shows me that although at the beginning the account worked, my table did not have all the possible combinations. And the more possibilities there were, the more my accounts failed, with the odds getting smaller as more dice were added to the roll.
This is the combinations I was considering, in this example on 7d6 case.
In the first code the account was:
successes**factor *failures**factor *d
The mistake is in assuming that the number of possible combinations was equal to d (which is a coincidence up to 3 dice for the tests I did before thanks to factorials of 1 = 1 and factorial of 2 = 2).
Now notice that, in 7d6 example, in the exact 3 block there are some missing possible combinations in yellow:
The correct account for this term of the equation is:
factorial(d) / factorial (failures) * factorial (successes)
With this account we can find out what the chance of exactly n faces rolling is, and then if we want, for example, to know the chance of at least once getting the number 1 in 3d6, we just need to add the chances of getting exactly 1 time, 2 times and 3 times. What the code already did well.
Finally, let's get to the code:
Daniel-Eber solution:
def dado(fs=1,ft=6,d=1,ns=1,exato=False):
'''
fs = faces sucesso
ft = faces totais
d = n de dados
ns - n de sucessos esperados modificados por exato
exato = True: chance de exatamente ns ocorrerem, False: chance de pelo menos ns ocorrerem
'''
from math import factorial
s = fs/ft
f = (ft-fs)/ft
d = d
ns = ns
p_max = s**d
falhas = 1
po_soma = 0
if exato == False:
for i in range(d-1):
po = ( (s**(d-falhas)) * (f**(falhas))) * (factorial(d)/(factorial(falhas)*factorial((d-falhas))))
po_soma += po
falhas += 1
else:
p_max = 0
falhas = d-ns
po_soma = ( (s**(d-falhas)) * (f**(falhas))) * (factorial(d)/(factorial(falhas)*factorial((d-falhas))))
return f'{(p_max + po_soma)*100:.2f}%'
print(dado(1,6,6,1))
Raibek solution:
from scipy.special import comb
def calc_dados(f_sucessos=1, faces=6, n_dados=1):
assert f_sucessos <= faces
outcomes_total = faces ** n_dados
outcomes_success = 0
f_fail = faces - f_sucessos
for i in range(1, n_dados + 1):
one_permutation = (f_sucessos ** i) * (f_fail ** (n_dados - i))
n_permutations = comb(n_dados, i)
outcomes_success += one_permutation * n_permutations
p = outcomes_success / outcomes_total
return f'{(p)*100:.2f}%'