Example 8.36: Quadratic equation with real roots

April 29, 2011
By

(This article was first published on SAS and R, and kindly contributed to R-bloggers)

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.

To leave a comment for the author, please follow the link and comment on his blog: SAS and R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags: , , , , , , ,

Comments are closed.