probability of a difference between two sampling means of two populations If I have two normally distributed populations $\mu_1 = 7, \sigma_1 = 0.5$ and $\mu_2 = 6.6, \sigma_2 = 0.5$ and I sample each of the populations say with 10 samples each. If I want to work out the probability that $(\overline X_1 - \overline X_2) > 0.6$, would I use the following formula? 
$$Z = \frac{(\overline X_1 - \overline X_2)-(\mu_1 - \mu_2)}{\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}} = \frac{0.6-0.4}{\sqrt{\frac{0.5^2}{10}+\frac{0.5^2}{10}}}= \frac{0.2}{0.2236}=0.8945$$
Therefore, 
$$P((\overline X_1 - \overline X_2)>0.6)=P(Z>0.8945) = 1 - P(Z<0.8945) = 0.1855$$ 
That's my theoretical answer. However, I believe the distribution of (\overline X_1 - \overline X_2) should be normal with $\mu_{\overline X_1 - \overline X_2} = \mu_1-\mu_2 = 0.4$ and the standard deviation of this distribution should be: 
$$\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}=\sqrt{\frac{0.5^2}{10}+\frac{0.5^2}{10}} = 0.2236$$
My question is, that if I run the code to build the sampling distributions of the samples, and replicate the sampling for each population 1000 times, and then run histograms, as below. Can I establish this probability from the histograms instead of the theoretical way above? i.e. the probability that the difference between the sampling means $\overline X_1 - \overline X_2 > 0.6$
 simulate.diff.norm = function(n, R, mu, sigma){
      x_bar = rep(NA, R)
      for(i in 1:R){
        data = rnorm(n, mean=mu, sd=sigma)
        x_bar[i] = mean(data)
      }
      return(x_bar)
    }

SampleDist1 = simulate.diff.norm(n=9, R=10000, mu=7, sigma=0.5
SampleDist2 = simulate.diff.norm(n=9, R=10000, mu=6.6, sigma=0.5
hist(SampleDist1, main = paste("Histogram of Pop A Sampling Mean"))
plot.new()
hist(SampleDist2, main = paste("Histogram of Pop B Sampling Mean"))

 A: I ended up changing that code to simulate 10000 experiments of the two sample distributions together and built another distribution consisting of the difference between the means of each of the nine-sample experiments, for each of the 1000 loops. This resulted in three vectors to feed into the histogram, the sampling distribution for population A, the sampling distribution for population B, and the sampling distribution for the difference between the means of each ith sample of the populations i:10000  
Then I could compare the theoretical result with the result of the simulation
simulate.diff.norm = function(n, R, mu1, mu2, sigma1, sigma2){
  SchoolA_x_bar = rep(NA, R)
  SchoolB_x_bar = rep(NA, R)
  X_bar1_minus_X_bar2 = rep(NA, R)
  for(i in 1:R){
    SchoolA = rnorm(n, mean=mu1, sd=sigma1)
    SchoolB = rnorm(n, mean=mu2, sd=sigma2)
    SchoolA_x_bar[i] = mean(SchoolA)
    SchoolB_x_bar[i] = mean(SchoolB)
    X_bar1_minus_X_bar2[i] = SchoolA_x_bar[i] - SchoolB_x_bar[i]
  }
  return(list(SA=SchoolA_x_bar, SB=SchoolB_x_bar, X1X2=X_bar1_minus_X_bar2))
}

# Theoretical model
mydata = rnorm(10000, mean=0.4, sd=0.2236)
theoretical.diff = sum(mydata>0.6)/10000
theoretical.diff
#start the simulation
Simulated_SA_Minus_SB_Sampling_Mean = simulate.diff.norm(n=10, R=10000, mu1=7, mu2=6.6, sigma1=0.5, sigma2=0.5)
hist(mydata, main = bquote(paste("Theoretical ",E(bar(XA)-bar(XB)))), xlab="Mean difference between sample School A and sample School B")
## now plot the histograms of the simulation
plot.new()
hist(Simulated_SA_Minus_SB_Sampling_Mean[["SA"]], main = "Simulated Sampling Mean for School A")
hist(Simulated_SA_Minus_SB_Sampling_Mean[["SB"]], main = "Simulated Sampling Mean for School B")
hist(Simulated_SA_Minus_SB_Sampling_Mean[["X1X2"]], main = bquote(paste("Simulated Mean Difference ",E(bar(XA)-bar(XB)))))
simulated.diff = sum(Simulated_SA_Minus_SB_Sampling_Mean[["X1X2"]]>0.6)/10000
simulated.diff    

with the following output: 
> # using a sample size of 9 run 10000 simulations of sampling of each of the two populations.
> # find the mean of each sample and find the difference between these two means for each i:10000 sampling experiments
> simulate.diff.norm = function(n, R, mu1, mu2, sigma1, sigma2){
+   SchoolA_x_bar = rep(NA, R)
+   SchoolB_x_bar = rep(NA, R)
+   X_bar1_minus_X_bar2 = rep(NA, R)
+   for(i in 1:R){
+     SchoolA = rnorm(n, mean=mu1, sd=sigma1)
+     SchoolB = rnorm(n, mean=mu2, sd=sigma2)
+     SchoolA_x_bar[i] = mean(SchoolA)
+     SchoolB_x_bar[i] = mean(SchoolB)
+     X_bar1_minus_X_bar2[i] = SchoolA_x_bar[i] - SchoolB_x_bar[i]
+   }
+   return(list(SA=SchoolA_x_bar, SB=SchoolB_x_bar, X1X2=X_bar1_minus_X_bar2))
+ }
> 
> # Theoretical model
> mydata = rnorm(10000, mean=0.4, sd=0.2236)
> theoretical.diff = sum(mydata>0.6)/10000
> theoretical.diff
[1] 0.1866
> #start the simulation
> Simulated_SA_Minus_SB_Sampling_Mean = simulate.diff.norm(n=10, R=10000, mu1=7, mu2=6.6, sigma1=0.5, sigma2=0.5)
> hist(mydata, main = bquote(paste("Theoretical ",E(bar(XA)-bar(XB)))), xlab="Mean difference between sample School A and sample School B")
> ## now plot the histograms of the simulation
> plot.new()
> hist(Simulated_SA_Minus_SB_Sampling_Mean[["SA"]], main = "Simulated Sampling Mean for School A")
> hist(Simulated_SA_Minus_SB_Sampling_Mean[["SB"]], main = "Simulated Sampling Mean for School B")
> hist(Simulated_SA_Minus_SB_Sampling_Mean[["X1X2"]], main = bquote(paste("Simulated Mean Difference ",E(bar(XA)-bar(XB)))))
> simulated.diff = sum(Simulated_SA_Minus_SB_Sampling_Mean[["X1X2"]]>0.6)/10000
> simulated.diff
[1] 0.1882    

Very similar results throughout.
Thanks.
