Is the lottery a smart bet?

It is said that the lottery is a tax on people who are bad at math. What is intended by this is that the odds of winning the lottery as so miniscule that a rational cost-benefit analysis would recommend against participation. But is this always true? Since the jackpot “rolls over” and accumulates for successive drawings, is there some threshold at which the value of a lottery ticket exceeds the price?

Let’s consider the Megamillions lottery played Tuesday night. As has been well documented elsewhere, the probability of buying a winning number for a particular draw can be calculated thus:

Lottery Math

Number of balls in main cage:	75
Number of balls drawn from main cage:	5
Number of balls in Megaball cage:	15

Number of possible combinations of the six balls:

Probability of drawing correct number:

 

It follows that if the net present value of the lottery jackpot exceeds $N, then a $1 lottery ticket becomes a good investment. Correct?

First, consider that the net present value of the lottery jackpot is not the advertised jackpot. The one-time payout for last week’s $636M jackpot is only $340M, whose value is then diminished by federal, state, and local income taxes at (roughly) 45%. So the winner’s after-tax take is only $187M.

Yet not even this is the net present value. We must also consider the probability that, as happened last night, the jackpot is shared between multiple winners. What is that probability?

Number of Tickets sold: T

The probability that at least one winning lottery ticket (i.e., w ≥ 1) would be sold: p(T)

Each successive ticket t purchased adds to the total probability of a winning ticket p(t) the probability P multiplied by the probability that none of the earlier tickets were successful, which is 1 – p(t-1). Thus:

Last week’s drawing had 336M participants. Thus:

,

for P as calculated above.

We included in that that 73% probability the chance that multiple winning tickets were selected. But since multiple tickets means splitting the prize, then the expected value of the prize is the jackpot divided by the expected (e.g., the average over multiple trials) number of winners. Can this be calculated?

In general, the formula for the expected value of random variable W is

,

i.e., the sum, over all possible values of w, of the product of w and the probability of w. In this case, W is the random variable describing the number of winners, w is a specific number of winners, and p(w) is the probability of that w winners.

Given P as calculated above, let’s shift variables and designate p(T,w) as the probability of exactly w winners given T tickets sold.

The probability of one winner when one ticket is sold: p(1,1) = P

The probability of one winner when two tickets are sold:

The probability of one winner when three tickets are sold:

In our application, the number of tickets T is presumably fixed and known. We are also interested in knowing how the probability changes with w, the number of winners. Each ticket is an independent event, so:

The probability of two winners when two tickets are sold: p(2,2) = PP = P²

The probability of two winners when three tickets are sold:

“Wait a second!” the astute among you are saying about now, “that’s starting to look like the binomial theorem!” Indeed it does. In fact, we can generalize that for exactly k winners among T tickets sold, the probability is:

 ,

i.e., the product of the probabilities of all winners and losers and the number of ways we can combine them.

In our calculation of the expected value, it would be nice if we could apply the binomial theorem,

 ,

but we can’t, because we must weight the terms by as written, because we must weight the terms by w:

Unfortunately, I’m not good enough at math enough to figure out a closed-form solution to the summation above without knowing the answer in advance. And the usual way of evaluating , with factorials as show above, won’t work when T is hundreds of millions of tickets sold; no computer generates factorials that large. But the Wikipedia article suggests using iteration and logarithms to calculate for large T, which I have done in Matlab here:

%%
% Parameters
n = 75; % Number of balls in main bin
v = 5;  % Number of balls drawn from main bin
m = 15; % Number of balls in MegaBin
N = m*nchoosek(n,v); % Number of possible lottery combinations
P = 1/N; % Probability of receiving the winning number
%
% Calculate the binomial coefficients with logarithms
%
T = N;   % Arbitrarily set the number of number of tickets sold to 239E6; we vary this number over multiple runs
i = 1:T; % A vector containing all integers between 1 and T
logi = log(i); % Calculate the logarithm of all integers between 1
% and T
logu = 0;
logl = 0;
for k = i,
logu = logu + logi(T-k+1);
logl = logl + logi(k);
logBC(k) = logu – logl; % the natural log of the binomial coefficient, i.e. log[(T k)]
end;
% Calculate the expected value
logP = log(P);
logQ = log(1-P);
% natural log of the product of k wins and the probability of k wins
logwpw = logi + logBC + i*logP + (T-i)*logQ;
EW = sum(exp(logwpw)); % Calculate the antilog and sum
display(EW);

 

Fortunately, on a Xeon X5672 at 3.2GHz, this algorithm can be executed in a few minutes; note, however, that I optimized it for speed by minimizing the contents of the for loop.  The tradeoff is that I have several N-length floating-point vectors that managed to take up 18GB of memory.  Proceed at your own risk!

I ran the program for several values of T, and convinced myself that E[W] = T/N = PT. Now that I know the answer, I’d like a second crack at a closed form solution:

I will now resort to a bit of hand-waving:  given values of T in the hundreds of millions, the value for E[W] is not measurably diminished by conditioned on a given number of winners.  This is analogous to the probability of rolling double-sixes given that you already rolled one six.  (Answer:  1/6, not 1/36.)  So given that you won the jackpot, the number of additional expected winners is still PT.

