To me ~0.29% seems more realistic than 35.3028752%. So assuming the bet was a max bet (which it probably wasn't but was close) the house overbet by a factor of roughly 120x. Good job!
Kelly tells us to maximise the sum of the probabilities times the log of (1 + b.x) where b is the profit divided by the stake, and x is the proportion of the bankroll to risk.
Plotting that function against x I see:
Find the x value where that peaks, and we have the proportion of the bankroll to risk. It looks like it peaks at around 0.3, meaning we risk 30%.
Let's zoom in:
So it peaks very close to 0.353, meaning we risk 35.3% of the bankroll.
Thanks Kelly!
Is your Haskell code remembering that we aren't risking the stake which is returned to the player? On the potential profit (payout - stake) is "risk" to the bankroll.
Edit: I should post the equation of the curve I'm plotting so you can check I didn't do anything stupid. I'm looking at the bet 'backwards', ie. from the house's point of view. I am considering the house bet size to be the most it can lose on the bet. So the 'b' for the jackpot is -1 (a profit of -1 times the bet size). 'b' for the 2nd highest win is -46/120 (we only lose 46/120ths as much as we were willing to risk because the highest two payouts are 121 and 47). Etc:
plot [0.352:0.354] \
1/65536.0 * log(1 + x * -120/120) + \
16/65536.0 * log(1 + x * -46/120) + \
120/65536.0 * log(1 + x * -12/120) + \
560/65536.0 * log(1 + x * -4/120) + \
1820/65536.0 * log(1 + x * -2/120) + \
4368/65536.0 * log(1 + x * -0.4/120) + \
8008/65536.0 * log(1 + x * 0/120) + \
11440/65536.0 * log(1 + x * 0.5/120) + \
12870/65536.0 * log(1 + x * 0.7/120) + \
11440/65536.0 * log(1 + x * 0.5/120) + \
8008/65536.0 * log(1 + x * 0/120) + \
4368/65536.0 * log(1 + x * -0.4/120) + \
1820/65536.0 * log(1 + x * -2/120) + \
560/65536.0 * log(1 + x * -4/120) + \
120/65536.0 * log(1 + x * -12/120) + \
16/65536.0 * log(1 + x * -46/120) + \
1/65536.0 * log(1 + x * -120/120) \
title "maximise this"
Edit2: I tried reading the Haskell code. (I don't know Haskell. I don't even know how many L's it has):
table =
[ ( 65536 / (2^32) , 1 - 121.0 ), ( 749731840 / (2^32) , 1 - 0.5 )
[...]
fun x = sum [ p * log (1 + b * x) | (p,b) <- table ]
It looks like you're using newton's method to maximize the sum of the products - but you're using the payout multipliers from the player's point of view. You need to look at it from the other side. The house never wins 121x. The player's big wins are the house's big losses. We need to calculate each house profit or loss as a multiplier of the amount the house is risking.