Tentative Lotto Odds Calculations

from __future__ import division
from math import *
import random
import os

# Pretty

def MaxDP(x, n):
  s = '%.*f' % (n, x)
  if '.' in s:
    while s[-1:] == '0':
      s = s[:-1]
    if s[-1:] == '.':
      s = s[:-1]
  return s

# n choose k

def C(n, k):
  if k >= 0 and k <= n:
    if k == 0 or k == n:
      Result = 1
    else:
      if n < 4:
        Result = n
      else:
        if k <= (n // 2):
          j = k
        else:
          j = n - k
        Numerator = 1
        Divisor = 1
        i = 1
        while i <= j:
          Numerator *= (n - j + i)
          Divisor *= i
          i += 1
        Result = Numerator / Divisor
  else:
    Result = 0
  return Result

# Permutations of k items drawn from n items

def P(n, k):
  Result = 1
  Limit = n - k
  x = n
  while x > Limit:
    Result *= x
    x -= 1
  return Result

# Combinations of k items drawn from n items which include
# exactly j specified matches.

def J(n, k, j):
  return C(n - j, k - j)

# (Faulty) Lotto odds

def B(n, k, j, b):
  return (C(k, j) * C(1, b) * C(n - k - 1, k - j - b)) / C(n, k)

# >>> 1/B(40, 6, 6, 0)
# 3838380.0
# >>> 1/B(40, 6, 5, 1)
# 639730.0
# >>> 1/B(40, 6, 5, 0)
# 19385.757575757576
# >>> 1/B(40, 6, 4, 1)
# 7754.303030303031
# >>> 1/B(40, 6, 4, 0)
# 484.64393939393943
# >>> 1/B(40, 6, 3, 1)
# 363.4829545454545
# >>> 1/B(40, 6, 3, 0)
# 35.17576979472141

# Slosh's formula:

def Zb(n, k, j):
  return C(k,j)*C(n-k-1,k-j)/C(n,k)*1/(n-k)

def Zn(n, k, j):
  return C(k,j)*C(n-k-1,k-j)/C(n,k)*(1 - 1/(n-k))

# Simulation

def DieRoll(n):
  return 1 + int(floor(n * random.random()))

def Draw(n, k):
  # A boring sequence is statistically as good as any other to win even if
  # Lotto forbids the customer to make such a selection and risk sharing
  # their prize with many other people who think runs of numbers are great.
  Line = list(range(1, k + 1))
  LineBonus = k + 1
  Bag = list(range(1, n + 1))
  Drawn = []
  for i in range(k):
    Drawn.append(Bag.pop(DieRoll(len(Bag)) - 1))
  Bonus = Bag.pop(DieRoll(len(Bag)) - 1)
  NumMatches = 0
  for x in Line:
    if x in Drawn:
      NumMatches += 1
  return (NumMatches, Bonus == LineBonus)

def WinDist(n, k, Trials):
  H = [(0, 0) for j in range(k + 1)]
  i = 0
  while i < Trials:
    i += 1
    (m, b) = Draw(n, k)
    if b:
      H[m] = (H[m][0], H[m][1] + 1)
    else:
      H[m] = (H[m][0] + 1, H[m][1])
  Eta = 1/Trials
  N = [(Eta * x[0], Eta * x[1]) for x in H]
  return N

# More pretty

def OddsStr(x, DP=None):
  try:
    Odds = 1/x
    if DP is None:
      DP = max(0, 2 - int(floor(log10(Odds))))
    return '1:' + MaxDP(Odds, DP)
  except (ZeroDivisionError):
    return 'None'

def OddsH(H, DP=None):
  return [(OddsStr(x[0], DP), OddsStr(x[1], DP)) for x in H]

def PrintLottoH(H):
  PLose = sum(H[0]) + sum(H[1]) + sum(H[2]) + H[3][0]
  DivRecs = {
    1: (H[6][0] + H[6][1], '6 balls'),
    2: (H[5][1], '5 balls + bonus'),
    3: (H[5][0], '5 balls, no bonus'),
    4: (H[4][1], '4 balls + bonus'),
    5: (H[4][0], '4 balls, no bonus'),
    6: (H[3][1], '3 balls + bonus')
  }
  PWin = 0
  for d in sorted(DivRecs):
    (p, Desc) = DivRecs[d]
    PWin += p
    print 'Division %d (%s): %s' % (d, Desc, OddsStr(p))
  print 'Any division: ' + OddsStr(PWin)
  print 'P(No luck): ' + str(PLose)

# >>> H = WinDist(40, 6, 50000000)
# >>> PrintLottoH(H)
# Division 1 (6 balls): 1:7142857
# Division 2 (5 balls + bonus): 1:617284
# Division 3 (5 balls, no bonus): 1:19904
# Division 4 (4 balls + bonus): 1:16573
# Division 5 (4 balls, no bonus): 1:466
# Division 6 (3 balls + bonus): 1:1198
# Any division: 1:323
# P(No luck): 0.99690874
# >>> H
# [(0.34190612, 0.00847088), (0.42398538, 0.01088244), (0.17660452000000001, 0.00470074), (0.03035866, 0.00083488), (0.00214404, 6.034e-05), (5.024e-05, 1.62e-06), (1.4e-07, 0.0)]
# >>> OddsH(H)
# [('1:2.92', '1:118'), ('1:2.36', '1:91.9'), ('1:5.66', '1:213'), ('1:32.9', '1:1198'), ('1:466', '1:16573'), ('1:19904', '1:617284'), ('1:7142857', 'None')]

# >>> H = WinDist(40, 6, 1000000000)
> >> PrintLottoH(H)
# Division 1 (6 balls): 1:4032258
# Division 2 (5 balls + bonus): 1:621118
# Division 3 (5 balls, no bonus): 1:19390
# Division 4 (4 balls + bonus): 1:16413
# Division 5 (4 balls, no bonus): 1:469
# Division 6 (3 balls + bonus): 1:1197
# Any division: 1:324
# P(No luck): 0.996917296
# >>> H
# [(0.34188522200000004, 0.008488078000000001), (0.424051646, 0.010917859), (0.17652759, 0.00470585), (0.030341051, 0.000835703), (0.002132643, 6.0927e-05), (5.1573e-05, 1.61e-06), (2.43e-07, 5e-09)]
# >>> OddsH(H)
# [('1:2.92', '1:118'), ('1:2.36', '1:91.6'), ('1:5.66', '1:213'), ('1:33', '1:1197'), ('1:469', '1:16413'), ('1:19390', '1:621118'), ('1:4115226', '1:200000000')]

# >>> H1 = [(Zn(40,6,j), Zb(40,6,j)) for j in range(7)]
# >>> OddsH(H1)
# [('1:3.57', '1:118'), ('1:2.78', '1:91.6'), ('1:6.44', '1:213'), ('1:36.2', '1:1196'), ('1:499', '1:16478'), ('1:19973', '1:659116'), ('1:3954695', '1:130504920')]
#
# >>> PrintLottoH(H1)
# Division 1 (6 balls): 1:3838380
# Division 2 (5 balls + bonus): 1:659116
# Division 3 (5 balls, no bonus): 1:19973
# Division 4 (4 balls + bonus): 1:16478
# Division 5 (4 balls, no bonus): 1:499
# Division 6 (3 balls + bonus): 1:1196
# Any division: 1:339
# P(No luck): 0.847048647668