I get 5/36 + log(2)/6, not 1/9. David L. Reiner, PhD Head Quant Rho Trading Securities, LLC
-----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Alberto Monteiro Sent: Wednesday, August 20, 2008 7:56 AM To: Van Wyk, Jaap; r-help@r-project.org Subject: Re: [R] Simple estimate of a probability by simulation Jaap Van Wyk wrote: > > I would appreciate any help with the following. > Problem: Suppose A, B and C are independent U(0,1) random variables. > What is the probability that A(x^2) + Bx + C has real roots? I have > done the theoretical work and obtained an answer of 1/9 = 0.1111. > Now I want to show my students to get this in R with simulation. > Below are two attemps, both giving the answer to be about 0.26. > > Could anybody please help me with providing a more elegant way to do > this? (I am still learning R and trying to get my students to learn > it as well. I know there must be a better way to get this.) I must > be doing something wrong ? > Always think that R is vector-oriented, so you should think in terms of vectors. A simple solution for your problem would be: n <- 10000 # n <- 10000? Make n <- 1000000 ! n <- 1000000 a <- runif(n) # a vector of n unifs b <- runif(n) c <- runif(n) delta <- b^2 - 4 * a * c # a vector of deltas # The answer then is how many deltas are non-negative: sum(delta > 0) / length(delta) # Explanation b^2 is the vector of the squares of b's a * c does the pointwise product of vectors delta > 0 is a logical array with TRUE or FALSE sum(delta > 0) coerces TRUE to 1 and FALSE to 0 length(delta) is the length of delta (n) Alberto Monteiro ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
