CatanOdds

catan-odds.py

#!/usr/bin/env python

import operator as op
import sys
from decimal import Decimal

import decimal
decimal.getcontext().prec = 6

def ncr(n, r):
    r = min(r, n-r)
    if r == 0: return 1
    numer = reduce(op.mul, xrange(n, n-r, -1))
    denom = reduce(op.mul, xrange(1, r+1))
    return numer//denom

def mul(seq):
    return reduce(lambda m,n: n*m, seq, 1)

def possible_outcomes(dice):
    return range(len(dice), sum(dice)+1)

def roll_sum_variants(roll_sum):
    p = 0
    for i in range(1, 7):
        j = roll_sum - i
        if j > 6: continue
        if i <= j:
            p += 2 if i != j else 1
        else:
            break
    return p

def main(roll_sum, roll_counts, roll_total):
    x = roll_sum_variants(roll_sum)
    p = 0
    for roll_count in roll_counts:
        p += mul([
            ncr(roll_total, roll_count),
            pow(Decimal(x) / 36, roll_count),
            pow(Decimal(36 - x) / 36, roll_total - roll_count)
        ])

    return p

if __name__ == '__main__':
    rolls = {
        2: 1,
        3: 5,
        4: 5,
        5: 10,
        6: 11,
        7: 16,
        8: 4,
        9: 6,
        10: 3,
        11: 4,
        12: 6
    }
    roll_total = sum(rolls.values())

    print("After %d rolls using 2d6:" % roll_total)
    for roll_sum in sorted(rolls.keys()):
        roll_count = rolls[roll_sum]
        rsv = roll_sum_variants(roll_sum)

        expected = int(Decimal(rsv) / 36 * roll_total)
        if abs(expected - roll_count) > 3: #todo: 3 is magic, how to make dynamic?
            print(" ! Rolled it %d times, expected ~%d times" % (
                roll_count, expected
            ))
            if expected < roll_count:
                print "   + P(%d) = %0.1f%% [%d / 36]; Rolling that %d or more times is only %0.3f%%" % (
                    roll_sum, 100 * Decimal(rsv) / 36, rsv,
                    roll_count, 100 * main(
                        roll_sum,
                        range(roll_count, roll_total + 1),
                        roll_total
                    )
                )
            else:
                print "   - P(%d) = %0.1f%% [%d / 36]; Rolling that %d or less times is only %0.3f%%" % (
                    roll_sum, 100 * Decimal(rsv) / 36, rsv,
                    roll_count, 100 * main(
                        roll_sum,
                        range(roll_count + 1),
                        roll_total
                    )
                )

$ ./catan-odds.py
After 71 rolls using 2d6:
 ! Rolled it 16 times, expected ~11 times
   + P(7) = 16.7% [6 / 36]; Rolling that 16 or more times is only 12.335%
 ! Rolled it 4 times, expected ~9 times
   - P(8) = 13.9% [5 / 36]; Rolling that 4 or less times is only 2.388%
 ! Rolled it 6 times, expected ~1 times
   + P(12) = 2.8% [1 / 36]; Rolling that 6 or more times is only 1.413%