On Wed, Aug 20, 2008 at 7:56 AM, Alberto Monteiro <[EMAIL PROTECTED]> wrote: > > 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)
Which can be even more succinctly expressed as mean(delta > 0) It's also often informative to plot the running estimate so you have some idea whether or not you've done enough samples: plot(cumsum(delta > 0) / seq_along(delta), type="l") Hadley -- http://had.co.nz/ ______________________________________________ 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.
