09-14-2011 , 12:07 AM
I think I’ve discovered a new algorithm to accurately calculate ICM equities for large numbers of people. I’ve done some searching and wasn’t able to find any evidence that this has been used before, so I think I’m the first. If someone else has come up with this before, then my apologies…

First, why is this important? Traditional ICM calculations require you to do a number of calculations proportional to the factorial (n!) of the number of players remaining in a tournament. This is why you never see ICM calculators that can handle more than 10 players – it can take hours to calculate for 15 people and days to years to calculate more than that. My new algorithm can produce very accurate numbers in less than a minute, even if you have over 50 players.

Software developers may find this useful and do what was previously considered impossible, perhaps incorporating some concepts around the bubble factor from Kill Everyone. I will personally be doing some research of my own now that this door is opened.

This algorithm works on the principle of ICM where you imagine that all of the chips in the tournament are randomly given numbers and each player’s place in the tournament is based on their highest scoring chip. However we can use some math tricks to make the calculations easier. If each chip is given a random number between 0 and 1, and we have n chips, then the probability that all of our chips will score x or less will be x^n. We can also scale things so that we can work with smaller stacks since fractional chips are okay in this formula. This has the advantage that we only have to generate 1 random number for each player rather than 1 number for each chip.

Here is the algorithm:

Given n players
Given stack sizes S = S1, S2, … Sn
Given prizes Z = Z1, Z2, … Zn
Get normalized stack sizes Q = Q1, Q2, …Qn by dividing S by average stack

Simulate 1 tournament:
For each player i, generate a “place” for each player (Pi) by generating a random number between 0 and 1 and raising it to the power of 1/Qi. Pi = rand ^ (1.0 / Qi)
Sort players by their place P – higher is better
Award prizes Z according to P

Repeat for as many tournaments as desired and calculate the average prize value for each player. The error will be proportional to 1/sqrt(iterations) so the more iterations you run, the more accurate it will be.

I’ve tested this out and compared them to traditional ICM calculations and they converge to the same numbers.

Tysen
09-14-2011 , 05:18 AM
Interesting approach. Have you thought about solving the place % analytical, instead of sampling?

Fwiw, with some optimizations it is perfectly possible to calculate up to ~30 players of ICM in reasonable time, 20 players / payouts can be done within seconds. e.g: http://www.holdemresources.net/hr/sn...alculator.html

Edit: Meh, i guess we just get something similar to the original ICM algorithm when solving analytical.

Last edited by plexiq; 09-14-2011 at 05:40 AM.
09-14-2011 , 08:17 AM
Nice work!

Do you have a name for the algorithm? I have implemented this in my poker library but don't know what to call it.

I'm mostly into Heads Up now, but if/when I return to MTTs (probably 16-180 man) this will be very useful. Thanks!
09-14-2011 , 11:49 AM
Yeah, I've been thinking of a good name but haven't come up with anything better than Simulated Independent Chip Model or SICM. Or maybe I should try to attch my name to it by calling it Streib's ICM. I guess that's still SICM...
09-14-2011 , 02:16 PM
I like the name, Sic'em.
09-14-2011 , 05:10 PM
^nah... the tysen model
09-14-2011 , 05:24 PM
nice level OP, they beleive you are serious with this. If this isn't a level... next semester my brother is teaching a statistics class at Boston College. PM me if you need deatils, as I think these forums should offer the most assitance they can. Seriously, I want to help as many ppl as I can on the forums
09-14-2011 , 08:15 PM
Quote:
Originally Posted by raiseya
nice level OP, they beleive you are serious with this. If this isn't a level... next semester my brother is teaching a statistics class at Boston College. PM me if you need deatils, as I think these forums should offer the most assitance they can. Seriously, I want to help as many ppl as I can on the forums
k lol is this a level post yourself? OP is author of Kill Everyone and Raisers Edge!
09-14-2011 , 08:41 PM
With a couple of very simple test cases, I do see identical results between your new algorithm and the traditional algorithm.

Can you explain more how you derive this approach? I understand this:

Quote:
the probability that all of our chips will score x or less will be x^n
But I can't see how you get from there to your algorithm.
09-14-2011 , 09:10 PM
A nerdy comment.

