[R] Simple estimate of a probability by simulation

hadley wickham h.wickham at gmail.com
Wed Aug 20 15:57:17 CEST 2008


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

-- 
http://had.co.nz/



More information about the R-help mailing list