# Weighted lottery

A charter school assigns seats by conducting a lottery. In the past, each child was given a number, and balls were selected randomly. This year, they are starting a new policy of giving low-income students 3 balls, and other kids 2 balls.

They have 10 seats to assign, and 31 applications. I don't know how many students are low-income, but let's assume 12.

How is the number of low-income students likely to be changed by the new policy? Is there a way to create a formula where I could plug in the number of low-income students and see the probability distribution?

• Is this a question from a course or textbook? If so, please add the [self-study] tag & read its wiki. Commented Apr 29, 2016 at 19:11
• How many students will be allowed to enrol ? Commented Apr 29, 2016 at 22:14

Let each low-income child (type $P$) receive $p$ balls and each other child (type $Q$) receive $q$ balls. When one of a child's balls is drawn, all the others are then effectively removed (because that child cannot receive more than one place). Let there be balls for $n$ type $P$ children in the lottery (contributing $p\, n$ balls) and balls for $m$ type $Q$ children (contributing $q\, m$ balls).

Let's keep track of the lottery as it is conducted. Let $(i,j)$ denote the state after $i+j$ balls are drawn and $i$ were of type $P$ and $j$ of type $Q$. There will consequently be $p(n-i)$ balls of type $P$ and $q(m-j)$ balls of type $Q$ available.

The initial state is $(0,0)$. The chance of arriving at state $(i,j)$ is the sum of chances of arriving at states $(i-1,j)$ and $(i,j-1)$, each multiplied by the chances of drawing a type $P$ ball (in the first case) or a type $Q$ ball (in the second). If we write those chances as $f(i,j)$, then by counting the balls of each type that are left and assuming each is equally likely to be drawn, we obtain

$$f(i,j) = f(i-1,j) \frac{(n-i+1) p}{(n-i+1)p + (m-j)q} + f(i,j-1) \frac{(m-j+1)q}{(n-i)p + (m-j+1)q}.$$

This is a dynamic program for $f$. It finds the full distribution of results for $k$ total draws with just $k(k+1)/2 - 1$ total calculations. This makes it feasible to find exact answers for moderately large lotteries (up to a few thousand draws).

When $p=3$, $q=2$, $n=12$, $m=19$, and $k=10$, the probabilities are rational fractions with up to 21 digits in the denominator, so let's use floating point arithmetic. Up to three digits accuracy, the probability distribution for the number of type $P$ children is

      Count:     0     1     2     3     4     5     6     7     8     9    10
Probability: 0.000 0.005 0.037 0.135 0.265 0.294 0.186 0.065 0.012 0.001 0.000


The expectation is obtained by multiplying the counts by the probabilities and adding those all up; it equals $\frac{271557487575886031239}{57968574550829433400}\approx 4.684564$.

We may compute the probabilities for all possible states and plot them in a bubble plot. I have drawn lines corresponding to total draws $k=1, 2, 12+19$ in graduated colors. The total area of bubbles crossing any given line is unity. for example, the preceding figure plots the bubble areas (as bar heights) for the line $k=10$ proceeding from $(0,10), (1,9), \ldots, (10,0)$. (Notice that, as with matrix notation, the first index of the state gives the row $i$ while the second gives the column $j$.)

Note how the center of the distribution rises from $(0,0)$ at roughly a unit slope before tapering off. The $3:2$ weights thereby cause approximately equal numbers of each type to be chosen until around $10-11$ children are drawn, after which the ratios (necessarily) have to return to the $12:19$ ratio within the population as more and more seats are filled.

It may be of interest to observe that this analysis applies equally well to the case where each child has a single ball available and the relative chances of drawing the balls are $p$ and $q$: those values do not have to be integers. It should also be apparent that the same solution method applies when there are more than two types of balls (although the complexity of the solution scales badly, growing from quadratic to cubic to quartic, etc.)

The R code that produced these results was implemented exactly as given in the formulas above. It uses an auxiliary array X (indexed, according to R conventions, starting at $1$ instead of $0$) to store the chances as they are computed.