Applying this to the problem at hand, it follows that the after-tax value of last night’s jackpot must be divided through by 1 (i.e., you) + 336/259 (the other winners), reducing the value of your jackpot share to $81M.

Now, it is theoretically possible that, given a high enough jackpot, that $81M could increase to the $259M necessary to make the expected value of the lottery ticket greater than the dollar you’re paying for it.  But remember that the number of tickets T for any given drawing will probably go up with the increasing jackpot, diminishing its value to an individual player.  And even if it didn’t, it wouldn’t follow that buying multiple tickets would be a good idea.  I will leave working the math out to the reader, but you would be assured of buying duplicate tickets, and thereby playing against yourself, after a few tens of thousands of tickets purchased, and the model would need to be extended to account for this.

(NOTE:  The foregoing was for entertainment purposes only.  Those wishing to invest their life savings in lottery tickets should consult a qualified dissipation professional.)

Radiometric Error Evaluation, Part V: Results

In the last post, I promised to evaluate the integral:

This is only partially true.  I have no idea how to solve an equation like this analytically, although I know the correct answer.  What I can do, however, is solve it numerically, setting our constant radiance = 1000, the radius of our emitting sphere = 1:

function [R,E]=intE()
    L = 1000;
    re = 1; 
    l1 = 0; % Sensor nadir
    figure;
    for i = 1:200,
        R = 10.^(-1+.02*i);
        l2 = acos(re./(R + re)); % Limit of field of view;
        E(i) = quad(@irradiance,l1,l2);
    end;
    function E = irradiance(phi)
        Rprime = sqrt((R + re)^2 + re^2 - 2*(R+re)*(re).*cos(phi));
        gamma = asin(((R + re)./Rprime).*sin(phi));
        theta = asin((re./Rprime).*sin(phi));
       
        E = (2*pi*L*re^2./Rprime.^2.*sin(phi).*cos(gamma).*cos(theta));
    end
end

I graphed the value of E vs. R, along with the value of our estimate from two lessons ago:

The estimate appears to be pretty good for R >> r_e.  But before I graph the error, I decided to add to the estimate E = π L r_e²/R² the estimate E = π L r_e²/(R +r_e)²:

Well, how about that.  It turns out E = π L r_e²/(R +r_e)²: was so good an estimate it completely overwrote the true value across all values of R.

Let’s plot the error.

By the looks of it, the error of E_s2 (taken as an absolute value to graph on a log scale) is pretty much limited to accumulated random rounding errors, although these errors might be getting larger the farther away my receiving aperture is, and the smaller the emitting surface appears.

Basically, this is saying that for purposes of analysis,

Investigating Irradiance Error, Part IV: True E

Let’s imagine a differential emitting area dA_s on a spherical emitting surface of radius r_e at a distance R from a differential receiving surface dA_R.  dA_s is at distance R’ and aspect angle γ from dA_R; dA_R is at aspect angle θ to dA_s.  Also dA_s is at an angle φ from the point of nadir to dA_R.  In summary,

Applying the Law of Sines gives us:

allowing us to solve for both θ and γ in terms of φ:

The Law of Cosines gives us:

These relationships will be important later.

Radiance is the amount of power emitted into each differential solid angle from each differential area perpendicular to that angle.  For a lambertian surface, it is by definition a constant:

In units of steradians, a differential solid angle dΩ measures the ratio between the surface of a sphere subtended by that angle and the square of the sphere’s radius, much in the way an angle in radians measures the ratio of the length of the arc it subtends and the radius of the circle.

We must also bear in mind that the perpendicular area (i.e. the cross-sectional area as it appears from the direction of observation) bears the relationship to the true area shown in the figure below.

Thus:

Substituting into our equation for radiance gives us:

The quantity in which we are interested is irradiance E, the power per unit area entering our differential receiving area, dA_R.  Solving for this quantity yields:

The differential surface area of a sphere, dA_s, is the standard calculus formula:

where dα integrates around the circle of the visible surface to 2π, and the interior body angle ranges from 0 to the point where  γ is a right angle.  A simple application of the definition of a cosine yields:

Substituting in these values gives us our final integral:

We will evaluate this integral in our next lesson.

Investigating Irradiance Error III: Estimated E

In the last two posts, we set out to find the relationship between four radiometric quantities:

• Exitance (M): the power density emitted by a surface, measured in Watts / m². We will reserve B for the per-micron un-integrated quantity.
• Radiance (L): the power density emitted by a surface into each solid angle, measured in Watts / m² / steradian.
• Intensity (I): the power emitted by a point source into each solid angle, measured in Watts / steradian.
• Irradiance (E): the power incident upon a receiving surface, measured in Watts / m².

We discovered that the relationship between exitance and radiance followed the simple formula for a lambertian surface:  M = πL.  We will now look at the relationship between these and  irradiance.  If anything, properly calculating this last is the most critical; as its name implies, radiometry is about using the power incident on our measuring device to infer the power being emitted from the surface of interest.

