1 Introduction

One of the best places on the web where BloodBowl can be played quite comprehensively and which admirably reflects the latest rules published by GW is FUMBBL

And one of the many beauties of this site is that you can spectate the games of other coaches.

Often we see matches on FUMBLL where coaches, true champions of whining, go into a frenzy when they roll a double skull, re-roll, and get again double skull again. ‘Again! But how is this possible! I've been seeing double skulls all week! Ah D6 face die suck! BloodBowl rules suck! This game suck!’.

The problem we pose in this treatise - between the serious and the facetious - is to calculate the probability of having an indented double skull i.e. followed by a subsequent double skull if we roll two block dices (of six faces) N times.
And, as a corollary, try to understand how often this punch in the face can happen to you during your weeks of play.

We will provide a parametric way of calculating this probability P(r,N) as a function of the number N of throws and also in function of the sequence of r consecutive occurrences.

We will give a numerical example for N=70 and then provide a way to graph this probability as N increases.

70 is a number we have defined from experience, which roughly corresponds to the average number of 2D blocks one coach make during a game (we are neglecting the fact that we can make 1D blocks).


2 Methodology

The first step is to reduce the problem to something simple and apply classical definition of Probability: for example, if I flip a coin 3 times, what are my chances of getting two consecutive heads? In probability theory if we use the classical definition of the probability of an event A like

(1)

with NA number of favorable events we have, and N number of total outcomes. So in this case of {"the head twice in a row out of 3 launches"} we have to list all possible outcomes, taking only the favorable ones, dividing them by the total and thus obtaining the probability sought.

In Table 1 we have indicated in red the row where the favorable event (two consecutive heads where head=1) happens. The other row all blacks are events where the two consecutives heads do not happen.


So referring to Table 1 we see that all possible outcome are 8 (all the combination of 3 bit 1,0 from a to h), success happen 3 times (event a, b and e) so based on Eq.1 probability is 3/8=37,5%.

2.1 A possible (wrong) method to calculate Probability of 2 doubleskulls run by using Bernoulli's formula trials

Before looking at the two methods that are correct, we propose a method that is conceptually wrong in solving this problem and which makes use of the Bernouilli Trials formula, quite known in the literature and can be founded for instance in chapter 3-2 of [1].


In that chapter is demonstrated the following fundamental theorem:

THEOREM 1: if we define an event A with a probability p to happen [and complementary probability q=(1-p) not to happen], then the probability that A occurs k times in any order over a number of trials n is

(2)

