[R] Simple estimate of a probability by simulation

[Alberto Monteiro]

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