Hallo
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 ?
n <- 10000
### Method 1 ###
cnt <- rep(0,n)
for (i in (1:n)) {
a <- runif(1)
b <- runif(1)
c <- runif(1)
cnt[i] <- ifelse(((b^2)>4*a*c),1,0)
}
sum(cnt)/n
### Method 2 ###
one.s <- function(x) {
ret <- ifelse(((x[2]^2) > 4 * x[1] * x[3]),1,0)
ret
}
m <- cbind(runif(n),runif(n),runif(n))
sum(apply(m,1,one.s))/n
######
THANK YOU.
Jacob L van Wyk
Dept of Statistics
University of Johannesburg
P O Box 524
Auckland Park 2006
South Africa
Tel: +27-11-559-3080
Fax: +27-11-559-2832
Cell: +27-82-859-2031
[[alternative HTML version deleted]]
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