Since the simulation will likely involve many thousands (millions?) of calculations, it will be a bit faster by noting that 1/Qi is equal to Average Stack/Si = Vi, say. Then replace

rand^(1/Qi) with rand^Vi
09-15-2011 , 04:15 AM
That stack-normalization is not necessary for the correctness, using rand^(1/stack) will give just the same results in theory.

In practice, the normalization is probably a good idea though. It does not make any difference to runtime, you can simply calculate the final exponent for every stack once during initialization.
09-15-2011 , 05:41 AM
Here is a quick & dirty Java implementation of the algorithm.

Code:
```import java.util.Arrays;
import java.util.Comparator;
import java.util.Random;

public class SimulatedIcm {

private Random rnd = new Random();

public double[] getEquities(double[] payouts, double[] stacks,
int maxIterations) {
if (payouts.length > stacks.length)
payouts = Arrays.copyOf(payouts, stacks.length);

double[] equities = new double[stacks.length];
double[] exp = getExponents(stacks);
final double[] r = new double[stacks.length];

Integer[] ids = new Integer[stacks.length];
for (int i = 0; i < ids.length; i++)
ids[i] = i;

Comparator<Integer> comp = new Comparator<Integer>() {
@Override
public int compare(Integer i1, Integer i2) {
if (r[i1] < r[i2])
return 1;
if (r[i1] > r[i2])
return -1;
return 0;
}
};

for (int iteration = 0; iteration < maxIterations; iteration++) {
for (int i = 0; i < r.length; i++)
r[i] = Math.pow(rnd.nextDouble(), exp[i]);
Arrays.sort(ids, comp);
for (int i = 0; i < payouts.length; i++)
equities[ids[i]] += payouts[i];
}

for (int i = 0; i < equities.length; i++)
equities[i] /= (double) maxIterations;

return equities;
}

private static double[] getExponents(double[] stacks) {
double t = 0;
for (double s : stacks)
t += s;
t /= (double) stacks.length;
double[] exp = new double[stacks.length];
for (int i = 0; i < stacks.length; i++)
exp[i] = t / stacks[i];
return exp;
}

public static void main(String[] args) {
SimulatedIcm sicm = new SimulatedIcm();

double[] payouts = new double[] { 50, 30, 20 };
double[] stacks = new double[] { 1000, 2000, 3000, 4000, 5000, 6000,
7000, 8000, 9000 };
double[] eqs = sicm.getEquities(payouts, stacks, 1000000);
System.out.println(Arrays.toString(eqs));
}
}```
That is entirely unoptimized, obviously. Especially if the number of payouts is significantly lower than the # of players, the sorting step can be done a lot more efficiently.

Last edited by plexiq; 09-15-2011 at 05:48 AM.
09-15-2011 , 06:03 AM
Quote:
Originally Posted by trojanrabbit
I’ve tested this out and compared them to traditional ICM calculations and they converge to the same numbers.
Can you prove it, or at least explain why the result would be expected?

Edit: Took about 20 minutes to cook up a Python implementation

Code:
```from random import random

def trial(weights, payouts):
scores = sorted((random() ** weight, i) for i, weight in enumerate(weights))
results = [0] * len(payouts)
for payout, score in zip(payouts, scores): results[score[1]] = payout
return results

def sicm(stacks, payouts, trials):
avg = sum(stacks) / float(len(stacks))
payouts += [0] * (len(stacks) - len(payouts))
payouts.sort()
weights = [avg / s for s in stacks]
return [sum(player) / float(trials) for player in zip(*(
trial(weights, payouts) for i in xrange(trials)
))]

# I guess I should add an equivalent test to be fair to the Java implementation...
if __name__ == '__main__':
print sicm([1000 * n for n in range(1, 10)], [50, 30, 20], 1000000)```

Last edited by halftilt; 09-15-2011 at 06:31 AM.
09-15-2011 , 07:04 AM
Damn, i need to learn Python sometime. Implementation is like 4x slower though :P
09-15-2011 , 04:36 PM
Yeah, python = optimizing for elegance ldo. I could probably get it faster using Numpy or something though.
09-17-2011 , 08:48 PM
Quote:
Originally Posted by plexiq
That stack-normalization is not necessary for the correctness, using rand^(1/stack) will give just the same results in theory.