prob <- function(np=3, nq=2, n=12, m=19) {
X <- matrix(NA, nrow=n+1, ncol=m+1)
X[1,1] <- 1
f <- function(i, j) {
if (is.na(X[i+1, j+1])) {
x <- y <- 0
if (i > 0) x <- f(i-1, j) * (n-i+1)*np / ((n-i+1)*np + (m-j)*nq)
if (j > 0) y <- f(i, j-1) * (m-j+1)*nq / ((n-i)*np + (m-j+1)*nq)
X[i+1, j+1] <<- x+y
}
return(X[i+1, j+1])
}
f(n, m)
return(X)
}
X <- prob()


To extract the distribution for $10$ draws I calculated

p <- X[cbind(1:11, 11-(0:10))]
names(p) <- 0:10
(round(p, 3))


The barplot of that distribution was carried out with the formula barplot(p, main="Distribution").

### Edit

In case you are concerned about doubts raised in other posts in this thread, compare these theoretical results to a (reproducible) simulation in R:

set.seed(17)
n <- 1e6
fill <- function(i,j,n,m) c(rep(i,n), rep(j,m))
sim <- table(replicate(n, sum(sample(fill(1,0,12,19), 10, prob=fill(3,2,12,19)))))/n
plot(sim)


This one takes ten seconds to simulate a million instances of the lottery. In the graphic below, the estimated probabilities are plotted invisibly. They are surrounded by red vertical bars extending two standard errors above and below the estimates. (They may be difficult to see due to their shortness!) Light gray bars extend from the baseline up to the bottoms of the red bars. The values as computed using the preceding theory are marked with horizontal black lines. Thus, we hope that these horizontal lines transect the red sections of the bars and we will look closely at any discrepancies as evidence of a possible error in the theory.

A $\chi^2$ test of the fit gives a statistic of $6.2468$. With $11-1=10$ degrees of freedom, the p-value is $0.7941$, which is so large that it supports our visual perception that the theoretical answers agree with the simulation.

For those who would like to experiment further, here is the post-processing code I used. It includes the exact solution (in the array y), runs the $\chi^2$ test, computes the standard errors, and produces the plot.

y <- c(247/931364,
1887782955567/366743378307796,
491014359017/13097977796707,
187298940516675/1388385646450942,
4779993416332215/18049013403862246,
5314336834283121/18049013403862246,
41911565382746888/225612667548278075,
5715267924185665329/88440165678925005400,
151370444683499943/13073763622101957320,
226164871853/243347902258396,
315675954/13225429470565)
names(y) <- 0:10

x <- rep(0, 11); names(x) <- 0:10
x[names(sim)] <- sim*n
chisq.test(x, p=y, simulate.p.value=(min(y)*n < 5), B=min(1e4, floor(1e7 / n)))

se <- sqrt(sim * (1-sim) / n)
mu <- sum(as.numeric(names(sim)) * sim)
plot(0:10, y + 2*max(se), type="n", bty="n",
main="Simulation vs. Theory",
ylab="Probability",
xlab=paste("Mean =", format(mu, digits=5)))
lines(rowSums(expand.grid(c(-1,1,NA)/2.5, 0:10)), rep(y, each=3))
lines(sim + 2*se, col="#a02020", lwd=5, lend=1)
lines(sim - 2*se, col="#f0f0f0", lwd=5, lend=1)


## As whuber remarks below, this is only an approximate answer

Let

• $n_{low}$ be Number of kids applying
• $n_{high}$ be Number of low income kids applying
• $n_{all}$ be Number of high income kids applying
• $b_{all}$ be Number of kids applying
• $P_{low}$ be Probability of an individual low income kids receiving a seat
• $E_{low}$ be Expected number of low income kids receiving a seat

Then

• $b_{all} = 3 . n_{low} + 2 . n_{high} = n_{low} + 2 . n_{all}$
• $P_{low} = 3 / b_{all}$
• $E_{low} = n_{low}.P_{low} = n_{low} . 3 / b_{all} = \frac{3 . n_{low}}{n_{low} + 2 . n_{all}}$

So if $n_{low} = 12$, then $E_{low} = \frac{3 . 12}{12 + 2 . 31} = \frac{36}{74} = 48.65\%$

