# jackpot amounts (as of 12/16/2011)
# http://www.lottostrategies.com/script/jackpot_history/draw_date/101
# http://www.usamega.com/powerball-jackpot.asp
jackpots <- read.csv('http://biostatmatt.com/csv/jackpot.csv')
average_jackpot <- mean(jackpots$PrizeMillions)

# for convenience
ch <- choose

# old rules
# choose 5 from 59 white balls, 1 from 39 red balls

# probabilities of winning
powerball_old <- function(w, r)
    ch(5,w)*ch(54,5-w)/ch(59,5)*ch(1,r)*ch(38,1-r)/ch(39,1)

# 9 ways to win
p_old <- vector("numeric", length=9)
p_old[1] <- powerball_old(5, 1) # five white + powerball_old
p_old[2] <- powerball_old(5, 0) # five white
p_old[3] <- powerball_old(4, 1) # four white + powerball_old
p_old[4] <- powerball_old(4, 0) # four white
p_old[5] <- powerball_old(3, 1) # three white + powerball_old
p_old[6] <- powerball_old(3, 0) # three white
p_old[7] <- powerball_old(2, 1) # two white + powerball_old
p_old[8] <- powerball_old(1, 1) # one white + powerball_old
p_old[9] <- powerball_old(0, 1) # powerball_old

# winnings
w_old <- vector("numeric", length=9)
w_old[1] <- 1000000 * average_jackpot
w_old[2] <- 200000
w_old[3] <-  10000
w_old[4] <-    100
w_old[5] <-    100
w_old[6] <-      7
w_old[7] <-      7
w_old[8] <-      4
w_old[9] <-      3

# expected winnings (loss)
# cost to play is $1
expected_winnings_old_rules <- -1 * (1 - sum(p_old)) + sum(p_old * (w_old - 1))

# new rules
# choose 5 from 59 white balls, 1 from 35 red balls
# probabilities of winning powerball
powerball_new <- function(w, r)
    ch(5,w)*ch(54,5-w)/ch(59,5)*ch(1,r)*ch(34,1-r)/ch(35,1)

# 9 ways to win
p_new <- vector("numeric", length=9)
p_new[1] <- powerball_new(5, 1) # five white + powerball_new
p_new[2] <- powerball_new(5, 0) # five white
p_new[3] <- powerball_new(4, 1) # four white + powerball_new
p_new[4] <- powerball_new(4, 0) # four white
p_new[5] <- powerball_new(3, 1) # three white + powerball_new
p_new[6] <- powerball_new(3, 0) # three white
p_new[7] <- powerball_new(2, 1) # two white + powerball_new
p_new[8] <- powerball_new(1, 1) # one white + powerball_new
p_new[9] <- powerball_new(0, 1) # powerball_new

# winnings
# jackpots start $20M larger than before
# second prize is $1M instead of $0.2M
# prize for powerball only is $4 instead of $3
w_new <- vector("numeric", length=9)
w_new[1] <- 1000000 * (average_jackpot + 20)
w_new[2] <- 1000000
w_new[3] <-  10000
w_new[4] <-    100
w_new[5] <-    100
w_new[6] <-      7
w_new[7] <-      7
w_new[8] <-      4
w_new[9] <-      4

# expected winnings (loss)
# cose to play is $2 instead of $1
expected_winnings_new_rules <- -2 * (1 - sum(p_new)) + sum(p_new * (w_new - 2))