# bivariate normal distribution probability

We have two genes X and Y. Let $(X,Y)\sim N(\mu_x=9,\mu_y=10,\sigma^2_x=3,\sigma^2_y=5,\rho\sigma_x\sigma_y=2)$. To find $P(X+0.5<Y)$ the probability that the sample mean for the second gene exceeds the sample mean of the first gene by more than 0.5

Which can convert to $P(X+0.5<Y)$ => $P(Y-X>0.5)$

From my understanding, i probably can determine the mean and variance of the distribution Z:Y-X and then use the CDF of a univariate normal distribution?

UPDATE:

Thanks whuber and Glen_b pointed out my mistakes.

Since (X,Y) is normally distributed then $Y-X\sim N(\mu_y-\mu_x, \sigma_y^2-2\rho \sigma_x\sigma_y+\sigma_x^2)\sim N(1, 25-4+9)$ => $\sim N(1, 30)$

Using R we get,

 > 1- pnorm(0.5,mean=1,sd=sqrt(30))
[1] 0.5363678


Additionally, how do i use Monte Carlo simulation to approximate this probability, providing a 95% confidence interval for your estimation?

• @Xi'an oh sorry! just want to know if i'm on the right track solving this problem. – ads27 Feb 5 '15 at 21:58
• @Xi'an also, i'm wondering how to use Monte Carlo simulation to approximate this probability. – ads27 Feb 5 '15 at 22:01
• If you generate $n$ pairs $(X,Y)$ from the bivariate normal you can count how many times $Y-X$ is larger than $0.5$. – Xi'an Feb 5 '15 at 22:04
• Normally in this context "$\rho$" is a correlation, which must lie between $-1$ and $1$. Obviously it is not, because you claim it is $2$. What exactly does it mean then? You also seem to assume $\sigma_y^2 = 25$ and $\sigma_x^2=9$ in your calculation, but these values differ from the ones you specified originally. – whuber Feb 5 '15 at 22:08
• How did you decide what the covariance between X and Y was? Or was this information that was given but is (still) not included in your question? – Glen_b Feb 6 '15 at 0:18

## 1 Answer

For the Monte Carlo simulation,

 # random 50 sample data from multinomial
> gene <- rmvnorm(50,mean=c(9,10),sigma=matrix(c(3,2,2,5),nrow=2))
# sample means in each of the 50 runs
> mean.gene <- apply(gene,1,mean)
# sample variances in each of the 50 runs
> var.gene<- apply(gene,1,var)
> mean(mean.gene) + c(-1,1)*1.96*sqrt(var(mean.gene)/1000) #95%CI for mean
> mean(var.gene) + c(-1,1)*1.96*sqrt(var(var.gene)/1000) #95%CI for variance