Playing the lottery

[bokehman]

Hola amigos! I created a lottery class for fun but it has a problem. What it does is choose 6 random numbers between 1 and 49. It then plays those numbers every week until they win the jackpot. It then returns the date on which the numbers won. In theory there are 15,890,700 possible combinations but when I run the script it normal hits the jackpot on average after only 500,000 random numbers. What's wrong? Is the problem with the code or with PHP?

<?php

class lottery
{
	function lottery()
	{
		static $seed = true;
		if($seed)
		{
			srand((double)microtime()*1000000 );
			$seed = false;
		}
		for($i = 1; $i < 50; $i++)
		{
			$this->numbers[] = $i;
		}
	}

function gen_numbers($pulls)
{
	for($i = 0; $i < $pulls; $i++)
	{
		$this->pulls[] = $this->pull();
	}
}

function pull()
{
	shuffle($this->numbers);
	$pull = array_slice($this->numbers, 0, 6);
	sort($pull);
	return $pull;
}

function get_weeks()
{
	$this->weeks = 1;
	while(($this->pull = $this->pull()) and (!in_array($this->pull, $this->pulls)))
	{
		$this->weeks++;
	}
	return $this->weeks;
}

function result()
{
	$output = 'After '.($c = count($this->pulls)).' attempt'.(($c == 1)?null:'s').' per week your numbers ('.implode(', ', $this->pull).') came up on ';
	$ts = (($this->weeks % 52.1785)*7*24*60*60) + time();
	$output .= date('F jS ', $ts) . floor(date('Y', $ts) + ($this->weeks / 52.1785));
	return $output;
}
} 

$bad_news = new lottery;
$bad_news->gen_numbers(1);
echo $bad_news->get_weeks().' weeks<br>';
echo $bad_news->result();

?>

[NogDog]

Might be interesting to come up with a custom shuffle function that uses mt_rand(), as shuffle() depends on the same random number generator as rand(). Other than that, nothing jumps out at me as wrong in the code.

[Weedpacket]

One thing that does jump out at me:

srand((double)microtime()*1000000);