Let us assume that a spherical object of radius R_obj is sufficiently far away from our sensor that we can treat it as a point source of intensity I.  The total power Φ it emits is equal to the product of its exitance M and surface area A; its intensity I is that power divided by the total surrounding solid angle, 4π.  Thus:  

You will remember that intensity is in units of Watts / steradian.  Thus, if we assume that the sensor is looking directly at the emitting object, we can find the irradiance at the aperture of our receiver by multiplying the intensity by the solid angle Ω subtended by the aperture and dividing by the area A_rcvr of the aperture.  Remembering that the definition of solid angle is the ratio of the area of a sphere subtended by that angle divided by the square of the radius gives us:

where R_path is the distance from the emitting object to our sensor.

Alternatively, we can treat our emitting object as a disk of radius R_obj.  It’s radiance L is the exitance M divided by π as found previously.  The irradiance on a sensor at distance  is the product of the radiance, the solid angle Ω subtended by the sensor’s aperture, and the disk area A_obj = πR_obj², divided by the area of the sensor aperture:

Thus, we can see that whether we approximate our emitting object of radius R_obj as a point source or a disk, we get the same result for the amount of radiation incident on a sensor at distance R_path.  (The textbook irradiance of, for instance, the sun at the top of the earth’s atmosphere is 1480W/m², for instance.)

However, we must point out the obvious:  an emitting sphere is not a disk, and it may not appear as a point source.  For instance, the sun’s dimensions are apparent to an observer on the earth, 93M miles away.

In our next post, we will calculate the exact irradiance on the aperture of a sensor from a spherical emitting object.

(Note:  I apologize for the crappy looking equations.  I have a procedure for cleaning them up, but it’s tedious, and I’m tired.)

Investigating Irradiance Error, Part II: Relationships

In my last lesson, I defined four radiometric quantities:

• Exitance (M): the power density emitted by a surface, measured in Watts / m². We will reserve B for the per-micron un-integrated quantity.
• Radiance (L): the power density emitted by a surface into each solid angle, measured in Watts / m² / steradian.
• Intensity (I): the power emitted by a point source into each solid angle, measured in Watts / steradian.
• Irradiance (E): the power incident upon a receiving surface, measured in Watts / m².

Radiance is the amount of power emitted into each differential solid angle from each differential area perpendicular to that angle. For a lambertian surface, it is by definition a constant:

where Φ is the power in Watts.

In units of steradians, a differential solid angle dΩ measures the ratio between the surface of a sphere subtended by that angle and the square of the sphere’s radius, much in the way an angle in radians measures the ratio of the length of the arc it subtends and the radius of the circle.  Note that the subtended surface area of a sphere is only an approximation of the subtended area of a flat surface; however, the approximation is good enough at small angles, and exact at differential angles, which is what we are intending here.

Did you notice the little perpendicular sign after the differential sending area dA_s?  This is to remind us that the perpendicular area (i.e. the cross-sectional area as it appears from the direction of observation) bears the relationship to the true area shown in the figure below.

This relationship applies to both the receiving area dA_R and the sending area dA_s, as shown here:

Using these definitions and basic trigonometry gives us the following relationships:

We now have enough equations to solve for exitance M in terms of radiance L by integrating over the surface of a hemisphere defined by the differential surface area dA_R:

First, we point out that the hemisphere is always facing the center, hence θ_R =0°, so cosθ_R =1.  Next, we replace dA_R with the area element for integrating over the surface of a sphere in polar coordinates:

dα conveniently integrates out to 2π, and the R²’s cancel.  We’re only integrating over a hemisphere, so θ_s ranges from 0 to π/2.  Thus:

This relationship, M = πL for a lambertian surface, is sufficiently rigorous that Planck’s formula itself is often expressed in terms of radiance rather than exitance.  However, relating these quantities to irradiance and intensity often make use of assumptions which are less than intuitive.  We will examine these assumptions in the next post.

Math Bait II: The Spring-Mass System with Laplace Transforms

In my last post, we solved the differential equation

with traditional methods.  Today, we will examine the same problem using Laplace Transforms.

Let’s briefly review the definition of the Laplace Transform and a few of its properties:

Returning to our equation, we can apply the properties to rewrite it in the s domain:

and solve for X(s):

In order to find the inverse transform of X(s) and solve the problem, we must decompose the denominator of the first term using partial fraction expansion:

Much as we matched the coefficients for the sine and consine terms in the method of undetermined coefficients in the previous problem, we here must match the coefficients for the powers of s, yielding a system of equations and their solution:

Armed with these coefficients, we can now rewrite X(s) in terms of our expanded fractions:

In order to calculate the inverse Laplace transform, it would be easier if our terms matched those of the properties provided above. Rewriting:

It is now a trivial matter to find the final solution:

Resonant Frequency

Finding the result for the resonant frequency, while algebraically intense, is relatively straightforward with Laplace transforms. We begin by replacing ω with √(k/m) to obtain:

Before proceeding, and since we are expecting the same answer we obtained using traditional methods, let's apply the multiply-by-t property above to the sine function:

Now let's rewrite X(s) above so that it's terms match the formula:

Here again, it is a trivial matter to find the solution: