Determining when the first point in a simulation will exceed a certain value

Question

I have the following question about determining when the first point in a simulation will exceed a certain value:

Suppose you have two random variables "a" and "b" - let's say that both of these random variables have a Normal Distribution, N ~ (10,1). If I were to simulate points from "a" and "b", I could expect the following :

• When simulating points belonging to "a" and "b", each new simulation has absolutely no dependence on previous simulations.

• However, points simulated from the Normal Distribution are more likely to belonging to the "central regions" of the distribution (i.e. closer to the mean).

Question: Suppose I take my two variables "a" and "b" - on average, how many simulations would be required before "a" and "b" are both greater than 12?

Using the R programming language, I tried to write a simulation to answer this question. The simulation below goes as follows:

Step 1: Keep generating two random numbers "a" and "b" until both "a" and "b" are greater than 12

Step 2: Track how many random numbers had to be generated until it took for Step 1 to be completed

Step 3: Repeat Step 1 and Step 2 100 times

Here is the R code:

res <- matrix(0, nrow = 0, ncol = 3)    

for (j in 1:100){
    a <- rnorm(1, 10, 1)
    b <- rnorm(1, 10, 1)
    i <- 1
    while(a < 12 | b < 12) {
        a <- rnorm(1, 10, 1)
        b <- rnorm(1, 10, 1)
        i <- i + 1
    }
    x <- c(a,b,i)
    res <- rbind(res, x)
}

res = data.frame(res)

We can plot a histogram that details the results of the simulation (I repeated it 10,000 times):

hist(res$X3, main = "Results of Simulation", breaks = 100)

And since this is an "irregular looking distribution", there might exist some "better" statistics to summarize the results of this simulation:

#Define Mode Function
Mode <- function(x) {
    ux <- unique(x)
    ux[which.max(tabulate(match(x, ux)))]
}

> mean(res$X3)
[1] 1891.492

> median(res$X3)
[1] 1296.5

> Mode(res$X3)
[1] 929

For instance, based on the results of the simulations - on average it takes 1,891 iterations for both "a" and "b" to greater than 12.

But is there some mathematical formula that could have been used to analytically calculate these same statistics without running the simulations? Perhaps doing this would allow for confidence intervals on these statistics?

Thanks!

Note: Although "a" and "b" are univariate and independent, they could be considered to have a bivariate normal distribution with the variance-covariance matrix being an identity matrix. I think that in principle, the above simulation would apply to a bivariate normal distribution with a non-identity matrix (e.g. use the 'mvtnorm' library in R to simulate points).

Answer 1

The chance of A being $>2\sigma$ from the mean is $2.28\%$.

pnorm(12, 10, 1, lower.tail = FALSE)  

Thus the chance of A & B (2 independent events both occurring) is equal to $A*B$ or $0.052\%$ of the time.
Therefore expected frequency is once every 1932 times (agrees with the simulation results reasonably well).

The medium and mode is a bit more difficult, but I suspect that the distribution you modeled is a normal distribution tail and thus you could calculate the median and mode accordingly.