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?


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?

  • $\begingroup$ @Xi'an oh sorry! just want to know if i'm on the right track solving this problem. $\endgroup$
    – ads27
    Feb 5, 2015 at 21:58
  • $\begingroup$ @Xi'an also, i'm wondering how to use Monte Carlo simulation to approximate this probability. $\endgroup$
    – ads27
    Feb 5, 2015 at 22:01
  • $\begingroup$ If you generate $n$ pairs $(X,Y)$ from the bivariate normal you can count how many times $Y-X$ is larger than $0.5$. $\endgroup$
    – Xi'an
    Feb 5, 2015 at 22:04
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber
    Feb 5, 2015 at 22:08
  • 1
    $\begingroup$ 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? $\endgroup$
    – Glen_b
    Feb 6, 2015 at 0:18

1 Answer 1


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.