How to find P(A > B) for two normal distributions in R?

Question

I'm attempting to calculate the probability that normal distribution A is greater than normal distribution B, given the following:

A ~ N(mu1, var1)

B ~ N(mu2, var2)

So far, my attempt has been as follows:

1. Find the intersections of the two normal distributions to get x1 (upper intersect) and x2 (lower intersect)
2. Calculate the probability within the interval of x1 and x2 for distribution A by: areaA <- pnorm(x1, mean(A), sd(A)) - pnorm(x2, mean(A), sd(A))
3. Calculate the probability within the interval of x1 and x2 for distribution B by: areaB <- pnorm(x1, mean(B), sd(B)) - pnorm(x2, mean(B), sd(B))
4. Subtract areaB from areaA to get the probability of A > B

Does this logic make sense? Is there any easier way to do this than the way I've done it?

Thanks!

Answer 1

You already assume $A$ and $B$ are independent, so $C=A-B$ is also normally distributed with mean $\mu=\mu_1-\mu_2$ and variance $\sigma^2=\sigma_1^2+\sigma_2^2$. Then, you can use pnorm to find $P(C>0)$.