We can think of simplifying the problem of calculating the event {"at least two double skulls over 70 2d6”} as that of finding probability of the event {“at least one quad skull launching 70 times 4 dices”}.

If we make the calculation of Eq.2 with n=70, k=1, p=(1/36)*(1/36), q=(1-p) we obtain 0.051211 i.e. 5,121%.

Then we say: but in the second case I am launching 70 times 4d6 dices for a total of 280 1d6, where in the first case i launched 70 times 2d6 for a total of 140 1d6. So would be better to calculate “at least one quad skull launching 35 times 4 dices”.

If we make the calculation of Eq. 2 with n=35, k=1, p=(1/36)*(1/36), q=(1-p) we obtain 0.0263071 i.e. 2,63%. Even worst!

if we make the calculation of Eq. 2 with n=69, k=1, p=(1/36)*(1/36), q=(1-p) we obtain 0.05051824 i.e. 5,052%.

Well.

All these results are wrong.

Are wrong in the sense these are results of “different” experiments not the one we defined (double skull followed bu a double skull).

The case of N=35 refers to one experiment where we roll 35 times 4d6 dices and search for probability to get a quad skull one time. Which is NOT the same as saying: I roll 70 times 2d6 and see if there is a double skull followed by a double skull.

Think about it but they are different experiments with respect the one we are analyzing.
Despite the last case which give 5,052% we will see is extremely close as a result to the probability we are looking for of the double-double-skull on 70 throws of 2D6.

We'll see this later.

2.2 Brute force-recursive methods

The brute force method to be used here is to calculate all possible combinations of 2d6 on 70 consecutive throws of the dice. Then calculate the number of situations in which a double skull follows a double skull and do the division. These numbers can be very large, so a calculation programme is necessary.

A very elegant C# code developed by Christer and with the support of Beltroniko calculates the number of combinations and repetitions for a certain number of events (70 events) and a certain die size (36 faces), using a recursive algorithm and memorisation (via a cached dictionary) to speed up calculations of the numerical sequence which gives the result for n=70 in a reasonable low time (see results in Fig. 1)


using System.Diagnostics;
using System.Numerics;
int numberOfEvents = 70;
byte dieSize = 36;
int successes = 1;
Dictionary<int, BigInteger> cache = new();
Console.WriteLine($"Calculating for {numberOfEvents} events and a probability of {successes}/{dieSize}");
Stopwatch stopwatch = new Stopwatch();
stopwatch.Start();
var stats = CalcRecursive(numberOfEvents, dieSize);
stopwatch.Stop();
Console.WriteLine($"Number of combinations: {stats.Count}");
Console.WriteLine($"Number of repeats: {stats.Repeats}");
Console.WriteLine($"Ratio: {(double)stats.Repeats / (double)stats.Count}");
Console.WriteLine($"Time to calculate: {stopwatch.Elapsed}");

Stats CalcRecursive(int number, byte dieSize)
{
if (number == 1)
{
return new(dieSize, 0);
}
var subStats = CalcRecursive(number - 1, dieSize);
var count = subStats.Count * dieSize;
var repeats = subStats.Repeats * dieSize; // Base number of repeats from the sub stats
repeats += Seq(dieSize - 1, number - 1);
return new(count, repeats);
}

BigInteger Seq(int factor, int n)
{
if (cache.ContainsKey(n)) { return cache[n]; }

if (n == 0) { return 0; }
if (n == 1) { return 1; }
var result = factor * (Seq(factor, n - 1) + Seq(factor, n - 2));
cache[n] = result;
return result;
}

record Stats(BigInteger Count, BigInteger Repeats);


Fig.1

The result finally is 5,05562242% rounded to the 8th decimal digit.

2.3 The Probability of a Run

It is not so easy to find in the literature a mathematical formula defining the probability of obtaining a certain event of probability p repeated r times consecutively, or at least not its demonstration. In the article [2] that can be found on the web, a demonstration of a suitable algorithm is given. We refer the reader to the article - which requires a fair amount of mathematical knowledge - by quoting here just the final formula.

let's define

1. 
   
2. n :number of trials
   
3. E :event which single occurrence has probability p
   
4. q :complementary probability of p=1-q.
   
5. r :number of run (consecutive occurrence of the event E
   
6. P_run :Probability of a r “run”, i.e probability that event E happens r consecutive times over n trials
   

In [2] is demonstrated that

(3)
where
(4)
where
(5)

here is the greatest integer contained in n/r+1 or the floor function (the nearest integer less than or equal to that fraction). We can define a MATLAB script to calculate Eq.3

************************************************************************************************

%Matlab script that implement Theorem 1 reported here https://arxiv.org/pdf/math/0511652v1

%Parameters are at the beginning.
% Parameters

p = (1/36); % probability of the event E here: in this case the outcome is "11" (skull skull) in a single 2D launch.
r = 2; % number of consecutive times the event E happen.
q = (1-p); %complementary probability N = 70; %number of trials
N = 70; % number of trials

sum1=Beta(N,r,p,q); %Beta function defined at the bottom
sum2=Beta(N-r,r,p,q); %Beta function defined at the bottom

% Visualization of Beta, Zn and final probability

disp(['p = ', num2str(p), ', q = ', num2str(q), ', N = ', num2str(N), ', r = ', num2str(r)]);
disp('values of BETA as calculated in (4) of https://arxiv.org/pdf/math/0511652v1 formula');
disp(['BETA_(n,r) is : ', num2str(sum1)]);
disp(['BETA_(n-r,r) is : ', num2str(sum2)]);
disp(['p^r is : ', num2str(p^r)]);

finalprobability=sum1-(p^r)*sum2; %formula (5) https://arxiv.org/pdf/math/0511652v1 %final results
disp(['Zn = [Beta_(n,r) -p^r*(Beta_(n-r,r)]: ', num2str(finalprobability)]);

%complementary probability
complementary=1-finalprobability;
disp(['final complementary probability is 1-Zn : ', num2str(complementary)]);

%definition of function BETA(p1,p2,p3,p4) as (4) in https://arxiv.org/pdf/math/0511652v1
% BETA FUNCTION % Note: do not use factorial in MATLAB because besides N= 170 return NaN
%nchoosek(a,b) can be used to calculate the binomial, but when number start becoming high
%start giving several warning and limited numerical precision.
% For this we used gammaln for the factorial calculation .
%gammaln calculate the ln( of gammaN) gammaN is (n-1)! so gammaln is equivalent to factorial log.

function result = Beta(p1, p2, p3, p4)
sum_1 = 0; % sum initialization to 0
a1 = -1; % this is the parameter of forumle (4) -1
% Beta n,r summatory calculation : formula (4) in https://arxiv.org/pdf/math/0511652v1

for l = 0:(floor(p1/(p2+1)))
% Verification that values are valid
if (p1 - l*p2) >= l && (p1 - l*p2 - l) >= 0
% Binomial calculation by using gammaln
binomial_coefficient = exp(gammaln(p1 - l*p2 + 1) - gammaln(l + 1) - gammaln(p1 - l*p2 - l + 1));
% Sum of current term to the summatory
sum_1 = sum_1 + ((a1^l) * binomial_coefficient * ((p4 * (p3^p2))^l));
else
% Debug: shows when (if) values are not valid
disp(['UNVALID values for binominal coefficient : N = ', num2str(p1), ', l = ', num2str(l)]);
end
end


% Managing exceptions or particular cases in case NaN (this is a debug can be removed)

if isnan(sum_1)
disp(['Somma finale è NaN per p1 = ', num2str(p1), ', p2 = ', num2str(p2)]);
end

result = sum_1;
end

************************************************************************************************

Beta function is defined- as usual - at the END of MATLAB script.

To avoid problems with very large numbers, we decided to use logarithms. MATLAB does not have a built-in function for directly calculating very large binomial coefficients with high precision, but it is possible to calculate the logarithms of the binomial coefficients and then perform exponentiation.

For this purpose gammaln function can be used to calculate the logarithms of the factorials.

This approach is more numerically stable than calculating the factorials and binomial coefficients directly.

At the end the result is the same of the brute force algorithm


Fig.2

We can add a piece of code to MATLAB to plot the value of Eq.3 as N increase.

We quote the MATLAB code below which should be added just before the definition of the function Beta


%***************************************************

% plot the probability diagram over N

figure;
plot(N_values, complementary_values, 'b-', 'LineWidth', 2);
xlabel('N');
ylabel('Total_probability');
title('Probability evolution in function of N: number of 2D rolls ');
xlim([50 1000]); % Range of X axis from 50 to 1000
ylim([0 1]);
grid on;
% ***************************************************


Which give the following results in Fig.3



Figure 3

Figure 3 shows that we reach 50% of probability to see a double skull followed by a double skull after 900 2D6rolls.

Considering in average 60-70 2D6 rolls per each match by each players (for a total of 120-140 2D6 rolls) this means that the "snake" starts to become likely after 9-10 matches.

Finally we note that the P_run calculated for n=70 with this algorithm is not so far away from the original “rough” assumption of the quad skull by using Eq.2 with n=69, r=1, p=(1/36)*(1/36) which resulted in 5.052%, which is not that far from 5.055% of Fig.1. To the brave and enterprising we leave this dilemma as an exercise in the end

3 Conclusions

In conclusion, as always, the more times we roll the dice, the faster we will see a snake arrive. But if you don't buy a lottery ticket you are only sure not to win, but if you want to win sometimes you have to spend that risk money. So: Buy a ticket.!

It's all about knowing when it's time to take a risk.

A curiosity. If you put r=0 in the formula, the result is 100%. That is, it is certain that any event E of probability P repeated any number of times has 0 consecutive repetitions. Think about it, it seems obvious... But is it? Again:More exercises to pass the sleepless night!

4 Exercises

EXERCISE 1
Explain why if we apply the Eq.2 with n=69, k=1, p= (1/36)*(1/3) the result = 0,052% (calculated in subsection 2.1) is very similar (though not exactly the same) to the one calculated by the algorithms in subsection 2.2

exercise 1.1
Is there a way with the Eq.(2) to obtain the same probability calculated in subsection 2.2?

EXERCISE 2
Explain using only logic without mathematical formulae why the probability of a succession of 0 consecutive double skulls has a 100% chance of occurring .


EXERCISE 3
Calculate the probability of obtaining a triple skull followed by a triple skull during a game in which 146 block dice are rolled and explain the result. If the problem seems misplaced, explain why.


Acknowlegment

This paper is a summary of a nice discussion we had in Fumbbl's Discord group with Christer and Beltroniko whom I thank for their patience with my poorly reasoned objections.
The C# code was developed by Christer, who also pointed out my numerous logical errors.

--------------------------------
References

[1] Athanasios Papolis, Probability, Random Variables, and Stochastic Processes International edition, 1991

[2] Villarino Mark B., The Probability of a Run, 2008, https://arxiv.org/pdf/math/0511652v1