In practice, the normalization is probably a good idea though. It does not make any difference to runtime, you can simply calculate the final exponent for every stack once during initialization.
Correct, you don't need to do the normalization step, but I included it so that you don't get any problems trying to raise a small rand to a power of thousands or millions.

Tysen
09-17-2011 , 09:16 PM
Quote:
Originally Posted by halftilt
Can you prove it, or at least explain why the result would be expected?
It's kind of hard to explain, but let me take a crack at it. Okay, one of the properties of ICM is that the odds that any player will finish ahead of any other player is proportional to their stacks. Let's say that you and I want to make a last-longer bet, I have a stack of 20k and you have 10k. Skill aside, the odds that I will place higher than you are 2:1 in my favor. We could simulate this by having you generate 1 random number and I generate 2. Whoever gets the highest random number wins. But I don't want to generate 1 random number for each unit of stack, I just want 1 per player, so I use my formula.

Why does my formula work? I don't want to get into calculus, so I'll try to just illustrate an example and hope it makes sense. The probability that you with your 1 random number gets x or less is simply x. 50% of the time you'll have 0.5 or less and 25% of the time you'll have 0.25 or less. If I get 2 random numbers, the chance I'll have x or less is x^2. So there's only a 25% chance both my numbers are below 0.5, a 9% chance both my numbers are below 0.3, etc. See how that' s like probability distribution that follows y = x^2? Except that it's backwards. I'm picking a random number y and then saying what score (x) gives me y probability of getting that score? So to convert my probability distribution I take the square root of both sides and get sqrt(y) = x. So I take the square root of my random number and get my score x.

You can see this holds true in my last-longer bet situation. Say in a specific example you generate 1 random number and get a 0.300. If I generate 2 numbers, you only have a 9% chance of beating me (both my numbers have to be below 0.3, which is 0.3*0.3 = 9%). Or I could just generate 1 random number and take the square root. If that's >0.3 then I win. And now you can see that as long as the number I generate is at least 0.09 then after I take the square root then I'll have a higher "score."

Hope that's clear enough.

Tysen
09-18-2011 , 04:08 AM
Quote:
Originally Posted by trojanrabbit
It's kind of hard to explain, but let me take a crack at it. Okay, one of the properties of ICM is that the odds that any player will finish ahead of any other player is proportional to their stacks. Let's say that you and I want to make a last-longer bet, I have a stack of 20k and you have 10k. Skill aside, the odds that I will place higher than you are 2:1 in my favor. We could simulate this by having you generate 1 random number and I generate 2. Whoever gets the highest random number wins. But I don't want to generate 1 random number for each unit of stack, I just want 1 per player, so I use my formula.

Why does my formula work? I don't want to get into calculus, so I'll try to just illustrate an example and hope it makes sense.
Actually, I basically just read this far and kinda mentally pictured the hand-wavy calculus that fills in the rest.

Actually formalizing it seems kinda painful, though. In particular, it's hard to go from showing that it works for each pair of players to showing that it works for the whole set of players. Maybe some kind of freaky mathematical induction? :s Anyway, it seems to work, so cheers 8)
09-20-2011 , 01:33 PM
Here is an implementation in R, with reproducible examples, as well. It is worth noting that it always returns a warning message about the anticipated errors in the ICM figures.

