#!/usr/bin/env python
"""
 Craps simulation
 See for example http://en.wikipedia.org/wiki/Craps
 
 Chap 9, Exercise 7, pg 292 of Zelle's "Python Programming"

 Turns out that to get enough statistical power to really tell
 if you're more likely to win or lose you need around a million games.

 Jim Mahoney, Oct 24 2006 for the Intro Programming class

 --- output (with tests and all defaults) ----
 
   $ ./craps.py 
   Testing 5 games 

   === game 1 ===
       first roll is 11
       which wins outright.

   === game 2 ===
       first roll is 10
       roll 2 is 10
       which makes the point and wins.

   === game 3 ===
       first roll is 4
       roll 2 is 12
       roll 3 is 10
       roll 4 is 7
       which is a seven-out and loses.

   === game 4 ===
       first roll is 8
       roll 2 is 11
       roll 3 is 10
       roll 4 is 9
       roll 5 is 8
       which makes the point and wins.

   === game 5 ===
       first roll is 7
       which wins outright.

   Simulating 1e+04 games 
   and calculating the odds of winning at craps...
       number of games simulated  =    10000
       number of games won        =     4898
       odds of winning are        =   0.4898 +- 0.0100 at 95% confidence

 --- output (no tests, more simulated games) ----

   $ ./craps.py 
   Simulating 1e+06 games 
   and calculating the odds of winning at craps...
     number of games simulated  =  1000000
     number of games won        =   491834
     odds of winning are        =   0.4927 +- 0.0010 at 95% confidence

"""

from math import sqrt
from random import randint

def roll():
    """ Return the result of rolling one 6-sided die.       """
    return randint(1,6)

def game(verbose=False):
    """ Play a game of craps, 
        returning True or False for a win or a loss.
        If 'verbose' is set, print out a description as well. """
    first_roll = roll() + roll()
    if (verbose):
        print "  first roll is", first_roll
    if (first_roll in (2,3,12)):                          # We can lose right away ...
        if (verbose):
            print "  which craps out and loses." 
        return False
    elif (first_roll in (7,11)):                          # or win right away ...
        if (verbose):
            print "  which wins outright."
        return True
    else:                                                 # or keep rolling.
        roll_number = 1
        while True:
            roll_number = roll_number + 1
            next_roll = roll() + roll()
            if (verbose):
                print "  roll", roll_number, "is", next_roll
            if (next_roll == first_roll):
                if (verbose):
                    print "  which makes the point and wins."
                return True
            elif (next_roll == 7):
                if (verbose):
                    print "  which is a seven-out and loses."
                return False

def test_games(n_games=5):
    """ Run a bunch of verbose games  """
    print "Testing", n_games, "games "
    print
    for i in range(n_games):
        print "=== game", i+1, "==="
        game(verbose=True)
        print

def calculate_odds(n_games=1e5, verbose=True):
    """ Calculate the odds of winning by running a number of games.
        Return the measured odds as the fraction wins/n_games.
        If 'verbose' is True (the default), print the odds too. """
    if (verbose):
        print "Simulating %.0e games " % n_games
        print "and calculating the odds of winning at craps..."
    (wins, games) = (0, 0)
    n_games = int(n_games)
    while (games < n_games):
        games = games + 1
        if (game()):
            wins = wins + 1
    odds = 1.0 * wins / n_games
    # sigma_1 = estimated standard deviation of odds of a single game
    # sigma_n = estimated standard deviation of odds from n games
    # statistics say true value is in +-2*sigma with 95% probability.
    sigma_1 = sqrt( odds * (1.0 - odds) )
    sigma_n = sigma_1 / sqrt( 1.0 * n_games )
    if (verbose):
        print "  number of games simulated  = %8i" % n_games
        print "  number of games won        = %8i" % wins
        print "  odds of winning are        = " + \
              "%8.4f +- %.4f at 95%% confidence" % (odds, 2*sigma_n)
    return odds

test_games()        # Comment this out to avoid the verbose test runs.
calculate_odds()    # Pass in 1e6 or so for statistics to show odds < 0.5