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 '15 at 21:58
  • $\begingroup$ @Xi'an also, i'm wondering how to use Monte Carlo simulation to approximate this probability. $\endgroup$ – ads27 Feb 5 '15 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 '15 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 '15 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 '15 at 0:18

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, privacy policy and cookie policy

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