Code:
```SICM <- function(stacks, prizes, sims=1000) {
if(length(prizes) != length(stacks)) {stop("Number of prizes not equal to number of remaining stacks. You might need to add prizes of value 0.")}
normstacks <- 1/(stacks / mean(stacks))
res <- matrix(nrow=sims, ncol=length(stacks))
for(i in 1:sims) {
rnums <- runif(length(stacks))
places <- rnums^normstacks
res[i,] <- prizes[length(places)+1 - rank(places, ties.method="random")]
}
out <- data.frame("Start.Stack"=stacks, "ICM.EV"=colMeans(res))
warning(paste("Errors are proportional to", 1/sqrt(sims)))
return(out)
}```
Code:
```stacks <- c(rep(1500,10))
prizes <- c(50, 30, 20, rep(0,7))

# Ex 1
SICM(stacks, prizes)

# Ex 2
SICM(stacks, prizes, sims=10000) # changing the sims for more precise results

# Ex 3
SICM(stacks=c(5000, 3000, 2000), prizes=c(50, 30, 20))```
Code:
```Example output:

> SICM(stacks, prizes)
Start.Stack ICM.EV
1         1500  10.76
2         1500   9.43
3         1500  10.38
4         1500  10.26
5         1500   9.28
6         1500   9.97
7         1500  10.88
8         1500   9.16
9         1500  10.28
10        1500   9.60
Warning message:
In SICM(stacks, prizes) : Errors are proportional to 0.0316227766016838

> SICM(stacks, prizes, sims=10000)
Start.Stack ICM.EV
1         1500  9.838
2         1500  9.976
3         1500 10.060
4         1500 10.212
5         1500  9.542
6         1500  9.923
7         1500 10.652
8         1500  9.593
9         1500 10.113
10        1500 10.091
Warning message:
In SICM(stacks, prizes, sims = 10000) : Errors are proportional to 0.01

> SICM(c(5000, 3000, 2000), c(50, 30, 20), sims=1000)
Start.Stack ICM.EV
1        5000  37.81
2        3000  32.85
3        2000  29.34
Warning message:
In SICM(c(5000, 3000, 2000), c(50, 30, 20), sims = 1000) :
Errors are proportional to 0.0316227766016838```

Last edited by Sherman; 09-20-2011 at 01:47 PM. Reason: added a fail-safe for improper user input.
09-20-2011 , 02:20 PM
nice

is it possible/easy to alter/skew the rand? (to simulate skill differences maybe)
09-20-2011 , 05:33 PM
Oof, looks like the constant factor for the errors is pretty high...
09-20-2011 , 08:28 PM
Quote:
Originally Posted by sputum
nice

is it possible/easy to alter/skew the rand? (to simulate skill differences maybe)
I have always used an adjustment to the stack sizes to represent skill variation. If player 1 has a 10% edge on all the other players, then his adjustment factor is 1.1 while all others remain at 1.0. If player 9 has a -5% edge, than the adjustment factor for him is 0.95.

I suppose non-linear factors can be used but the function would just be a guess, so I stuck with the simple one.
09-21-2011 , 11:55 AM
Quote:
Originally Posted by statmanhal
I have always used an adjustment to the stack sizes to represent skill variation. If player 1 has a 10% edge on all the other players, then his adjustment factor is 1.1 while all others remain at 1.0. If player 9 has a -5% edge, than the adjustment factor for him is 0.95.

I suppose non-linear factors can be used but the function would just be a guess, so I stuck with the simple one.
This whole algorithm looks kinda similar to the theory of doubling up (to this maths amateur, anyway and assuming 50% win/lose.) That's what gave me the idea for a skill component, especially if the skill adjustments themselves are somehow less sensitive to tournament size than the TODU ones (where generating normal roi for different tournament sizes highlights the limitation of the doublingup chance being constant imo)

ty for the java plexiq and congrats tr sweet discovery
09-24-2011 , 02:25 AM
For implementation purposes, I think using a logarithmic transformation speeds up the process ( in R at least):

Code:
```> x <- runif(1)
> q <- runif(1,0,2)
> time <- proc.time()
> for(i in 1:10000000){
+   x^(1/q)
+ }
> proc.time()-time
user  system elapsed
23.05    0.03   23.14
>
> time <- proc.time()
> for(i in 1:10000000){
+   (1/q)*log(x)
+ }
> proc.time()-time
user  system elapsed
15.45    0.00   15.52```
09-24-2011 , 03:04 AM
Quote:
Originally Posted by Uitje
For implementation purposes, I think using a logarithmic transformation speeds up the process ( in R at least):
Thanks! I hadn't thought of this (been some time since I worked with logarithms). This gives about a 33% speedup i C# too.

m