# Buckets of dice - Binomial probabilities of either event occuring

I'm doing a bit of analysis on a dice rolling mechanic for a role playing game. This is a buckets of dice system where the result is based on the count of the number of dice rolling a given target value or below.

There is a twist - rolling a 1 one any of the dice is a critical success, which contributes 2 to the count.

The probabilities I'm trying to calculate are the probability of a given number of dice rolling two or more successes at a given probability, or at least one die rolling a 1. By Bayes's Theorem, P(AUB) = P(A) + P(B) - P(A∩B). In this case I've used scipy to calculate the probabilities of the individual events.

I'm trying to calculate the probability of two or more dice from a given pool of N (Say: D10s) rolling the target value or below, or one or more dice from the pool rolling a 1.

The script below uses scipy.stats to calculate the probabilities for the individual events. I'm trying to calculate the total probability of either occurring to see the overall probability of success with the dice roll.

One can trivially calculate P(A) and P(B) using a straight binomial survival function from the library. How does one go about computing P(A∩B) for these events?

# === BucketsOfDice.py ======================================
#

MIN = 8
MAX_DICE = 12
DIE = 10

import numpy
from scipy.stats import binom

# First compute the value of two or more dice rolling successes
# at the target value for N dice and a range of modifiers.
# The modifiers add or take one die from the roll. This is
# computed using the survival function of a binomial distribution
# parameterised with the number of dice and the probability of
# the target roll
#
for min in [2, 3, 4, 5, 6, 7, 8, 9, 10]:
prob = (DIE - min + 1) / float (DIE)
print
print 'Base: %d- (%0.1f),3,4,5,6,7,8,9,10,11,12' % (DIE - min + 1, prob)
for DM in [+2, +1, 0, -1, -2, -3, -4, -5, -6]:
csv=[]
csv.append('DM: %d' % (DM))
for skill in [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]:
if (skill + DM) < 1:
_sk = 1     # At least one die is always rolled.
else:
_sk = skill + DM
dist = binom (_sk, prob)
hit = dist.sf(1)
csv.append(',%0.2f' % (hit))
output = ''
for val in csv:
output = output + val
print output

# Now do the same computation for at least one die rolling a 1
# for a critical that counts as two successes (10% for a D10).
#
print
print 'Critical,3,4,5,6,7,8,9,10,11,12'

for DM in [+2, +1, 0, -1, -2, -3, -4, -5, -6]:
csv=[]
csv.append('DM: %d' % (DM))
for skill in [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]:
if (skill + DM) < 1:
_sk = 1
else:
_sk = skill + DM
dist = binom (_sk, 0.1)
hit = dist.sf(0)
csv.append(',%0.2f' % (hit))
output = ''
for val in csv:
output = output + val
print output


Adding a clarification for the comment. There are two things that can happen on the dice that are of interest. For some given number of dice:

• Two or more of the dice can roll under the target value (e.g. 3 or less).

• One or more of the dice can roll a 1, irrespective of the target value.

I can calculate the probability of either event using a binomial distribution probability calculator, but the events are not discrete (i.e. both can occur at the same time) so I am trying to calculate the probability of either event occurring.

• Are you asking for Python code, or how such probabilities are calculated mathematically? – gung - Reinstate Monica Jan 16 '16 at 13:54
• How to calculate the probabilities - I just used scipy to calculate PA and PB as it has a binomial calculator in the library. I'm not married to it, though; I'd be quite happy to implement the computation some other way. – ConcernedOfTunbridgeWells Jan 16 '16 at 14:12
• Based on the text--I didn't read the code--I haven't been able to figure out what probabilities you are trying to compute. Could you be a little clearer about that or perhaps offer an example? – whuber Jan 16 '16 at 15:00
• You could modify the the 1st part of the algorithm to count only 2s 3s up to the target value (i.e. specifically exclude 1s there). So, for a target of 4, the odds of making that with d10 would be treated as 0.3 instead of 0.4. I'm not familiar with scipy, so I can't really give you a code recommendation - sorry! – MikeP Jan 19 '16 at 14:21