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
dnorm(x, mean, sd)
or the custom
dnormal()
function, where \(\sigma^2 = sd^2\).# 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
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
\[ 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")