Simple estimate of a probability by simulation
On Wed, Aug 20, 2008 at 7:56 AM, Alberto Monteiro
<albmont at centroin.com.br> 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