• You appear to claim the expected number of low income kids is $36/74$, but that's obviously far too little. After all, it must exceed the expected number when all kids get equal chances, which is $12\times(10/31)\approx 3.9$. Since originally there are $3\times 12=36$ balls for low-income kids and $2\times 19=38$ balls for high-income kids, the expectation ought to be close to $10\times 36/(36+38)\approx 4.86$, but it should be a tiny bit less due to the fact that three balls are removed from the urn every time a low-income kid is chosen but only two balls otherwise. It's close to $4.86486$.
– whuber
Commented Apr 29, 2016 at 20:29
• Good remark. Maybe we can't find a closed formula answer and need a summulation. Commented Apr 29, 2016 at 21:03
• It's a worse than that in terms of obtaining formulas, but for getting answers a dynamic program is a good way to go. Regarding your edit, I am not trying to say your answer is approximate: as stated, it's completely wrong. Possibly you're trying to propose $10\times (36/74)$ as the approximation, but if that's your intention you should state it and work out (a) under what circumstances it's a reasonable approximation and (b) how much it is likely to err.
– whuber
Commented Apr 29, 2016 at 21:06

This is a SAS version of the answer of whuber. I created it because my whubers answer contradicted my simulation. As this program gives the same results as whuber's I now suspect my simmulation.

To be proceeded.

%macro enrole (n, m, p, q);
data take0;
&setLabels;
* you have a 100% chance of starting with no students of both income groups *;
i = 0;
j = 0;
f = 1;
run;

%do draw = 1 %to &n + &m;
data take&draw.(drop=old_: lag_: f_p f_q);
label
i = 'low income studens enrolled'
j = 'high income studens enrolled'
f = 'probability this happens'
;
set take%eval(&draw.-1) (rename=(i=old_i j=old_j f=old_f)) end=last;

f_p = old_f * (&n. - old_i) * &p. / ((&n. - old_i) * &p. + (&m. - old_j) * &q.);
f_q = old_f * (&m. - old_j) * &q. / ((&n. - old_i) * &p. + (&m. - old_j) * &q.);
label
f_p = 'probability of electing a low income student'
f_q = 'probability of electing a high income student'
;
if _N_ = 1 and f_p > 0 then do;
i = old_i + 1;
j = old_j;
f = f_p;
output;
end;
else do;
f = f_p + lag_f_q;
if f > 0 then do;
i = old_i + 1;
j = old_j;
output;
end;
end;
if last and f_q > 0 then do;
i = old_i;
j = old_j + 1;
f = f_q;
output;
end;

lag_f_q = f_q;
retain lag_:;
run;
%end;

data take_all;
set
%do draw = 1 %to &n + &m;
take&draw.
%end;
;
class = i + j;
run;

proc print data= take_all;
run;

proc means data= take_all min p25 mean p75 max;
by class;
var i;
weight f;
run;
%mend;

%enrole (n = 12, m = 19, p = 3, q = 2);


## Computer simulation in SAS

### Attention, I suspect my answer is wrong, because the average deviates from the theoretical one calculated by wubher and verified by me.

1000 trials to fill a class of 10 students from 31 applicants from which 12 are low income gave me an average of 4.067 low income students in the class. (minimum 1, p25 3, p75 5, maximum 8

The code is

%macro simulate(nAll,nLow,nClass,nTrials);
data scores;
format lowEnrolled 8.;
stop;
run;

%do trial = 1 %to &nTrials;
proc datasets noprint;
delete class;
run;

data bucket;
do applicant = 1 to &nAll;
if (applicant le &nLow) then income = 'low ';
else income = 'high';

output;
output;
if applicant le &nLow then output;
end;
run;

%do student = 1 %to &nClass;
data enrollment;
set bucket nobs=nBalls;
retain take;
if _N_ eq 1 then take = ceil(rand('uniform', 0, nBalls));
if _N_ eq take then do;
call symput('enroll', applicant);
output;
end;
run;

proc append data=enrollment base=class;
run;

proc sql;
delete * from bucket where applicant = &enroll;
quit;
%end;

proc sql;
insert into scores
select count(*) as lowEnrolled from class where income = 'low ';
quit;
%end;

proc sgplot data=scores;
histogram lowEnrolled;
run;

proc means data=scores min p25 mean median p75 max;
var lowEnrolled;
run;
%mend;

%simulate(31,12,10,100);

• As my results differ from those of whuber, one of us must have made an error. Commented Apr 30, 2016 at 3:55