Catan Is Broken. Unplayable

[Gulchen]

I finally ran an analysis on the game yfimc posted percentages from: ​ https://boardgamearena.com/table?table=445859988

The histogram of 55 rolls from that game has a p-value which is strictly between ​ ​ ​ 1 / 552 ​ ​ ​ and ​ ​ ​ 1 / 551 ​ .

Code: Select all

Python 3.8.0 (default, 
[GCC 8.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from math import floor,factorial,ceil
>>> from itertools import combinations_with_replacement
>>> from fractions import Fraction as frac
>>> #  This is after conditioning on the sum of the counts.
... def probability_of_specific_histogram_(binprobabilities,counts):
...     productsofar = 1
...     for n in counts:
...             productsofar = productsofar*factorial(n)
...     multinomialcoefficient = (factorial(sum(counts)))//productsofar
...     productsofar = frac(1,1)
...     for (p,n) in zip(binprobabilities,counts):
...             productsofar = productsofar*(p**n)
...     return multinomialcoefficient*productsofar
... 
>>> CDFfornumofhighrolls = []
>>> #  This loop results in CDFfornumofhighrolls using between 1GB and 2GB, and takes my laptop approximately 6 minutes.
... for sumofcounts in range(56):
...     distribforthissumofcounts = dict()
...     for cumulativecount in combinations_with_replacement(range(1+sumofcounts),4):
...             counts = (cumulativecount[0],cumulativecount[1]-cumulativecount[0],cumulativecount[2]-cumulativecount[1],cumulativecount[3]-cumulativecount[2],sumofcounts-cumulativecount[3])
...             p = probability_of_specific_histogram_((frac(5,15),frac(4,15),frac(3,15),frac(2,15),frac(1,15)),counts)
...             distribforthissumofcounts[p] = distribforthissumofcounts.get(p,frac(0,1))+p
...     pairs = sorted(distribforthissumofcounts.items())
...     cumulativeprob = pairs[0][1]
...     for i in range(1,len(pairs)):
...             x,probabilitymass = pairs[i]
...             cumulativeprob = cumulativeprob+probabilitymass
...             pairs[i] = (x,cumulativeprob)
...     CDFfornumofhighrolls.append(pairs)
... 
>>> min(len(entry) for entry in CDFfornumofhighrolls)
1
>>> def search_sorted_dictionary_(needle,sortedhaystack):
...     if needle <= sortedhaystack[0][0]:
...             if needle == sortedhaystack[0][0]:
...                     return (0,0)
...             else:
...                     return (-1,0)
...     low,high = 0,len(sortedhaystack)-1
...     if sortedhaystack[high][0] <= needle:
...             if sortedhaystack[high][0] == needle:
...                     return (high,high)
...             else:
...                     return (high,high+1)
...     while low+1 != high:
...             mid = (low+high)//2
...             if sortedhaystack[mid][0] == needle:
...                     return (mid,mid)
...             if sortedhaystack[mid][0] < needle:
...                     low = mid
...             else:
...                     high = mid
...     return (low,high)
... 
>>> all(entry[len(entry)-1][1] == frac(1,1) for entry in CDFfornumofhighrolls)
True
>>> def p_value_of_(histogram,digits=20):
...     if not ((type(histogram) == tuple) and (type(digits) == int)):
...             raise TypeError
...     if not ((len(histogram) == 11) and (0 <= digits < 999)):
...             raise ValueError
...     if not all(type(entry) == int for entry in histogram):
...             raise TypeError
...     if not (all(0 <= entry for entry in histogram) and (sum(histogram) < 56)):
...             raise ValueError
...     threshold = probability_of_specific_histogram_((frac(1,36),frac(2,36),frac(3,36),frac(4,36),frac(5,36),frac(6,36),frac(5,36),frac(4,36),frac(3,36),frac(2,36),frac(1,36)),histogram)
...     lowcumulator,highcumulator = frac(0,1),frac(0,1)
...     sumofhistogram = sum(histogram)
...     for cumulativecount in combinations_with_replacement(range(1+sumofhistogram),6):
...             coarsecounts = (cumulativecount[0],cumulativecount[1]-cumulativecount[0],cumulativecount[2]-cumulativecount[1],cumulativecount[3]-cumulativecount[2],cumulativecount[4]-cumulativecount[3],cumulativecount[5]-cumulativecount[4],sumofhistogram-cumulativecount[5])
...             probabilityofcoarsecount = probability_of_specific_histogram_((frac(1,36),frac(2,36),frac(3,36),frac(4,36),frac(5,36),frac(6,36),frac(15,36)),coarsecounts)
...             relevantCDF = CDFfornumofhighrolls[coarsecounts[6]]
...             low,high = search_sorted_dictionary_(threshold/probabilityofcoarsecount,relevantCDF)
...             if low == -1:
...                     continue
...             elif high == len(relevantCDF):
...                     lowcumulator,highcumulator = lowcumulator+probabilityofcoarsecount,highcumulator+probabilityofcoarsecount
...             elif low == high:
...                     highcumulator = highcumulator+(probabilityofcoarsecount*relevantCDF[high][1])
...                     if low != 0:
...                             lowcumulator = lowcumulator+(probabilityofcoarsecount*relevantCDF[low-1][1])
...             else:
...                     contribution = probabilityofcoarsecount*relevantCDF[low][1]
...                     lowcumulator,highcumulator = lowcumulator+contribution,highcumulator+contribution
...     if not (frac(0,1) <= lowcumulator < highcumulator <= frac(1,1)):
...             raise RuntimeError
...     lowstring = '0.'+(str(floor((10**digits)*lowcumulator)).zfill(digits))
...     if highcumulator == frac(1,1):
...             highstring = '1.'+('0'*digits)
...     else:
...             highstring = '0.'+(str(ceil((10**digits)*highcumulator)).zfill(digits))
...     return '[ '+lowstring+' , '+highstring+' ]'
... 
>>> percentages = (0,15,4,20,15,5,7,11,5,15,4)
>>> len(percentages)
11
>>> sum(percentages)
101
>>> max(abs(frac(p*55,100)-round(frac(p*55,100))) for p in percentages)
Fraction(1, 4)
>>> histogram = tuple(round(frac(p*55,100)) for p in percentages)
>>> sum(histogram)
55
>>> p_value_of_(histogram,10)    #  This line takes my laptop approximately 35 minutes.
'[ 0.0018128154 , 0.0018130603 ]'
>>> frac(1,552) < frac(18128,10**7) < frac(18131,10**7) < frac(1,551)
True
>>> 

has my code for those bounds. ​ The code is long and certainly not simple, so I also checked it by simulation:

Code: Select all

Python 3.8.0 (default, 
[GCC 8.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from math import factorial
>>> from random import randrange
>>> from fractions import Fraction as frac
>>> def probability_of_specific_histogram_(binprobabilities,counts):
...     productsofar = 1
...     for n in counts:
...             productsofar = productsofar*factorial(n)
...     multinomialcoefficient = (factorial(sum(counts)))//productsofar
...     productsofar = frac(1,1)
...     for (p,n) in zip(binprobabilities,counts):
...             productsofar = productsofar*(p**n)
...     return multinomialcoefficient*productsofar
... 
>>> def from_fiftyfive_rolls_():
...     counts = [0 for _ in range(11)]
...     for roll in range(55):
...             counts[randrange(6)+randrange(6)] += 1
...     return tuple(counts)
... 
>>> fortwodice = tuple(frac(6-abs(roll-7),36) for roll in range(2,13))
>>> fortwodice
(Fraction(1, 36), Fraction(1, 18), Fraction(1, 12), Fraction(1, 9), Fraction(5, 36), Fraction(1, 6), Fraction(5, 36), Fraction(1, 9), Fraction(1, 12), Fraction(1, 18), Fraction(1, 36))
>>> def simulation_based_p_value_(histogram,numberofsimulations):
...     threshold = probability_of_specific_histogram_(fortwodice,histogram)
...     lowcumulator,highcumulator = 0,0
...     for _ in range(numberofsimulations):
...             p = probability_of_specific_histogram_(fortwodice,from_fiftyfive_rolls_())
...             if p <= threshold:
...                     highcumulator += 1
...                     if p != threshold:
...                             lowcumulator += 1
...     return (lowcumulator/numberofsimulations,highcumulator/numberofsimulations)
... 
>>> percentages = (0,15,4,20,15,5,7,11,5,15,4)
>>> len(percentages)
11
>>> sum(percentages)
101
>>> max(abs(frac(p*55,100)-round(frac(p*55,100))) for p in percentages)
Fraction(1, 4)
>>> histogram = tuple(round(frac(p*55,100)) for p in percentages)
>>> sum(histogram)
55
>>> simulation_based_p_value_(histogram,10**7)    #  This line takes my laptop approximately 11 minutes.
(0.0018246, 0.0018252)
>>> 

[frogstar_A]

Can you perhaps explain what all that coding that means?!

[Gulchen]

I have started slowly putting together an explanation for that.

[Gulchen]

​
warmup:

Code: Select all

Let ROLLS be {2,3,4,...,11,12}, and consider the distribution
on ROLLS given by    roll 2d6 and add them together  .
That distribution results in each element of ROLLS having a probability.


For this warmup, I can easily list those probabilities:

2 -> 1/36    3 -> 1/18    4 -> 1/12    5 -> 1/9    6 -> 5/36    7 -> 1/6
8 -> 5/36    9 -> 1/9    10 -> 1/12    11 -> 1/18    12 -> 1/36


Let  f : ROLLS -> [0,1]  be given by
f(r) is the probability which that distribution assigns to r  .
For example,    f(2) = 1/36    and    f(6) = 5/36  .

This gives a way to compare the surprisingness of at least some pairs of elements of ROLLS:
For elements x and y or ROLLS, if  f(x) < f(y)  then x is more surprising than y.


Thus, for an element y of ROLLS, the p-value of y is

at least the probability of  f(x) < f(y)
for x chosen according to the distribution on ROLLS
and
at most the probability of  f(x) <= f(y)
for x chosen according to the distribution on ROLLS

.


To calculate those bounds, one needs the distribution of
f(x) for x chosen according to the distribution on ROLLS  .

abstractly:

Code: Select all

Let PAIRS be a non-empty set of ordered pairs, let Laxis be the set of
left entries of PAIRS, left Raxis be the set of right entries of PAIRS,
let DIST be a probability distribution on PAIRS that gives each
element nonzero probability, and let  f : PAIRS -> (0,1]  be given by
f((x,y)) is the probability that DIST assigns to (x,y)  .


For an element (c,d) of PAIRS, one wants

the probability of  f((a,b)) < f((c,d))
and
the probability of  f((a,b)) <= f((c,d))

,  each for (a,b) chosen according to DIST.


Let LDIST be the distribution on Laxis given by
x for (x,y) chosen according to DIST
,  and for each L in Laxis, let RDIST(L) be the distribution on Raxis given by
y for (x,y) chosen according to DIST, _conditioned on_ x=L.
https://en.wikipedia.org/wiki/Conditional_probability

By construction, DIST is equal to the distribution on PAIRS given by
choose L according to LDIST, choose R according to RDIST(L), output (x,y).

For each L in Laxis and R in Raxis, let Lprob(L) be the probability LDIST
assigns to L, and let Rprob(L,R) be the probability that RDIST(L) assigns to R.

By the equality from two sentences ago, one has
f((L,R))  =  Lprob(L) * Rprob(L,R)
for all L in Laxis and R in Raxis.



Now, for each element (c,d) of PAIRS, one has


the probability of  f((a,b)) < f((c,d))    for (a,b) chosen according to DIST

				      =

		   the probability of  f((L,R)) < f((c,d))
     for L chosen according to LDIST and R chosen according to RDIST(L)

				      =

	  the probability of    Lprob(L) * Rprob(L,R)  <  f((c,d))
     for L chosen according to LDIST and R chosen according to RDIST(L)

				      =

		  the sum, for L in Laxis, of  [ Lprob(L) *
	  [the probability of    Lprob(L) * Rprob(L,R)  <  f((c,d))
		   for R chosen according to RDIST(L)  ] ]

				      =

		  the sum, for L in Laxis, of  [ Lprob(L) *
	  [the probability of    Rprob(L,R)  <  f((c,d)) / Lprob(L)
		   for R chosen according to RDIST(L)  ] ]


,  where the division for the last step is valid because DIST gives each
element of PAIRS nonzero probability, so Lprob(L) is strictly positive.



For the  f((a,b)) <= f((c,d))  part, just
replace the    < s    with    <= s  .


For R chosen according to RDIST(L),

the probability of    Rprob(L,R)  <  f((c,d)) / Lprob(L)
and
the probability of    Rprob(L,R)  <=  f((c,d)) / Lprob(L)

come from the cumulative distribution function for
Rprob(L,R)  for R chosen according to RDIST(L).
https://en.wikipedia.org/wiki/Cumulative_distribution_function


(Just like in the warmup, the relevant thing here is the distribution
of the _probability_ Rprob(L,R), not the distribution of R itself.)

What has come so far, is just a huge rearrangement of terms,
and on its own would not speed up the computation.

However, the key point is that, for analyzing histograms, although RDIST(L) does
depend on L, that dependency is only on a small amount of information about L:
Given that information, RDIST(L) is independent of the rest of L.


Thus, one can compute a comparatively small number of CDFs
https://en.wikipedia.org/wiki/Cumulative_distribution_function
,  so that the subsequent computation becomes

the sum, for L in Laxis, of
Lprob(L) * the value of the relevant CDF at  f((c,d)) / Lprob(L)

,  or with the CDF evaluated at a slightly smaller input,
when one needs to be careful about  <  vs  <= .

abstract to concrete:

Code: Select all

Laxis is the set of functions  f : {2,3,4,5,6,7} -> {0,1,2,3,...,54,55}
such that the sum of the outputs is at most 55.

Raxis is the set of functions  f : {8,9,10,11,12} -> {0,1,2,3,...,54,55}
such that the sum of the outputs is at most 55.

PAIRS is the subset of the Cartesian product  Laxis x Raxis
https://en.wikipedia.org/wiki/Cartesian_product
consisting of the pairs where the sum of all 11 outputs is exactly 55.


DIST is the distribution on PAIRS given by

Roll 2d6 55 times, let  h : {2,3,4,...,11,12 -> {0,1,2,3,...,54,55}
be given by    h(n)  =  the number of times the sum was n
,  and output
(  restriction of h to {2,3,4,5,6,7}  ,  restriction of h to {8,9,10,11,12}  )

.


RDIST(L) depends on the the sum of L's outputs,
but does not depend on anything else about L:

The sum of R's outputs will be 55 minus the sum of L's outputs,
and once those sums are fixed - thereby giving _how many_ rolls
are high and _how many_ rolls are not high - the values of the
high rolls are independent of the values of the non-high rolls.

coding:

Code: Select all

One can sample from LDIST by generating a full
histogram and then dropping the high bins.


Since one knows the original number of rolls - 55 -
this is equivalent to merging the high bins together:

In that case, the new bin has probability    15 / 36  .


The original probabilities of the high bins are
5 / 36  ,  4 / 36  ,  3 / 36  ,  2 / 36  ,  1 / 36
,  so conditioning on landing in one of them gives probabilities
5 / 15  ,  4 / 15  ,  3 / 15  ,  2 / 15  ,  1 / 15    .

Thus, once one knows the sum of the counts of those bins
- this is 55 minus the sum of L's outputs - one can calculate
probability(R given L)  https://en.wikipedia.org/wiki/Conditional_probability
in the same way as one calculates the probability of any other histogram.

These histogram probability calcuations are done by my
Python function  probability_of_specific_histogram_ .


The  for sumofcounts in range(56)  loop computes the CDFs mentioned near the end
of the "abstractly:" section, and pushes them to the CDFfornumofhighrolls list.
Each of these is computed via the following steps:

Compute the  https://en.wikipedia.org/wiki/Probability_mass_function ,
as a  https://en.wikipedia.org/wiki/Dictionary_(data_structure) .
Sort the pairs in that dictionary:  Due to how tuple comparison works,
this sorts by keys, which are inputs to the PMF.
Replace the output entries in that list with their cumulative sums
https://en.wikipedia.org/wiki/Prefix_sum
,  to get the values of the CDF.


Each list in CDFfornumofhighrolls came from a dictionary, so
their encoded inputs are all unique, in addition to being sorted.
This simplifies the binary search code, which I have
in my Python function  search_sorted_dictionary_ .


My Python function p_value_of_ computes the sums from the end of
the "abstractly:" section, by calling search_sorted_dictionary_ on
the relevant list from CDFfornumofhighrolls to evaluate the CDFs:

threshold is f((c,d)) from the "abstractly:" section, and
the part of the function from  "if not (frac(0,1)" onwards
is just to display the results nicely, since lowcumulator and
highcumulator themselves will be fracs with huge numbers of digits.

​

[Wyliet]

The numbers in the game do seem a bit extreme. However, extreme results are normal in some games, and 55 rolls is not that many.

[pdy333]

Are you sure you didnt forgot to check the "luck draws" option in your player panel ? This should appears once you are premium.
But be cautious with this option, if several players around the table has the option enabled, the table explodes since you cant have rolls that are good for all the players.

Where is this “luck draws” option?

[Jellby]

Right next to the "disable sarcasm" option.

[yfimc]

Per example, the chance of rolling at least eight 3 in 56 rolls is 1,18% ( https://docs.google.com/spreadsheets/d/ ... 1846324187 )

If now you do my fast formula 1-((1-0.0118)^11) to try to multiply my formula by 11, you find 12%.
Which mean in 12% of the games you should have a game where a specific number is rolled more than what is expected compared to your number of 3s in that game.


(edit: forgot a 1, 18% => 12%)

thanks this is interesting