1. A variable was observed as 95 from a normal distribution. The mean and SD of the distribution are most likely to be:
par(mfrow=c(4,1), mar=c(3,3,2,3)) 

# Option A: Mean = 100, SD = 1
x = rnorm(10000, 100, 1)
plot(density(x), xlim=c(60,105), main="Mean=100, SD=1")
abline(v=95, col="red", lty=1)  # Vertical red dashed line at x=95

# Option B: Mean = 100, SD = 2
x = rnorm(10000, 100, 2)
plot(density(x), xlim=c(60,105), main="Mean=100, SD=2")
abline(v=95, col="red", lty=1)  # Vertical red dashed line at x=95

# Option C: Mean = 85, SD = 5
x = rnorm(10000, 85, 5)
plot(density(x), xlim=c(60,105), main="Mean=85, SD=5")
abline(v=95, col="red", lty=1)  # Vertical red dashed line at x=95

# Option D: Mean = 85, SD = 10
x = rnorm(10000, 85, 10)
plot(density(x), xlim=c(60,105), main="Mean=85, SD=10")
abline(v=95, col="red", lty=1)  # Vertical red dashed line at x=95

pnorm(-5, lower.tail = TRUE)  # Option A: Mean = 100, SD = 1, P = 0.0000002866 
## [1] 2.866516e-07
pnorm(-2.5, lower.tail = TRUE)  # Option B: Mean = 100, SD = 2, P = 0.0062 
## [1] 0.006209665
pnorm(2, lower.tail = FALSE)  # Option C: Mean = 85, SD = 5, P = 0.0228 
## [1] 0.02275013
pnorm(1, lower.tail = FALSE)  # Option D: Mean = 85, SD = 10, P = 0.1587 
## [1] 0.1586553

Empirical rules for normal distributions

# Set up the plot area with margins for labels
par(mar = c(5, 4, 4, 4))  # Bottom, left, top, right margins (in lines)

# Generate x values for the standard normal distribution (z-scores from -4 to 4)
x <- seq(-4, 4, length = 1000)
y <- dnorm(x, mean = 0, sd = 1)  # Density of standard normal distribution

# Plot the density curve
plot(x, y, type = "l", lwd = 2, col = "blue", 
     xlab = "Z-Score", ylab = "Density",
     main = "Standard Normal Distribution",
     ylim = c(0, 0.42),  # Adjust y-axis to fit labels
     axes = FALSE)  # Turn off default axes for custom ones

# Add custom axes
axis(1, at = c(-3, -2, -1, 0, 1, 2, 3), labels = c("-3", "-2", "-1", "0", "1", "2", "3"))
axis(2, las = 1)  # Vertical y-axis labels

# Add vertical lines at z = -3, -2, -1, 0, 1, 2, 3
abline(v = c(-3, -2, -1, 0, 1, 2, 3), col = "black", lty = 2)

# Add text annotations for percentages (positioned roughly as in the image)
text(0, 0.22, "68% of data", col = "black", cex = 1)  # Between -1 and 1
text(0, 0.17, "95% of data", col = "black", cex = 1)  # Between -2 and 2
text(0, 0.12, "99.7% of data", col = "black", cex = 1)  # Between -3 and 3

# Add arrows 
arrows(-1, 0.2, 1, 0.2, length = 0.1, code = 3, col = "black")  # 68%
arrows(-2, 0.15, 2, 0.15, length = 0.1, code = 3, col = "black")  # 95%
arrows(-3, 0.1, 3, 0.1, length = 0.1, code = 3, col = "black")  # 99.7%

Density function of uni-variate normal distribution

\[ f(x, \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{1}{2\sigma^2} (x - \mu)^2\right) \]

Explanation

# Both 'dnormal()' and 'dnorm()' should return identical results for the same inputs.  
dnormal=function(x=0,mean=0,sd=1){
  p=1/(sqrt(2*pi)*sd)*exp(-(x-mean)^2/(2*sd^2))
  return(p)
}
x=c(95,95,95,95)
mean=c(100,100,85,85)
sd=c(1,2,5,10)
dnormal(x,mean,sd)
## [1] 1.486720e-06 8.764150e-03 1.079819e-02 2.419707e-02
#Density function of normal distribution
dnorm(x,mean,sd)
## [1] 1.486720e-06 8.764150e-03 1.079819e-02 2.419707e-02
  1. Two variables were observed as 95 and 97 from a normal distribution. The mean and SD of the distribution are most likely to be:
mean=c(100,100,85,85)
sd=c(1,2,5,10)
x1=rep(95,4)
x2=rep(97,4)
p1=dnormal(x1,mean,sd)
p2=dnormal(x2,mean,sd)
p1*p2
## [1] 6.588916e-09 5.675558e-04 4.836409e-05 4.698734e-04

Density function of multi-variate normal distribution

  1. Three individuals with observations of 95, 100 and 70 have kinship as folloing. The population has mean of 90. Square root of genetic and residual variances are most likely to be:

\[ f(x; \mu, H) = \frac{1}{(2\pi)^{n/2} |H|^{1/2}} \exp\left(-\frac{1}{2} (x - \mu)^T H^{-1} (x - \mu)\right) \]

Explanation

dmnormal=function(x=0,mean=0,V=NULL){
  n=length(x)
  p=1/(sqrt(2*pi)^n*sqrt(det(V)))*exp(-t(x-mean)%*%solve(V)%*%(x-mean)/2)
  return(p)
}

x=matrix(c(100,95,70),3,1)
mean=rep(90,3)
K=matrix(c(1,.75,.25,.75,1,.25,.25,.25,1),3,3)

va=95
ve=5
V=2*K*va+ve
dmnormal(x,mean,V)
##              [,1]
## [1,] 6.897989e-06
va=5
ve=95
V=2*K*va+ve
dmnormal(x,mean,V)
##              [,1]
## [1,] 6.302263e-17
va=50
ve=50
V=2*K*va+ve
dmnormal(x,mean,V)
##              [,1]
## [1,] 3.255017e-06

Simulation with GAPIT

source("http://zzlab.net/GAPIT/gapit_functions.txt") 
myGD=read.table(file="http://zzlab.net/GAPIT/data/mdp_numeric.txt",head=T) 
myGM=read.table(file="http://zzlab.net/GAPIT/data/mdp_SNP_information.txt",head=T)
index1to5=myGM[,2]<6 

set.seed(99164) 
mySim=GAPIT.Phenotype.Simulation(GD=myGD[,c(TRUE,index1to5)],
GM=myGM[index1to5,],
h2=.7,
NQTN=40, 
effectunit =.95,
QTNDist="normal")

# CMLM in GAPIT
setwd("~/Desktop/CROP_545/Liang_Lab/Lab8/CMLM")
myGAPIT=GAPIT(Y=mySim$Y, GD=myGD, GM=myGM,
PCA.total=3, 
QTN.position=mySim$QTN.position, 
model="CMLM")