[R] Simple estimate of a probability by simulation

[davidr at rhotrading.com]

I get 5/36 + log(2)/6, not 1/9.

David L. Reiner, PhD
Head Quant
Rho Trading Securities, LLC

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of Alberto Monteiro
Sent: Wednesday, August 20, 2008 7:56 AM
To: Van Wyk, Jaap; r-help at 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 at 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.