SAS and R

We often simulate data in SAS or R to confirm analytical results. For example, consider the following problem from the excellent text by Rice:

Let U1, U2, and U3 be independent random variables uniform on [0, 1]. What is the probability that the roots of the quadratic U1*x^2 + U2*x + U3 are real?

Recall that for a quadratic equation A*x^2 + B*x + C to have real roots we need the discriminant (B^2-4AC) to be non-negative.

The answer given in the second and third editions of Rice is 1/9. Here's how you might get there:

Since, B^2> 4*A*C <=> B> 2*sqrt(A*C), we need to integrate B over the range 2*sqrt(a*c) to 1, then integrate over all possible values for A and C (each from 0 to 1).

Another answer can be found by taking y = b^2 and w = 4ac and integrating over their joint distribution (they're independent, of course). That leads to an answer of approximately 0.254. Here's how to calculate this in R:

f = function(x) {
 A = x[1]; B = x[2]; C = x[3];
 return(as.numeric(B^2> 4*A*C))
}
library(cubature)
adaptIntegrate(f, c(0,0,0), c(1,1,1), tol=0.0001, max=1000000)

which generates the following output:

$integral
[1] 0.2543692

$error
[1] 0.005612558

$functionEvaluations
[1] 999999

$returnCode
[1] -1

We leave the details of the calculations aside for now, but both seem equally plausible, at first glance. A quick simulation can suggest which is correct.

For those who want more details, here's a more complete summary of this problem and solution.

SAS

Neither the SAS nor the R code is especially challenging.

data test;
 do trial = 1 to 10000;
 u1 = uniform(0); u2 = uniform(0); u3 = uniform(0);
 res = u2**2 - 4*u1*u3;
 realroot = (res ge 0);
 output;
 end;
run;

proc print data=test (obs=10); 
run;

proc means data=test;
 var realroot;
run;

Leading to the result:

 The MEANS Procedure

 Analysis Variable : realroot

 N Mean Std Dev Minimum Maximum
 -----------------------------------------------------------------
 10000 0.2556000 0.4362197 0 1.0000000
 -----------------------------------------------------------------

R

numsim = 10000
u1 = runif(numsim); u2 = runif(numsim); u3 = runif(numsim)
res = u2^2 - 4*u1*u3
realroots = res>=0
table(realroots)/numsim

With the result

realroots
FALSE TRUE 
0.747 0.253 

The simulation demonstrates both that the first solution is incorrect. Here the simulation serves as a valuable check for complicated analysis.

Insights into where the 1/9 solution fails would be welcomed in the comments.

Smith College is a residential women's liberal arts college in Northampton, MA that is steeped in tradition. One such tradition is to give each student at graduation a diploma at random (or more accurately, in a haphazard fashion). At the end of the ceremony, a diploma circle is formed, and students pass the diplomas that they receive to the person next to them, and step out once they've received their diploma.

This problem, also known as the "hat-check" problem, is featured in Fred Mosteller's "Fifty Challenging Problems in Probability (with solutions)". Variants provide great fodder for probability courses.

The analytic (closed-form) solution for the expected number of students who receive their diplomas when they first receive one is very straightforward. Let X_i be the event that the ith student receives their diploma. E[X_i] = 1/n for all i, since the diplomas are uniformly distributed. If T is defined as the sum of all of the events X_1 through X_n, E[T] = n * 1/n = 1 by the rules of expectations. It is sometimes surprising to students that this result does not depend on n. The variance is trickier, since the outcomes are not independent (if the first student receives their diploma, the probability that the others will increases ever so slightly).

For students, the use of empirical (simulation-based) problem solving is increasingly common as a means to complement and enhance analytic (closed-form) solutions. Here we illustrate how to simulate the expected number of students who receive their diploma as well as the standard deviation of that quantity. We assume that n=650.

R
In R, we set up some constants and a vector for results that we will use. The list of students vector can be generated once, with the permuted vector of diplomas generated inside the loop. The == operator (section B.4.2) is used to compare each of the elements of the vectors.

numsim = 100000
n = 650 
res = numeric(numsim)
students = 1:n
for (i in 1:numsim) { 
 diploma = sample(students, n)
 res[i] = sum(students==diploma)
}

This generates the output:

> table(res)
 0 1 2 3 4 5 6 7 8 
36568 36866 18545 6085 1590 295 40 9 2 
> mean(res)
[1] 1.00365
> sd(res)
[1] 0.9995232

The expected value and standard deviation of the number of students who receive their diplomas on the first try are both 1.


SAS

In SAS we prepare by making a data set with all 650 * 100,000 diplomas by looping through the 100,000 simulations, creating 650 diplomas each time. We also assign a random number to each diploma. Then we sort on the random value within each simulation.

data dips;
do sim = 1 to 100000;
 do diploma = 1 to 650;
 randorder = uniform(0);
 output;
 end;
 end;
run;

proc sort data=dips; 
 by sim randorder; 
run;

Next, we create a student identifier based on position in the data set and mark whether the reordered diploma matches the student ID. Finally, we count how many times the diploma matches the student using proc summary (section 2.1). It would also be relatively easy to count this manually within the data step.

data match;
set dips;
 student = mod(_n_, 650) + 1;
 match = (student eq diploma);
run;

proc summary data=match;
 by sim;
 var match;
 output out=summatch sum=totalmatches;
run;

We can generate a table of the results using proc freq (section 2.3), and the mean and standard deviation using proc means.

proc freq data=summatch;
 tables totalmatches / nocum;
run;

totalmatches Frequency Percent
 0 36726 36.73
 1 36734 36.73
 2 18359 18.36
 3 6241 6.24
 4 1578 1.58
 5 301 0.30
 6 57 0.06
 7 4 0.00

proc means data=summatch mean std;
 var totalmatches;
run;

Analysis Variable : totalmatches
 Mean Std Dev
 1.0036200 1.0031734