I'd just leave this out. The only reason why you'd need to seed manually is if you want to get the same sequence on different runs. Doing the above is not only unnecessary, but damaging (since you're limiting yourself to only one million possible sequences out of billions).

Those C-library rand() functions are pretty cheap jobs, especially Windows; they tend to sacrifice statistical robustness for speed, and are only good enough for occasional use.

This has come up before: http://www.phpbuilder.com/board/showthread.php?t=10328472

As NogDog suggested, use a better RNG. For shuffling an array I wrote (I guess I should warn you that this is from memory):

function mt_shuffle(&$array)
{
	for($i=1, $count=count($array); $i<$count; $i++)
	{
		$j = mt_rand(0,$i);
		$t = $array[$i];
		$array[$i]=$array[$j];
		$array[$j] = $t;
	}
}

I'd also replace the loop that fills $this->numbers with a call to range().

Oh, and there are only 13,983,816 of choosing six numbers out of 49. If you want 15,890,700 you need to have 50 numbers to choose from.

[bokehman]

Weedpacket wrote:

Oh, and there are only 13,983,816 of choosing six numbers out of 49. If you want 15,890,700 you need to have 50 numbers to choose from.

Yes that right. When I first did the calulation I did it on an array that had a duplicated member. Heres the code I used for that calculation (just in case anyone is interested):

<?php

echo '<p>The number of possible lottery combinations of 6 numbers chosen from numbers between 1 and 49 is: ';
echo '<b>'.number_format(CombinationsWithoutRepetition(49, 6)).'</b>';

echo '<p>The of possible 5 card poker hand combinations possible to make from a 52 card deck is: ';
echo '<b>'.number_format(CombinationsWithoutRepetition(52, 5)).'</b>';

function CombinationsWithoutRepetition($total, $quantity)
{
	return((factorial($total))/(factorial($quantity)*factorial($total-$quantity)));
}

function factorial($in){
   if($in<2) return 1;
   for($f=2; $in-1>1; $f*=$in--){}
   return $f;
}

?>

I simplified my lottery code above to the following which uses mt_rand().

<?php

ini_set('max_execution_time', 100);

$i = 1;
$num = mt_rand(1, 13983816);
while($num !== mt_rand(1, 13983816))
{
	$i++;
}
$ts = (($i % 52.1785)*7*24*60*60) + time();
echo 'You won on ' . date('F jS ', $ts) . floor(date('Y', $ts) + ($i / 52.1785));

?>

The thing is neither of these pieces of code are very nice because they both rely on brute force to come up with an answer. I'd rather come up with some sort of probabilty equation that comes up with a similar realistic result (different every time it runs) but that only takes a couple of milliseconds to run because it doesn't rely on brute force (looping the RNG).

[NogDog]

Well, since we know that the odds of winning with any given set of 6 numbers is 1/13,983,816, there's really nothing gained by simulating the drawing of 6 numbers each week. Instead you could just do:

for($weeks = 1; mt_rand(1, 13983816) != 1; $weeks ++);
echo "It took $weeks weeks to win.";

It's still brute force, but it should be faster. My math skills aren't good enough to figure out a single equation for this, partiuclarly since theoretically it's possible to essentially approach an infinite number of tries without ever actually hitting the numbers. (The probability of not hitting within a given number of attempts approaches zero as the attempts approach infinity, but never actually reaches zero.) There's probably an equation that will calculate the expected average number of weeks, but my statistics skills aren't that good.

[bokehman]

NogDog wrote:

It's still brute force, but it should be faster.

To me on the surface looks pretty similar to the while loop posted above but for some reason the for loop is producing numbers that are too big. According to my calculations there is a 50% chance of having won by week 9.65 million but the majority of the time the for loop comes up with bigger numbers than this. Very strange! The only difference with the while loop is I started with a random number for the comparison instead of one.

Anyone got any idea to do this without brute force?

[bokehman]

Just in case anyone is interested I made a look-up table. 101 entries in the table from 0-100%. Obviously 100% would be an infinite number of weeks so I cut it at 100 million weeks which is a probability of having won of about 99.8%. 0% is week 1.

<?php

$lookup = array(0 => 1, 1 => 140543, 2 => 282511, 3 => 425936, 4 => 570848, 5 => 717276, 6 => 865255, 7 => 1014816, 8 => 1165994, 9 => 1318824, 10 => 1473343, 11 => 1629588, 12 => 1787599, 13 => 1947416, 14 => 2109080, 15 => 2272635, 16 => 2438126, 17 => 2605599, 18 => 2775102, 19 => 2946685, 20 => 3120399, 21 => 3296299, 22 => 3474439, 23 => 3654878, 24 => 3837675, 25 => 4022894, 26 => 4210599, 27 => 4400858, 28 => 4593742, 29 => 4789323, 30 => 4987678, 31 => 5188888, 32 => 5393035, 33 => 5600206, 34 => 5810493, 35 => 6023991, 36 => 6240799, 37 => 6461021, 38 => 6684767, 39 => 6912151, 40 => 7143294, 41 => 7378322, 42 => 7617368, 43 => 7860571, 44 => 8108079, 45 => 8360046, 46 => 8616638, 47 => 8878025, 48 => 9144392, 49 => 9415931, 50 => 9692848, 51 => 9975359, 52 => 10263696, 53 => 10558104, 54 => 10858843, 55 => 11166192, 56 => 11480449, 57 => 11801931, 58 => 12130978, 59 => 12467954, 60 => 12813252, 61 => 13167292, 62 => 13530530, 63 => 13903454, 64 => 14286598, 65 => 14680535, 66 => 15085893, 67 => 15503352, 68 => 15933659, 69 => 16377628, 70 => 16836157, 71 => 17310232, 72 => 17800945, 73 => 18309505, 74 => 18837262, 75 => 19385720, 76 => 19956570, 77 => 20551720, 78 => 21173329, 79 => 21823861, 80 => 22506139, 81 => 23223421, 82 => 23979493, 83 => 24778794, 84 => 25626566, 85 => 26529072, 86 => 27493868, 87 => 28530195, 88 => 29649516, 89 => 30866289, 90 => 32199117, 91 => 33672495, 92 => 35319598, 93 => 37186940, 94 => 39342644, 95 => 41892327, 96 => 45012938, 97 => 49036215, 98 => 54707013, 99 => 64402795, 100 => 100000000);

echo(mt_rand($lookup[($rand = mt_rand(0, 99)) - 1], $lookup[$rand+1])).' weeks';

?> 

[Weedpacket]

Bet on winning this week.

If I had to make a bet on which week I'd win the lottery, the smart money would be not to 🙂. But assuming I haven't got that option, I'd bet on winning it this week. Sounds counterintuitive, but it's the week on which I am most likely to win. This week I have a 1/13983816 chance of winning, and a 13983815/13983816 chance of losing. Next week I have a 1/13983816 chance of winning, but only if I lost this week (we're playing until we win once, after all). So the chance of winning next week is only 13983815/13983816 * 1/13983816 - which is less than 1/13983816. So to win on the tenth week not only do my numbers have to come up on the tenth week, but they must not come up on any of the nine previous weeks.

Winning this week is more likely than not winning this week and winning next week. Which in turn is more likely than not winning this week and not wining next week, and winning the week after that.....

Generally speaking the chance of winning on week n (counting this week as 1) is (1-p)^(n-1) * p, where p is of course the probability.

Try this

<?php
$res=array();
for($t=0; $t<1000000;++$t)
{
	$count=1;
	while(mt_rand(1,10)!=1)
		++$count;
	if(isset($res[$count]))
		++$res[$count];
	else
		$res[$count]=1;
}
ksort($res);
print_r($res);

It picks a number randomly between 1 and 10 inclusive, counting how many it has to pick before a 1 comes up, whereupon it tallies the count. It does this a million times, and then prints the final tally sheet.

What are the chances that it will pick 1 on the first try? One in ten. What are the chances that it will not? Nine in ten. What are the chances that it will pick 1 on the second try? Assuming that gets that far, one in ten. But it only gets that far nine out of ten times. So the total probability of picking 1 on the second try (but not the first) is 9/10*1/10 = 9/100 - slightly less than 10/100.

I just ran the above program, and the first few results were

    [1] => 99954    100000
    [2] => 90171     90000
    [3] => 81199     81000
    [4] => 72752     72900
    [5] => 65986     65610

And along side that I've put the relevant probability (9/10)^(n-1)*(1/10) ,

A direct simulation using your odds would be the same, but with 13983816 instead of 10.

All in all, this is a geometric distribution.

Now you don't want to run the experiment explicitly; you want to be able to generate these "waiting times" directly so that they'll have the correct distribution. Righty-ho. Integrating with respect to n, taking the inverse of the resulting function (and writing the result in PHP):

function get($p)
{
	if($p==1)
		return 1;
	while(0==($r=mt_rand())); // Having $r==0 below would be Bad.
	$r /= mt_getrandmax(); // $r is now a positive number (0,1].
	return log($r)/log(1-$p)+1;
}

Clever thing to do would be to make a class and set the probability $p as one of its properties. Should go without saying you'd want to use intval() to truncate the result.

[bokehman]

Weedpacket wrote:

Now you don't want to run the experiment explicitly; you want to be able to generate these "waiting times" directly so that they'll have the correct distribution.

That was the idea of the look-up table. Each step represents a 1 per cent chance of having already won by that point. Notice how the steps get wider a precisely according to your equation. The only entry that is not true is for 100 per cent as this would be infinate. 99.9% likelihood of having won is at 96,677,875 plays. The look-up table is near instateneous but the shortcoming is that there needs to be a look-up table for each probability.

Weedpacket wrote:

function get($p)
{
	if($p==1)
		return 1;
	while(0==($r=mt_rand())); // Having $r==0 below would be Bad.
	$r /= mt_getrandmax(); // $r is now a positive number (0,1].
	return log($r)/log(1-$p)+1;
}

Clever thing to do would be to make a class and set the probability $p as one of its properties. Should go without saying you'd want to use intval() to truncate the result.

I tried running that code but I got a division by zero message. In trying to track it down I echoed the following:

echo log(1-$p)+1;

This outputs:

-1.#IND

What is that? I've never seen anything like it before.

Lastly I think the look-up table was a good start but it's not dynamic. The 100 step array took PHP about 10 minutes to assemble and I made a 1000 step one also which took 1.5 hours. What dynamic solution could replace the look-up table and be dynamic (so the odds could change) and be almost as fast as the look-up tables (< 5 milliseconds)?

[Weedpacket]

-1.#IND

Not a Number; usually written as NaN. Sounds like the probability $p you supplied was greater than 1 (there's no such thing as log(negative number).

I just tried the function a million times with $p = 1/13983816 and it didn't barf once. Took six seconds to run, too. Gave me waiting times between 38 and 214,439,433 weeks. (So... between 9 June 2007 and 17 May 4111813)

bokehman wrote:

What dynamic solution could replace the look-up table and be dynamic (so the odds could change) and be almost as fast as the look-up tables (< 5 milliseconds)?

That's what I wrote.

[bokehman]

Yeah you're right, I was feeding it the number of combinations rather than the probability.

[bokehman]

Hi, thanks for the code. It perfect for my needs. I tried the following code with the function and it returns within 0.25% of perfect (in 4.5 seconds).

$j = 0;
$k = 0;

for($i = 0; $i < 1000000; $i++)
{
	if(get(1/13983816) < 9692848)
	{
		$j++;
	}
	else
	{
		$k++;
	}
}

echo round(($j/$k)*100, 2).'%';

The code to work out the probability is above but I don't know if this could be simplified. I suck at maths. Also in a case where there are 49 numbers, 6 options but for example only 2 are required for a minor win how do I work out the probability?

[Weedpacket]

Oh, for that you want the cumulative distribution function; the probability that your variate takes on a value <= some specified value.

Since the probability distribution function is pdf(n) = p(1-p)^n, the cdf is the integral of that from 0 up to the point you're interested in. Since it's a discrete distribution you can use a sum instead of a full-on integral.

So you're after

          n         k
CDF(n) = SUM p*(1-p)
         k=0

Which is 1-(1-p)^(n+1), or 1-pow(1-$p, $n+1).

Sure enough 1-pow(1/13983816, 9692848+1) = 0.5000002398

[bokehman]

Damn! I wish my maths skills were as good as yours! I think we are at cross purposes or I don't fully follow the relevance of your answer. Earlier I posted this code.

function CombinationsWithoutRepetition($total, $quantity)
{
    return((factorial($total))/(factorial($quantity)*factorial($total-$quantity)));
}

function factorial($in){
   if($in<2) return 1;
   for($f=2; $in-1>1; $f*=$in--);
   return $f;
}

This returns number of combinations without repetition. I don't know if it it more complicated than necessary, but what I wanted is to know is this: Let's take for example the National Lottery, there are 49 numbers and one line is a selection of 6 numbers. If only 2 numbers from your 6 match the you win a tenner. How do I get the odds for winning a tenner or for three of my balls being selected on any particular week?

[Weedpacket]

Well, you asked two questions I see. I answered the other one.

You have a bag containing 6 white balls and 43 black balls. You draw 6. What is the probability that you'll draw 2 white balls and 4 black balls? Well, let's say you draw them in the order wwbbbb. The probability of that is 6/49 for the first ball, times 5/48 for the second (once one white ball is taken out, there are only five left in a bag of forty-eight), times 43/47 for the third (forty-seven balls, of which forty-three are black), and 42/46, 41/45, 40/44 for the fourth, fifth and sixth balls.

Total probability 6/49 5/48 43/47 42/46 41/45 * 40/44 = 8815/998844 or just under 1%.

If you try the same process out for different orderings of the balls you'll find you get the same result; the denominators don't change of course (if you take a ball out of the bag, the number of balls in the bag decreases by one), and the numerators just get rearranged.

While you could work it out from combinatorics (the denominator is 49!/(49-6)!, the numerator is the product of 6!/(6-2)! and (49-6)!/(49-6-(6-2))! - 49 for the original number of balls in the bag, six are drawn, two are white) and that expression probably simplifies a bit (roll out the Pochhammer polynomial!), it's probably easier just to go with
(65)(43424140) / (494847464544)
Note that there are 6 terms in the denominator and 6 in the numerator. You want two white balls, so that's the first two factors in the numerator. And for the remaining terms in the numerator you start from 43 because if there are exactly six white balls in the bag of forty nine, there are exactly forty three that aren't.

// There is a bag containing $N balls, of which exactly $m are white. We want to
// select exactly $n white balls and ($m-$n) non-white balls.
function select($n,$m,$N)
{
	if($n>$m) return 0; // Can't select more white balls than there are.
	$prob=1;
	for($i=0; $i<$m; $i++)
	{
		$den=($N-$i);
		if($i<$n)
			$num = $m-$i;
		else
			$num = ($N-$m+$n)-$i;
		$prob *= $num/$den;
	}
	return $prob;
}

[bokehman]

QUOTE=Weedpacket(43424140) / (4948474645*44)[/QUOTE]Hi again Weedpacket, your clear logic has really made this easy for me to understand. But still a problem: I ran your function select(args) on 3, 4 and 5 ball wins and it returned:

5 balls: 1/325205
4 balls: 1/15485
3 balls: 1/1133

but according to the national lottery the odds are:

Jackpot - Match 6 main numbers : 1 in 13,983,816
Match 5 main numbers plus the bonus number: 1 in 2,330,636
Match 5 main numbers: 1 in 55,492
Match 4 main numbers: 1 in 1,033
Match 3 main numbers: 1 in 57

I tried the following code to get a brute force answer for 3 balls and it returns 57 each time. Where is the discrepancy coming from? Also can you shed any light on how the odds are worked out for the bonus ball?

for($k = $i = 0; $i < 100000; $i++)
{
	$numbers = range(1,49);
	shuffle($numbers);
	if(count(array_intersect(array(1,2,3,4,5,6), array_slice($numbers,0,6))) > 2)$k++;
}
echo round(100000/$k);