Lab8
Meijing_Liang
2025-02-26
- A variable was observed as 95 from a normal distribution. The mean and SD of the distribution are most likely to be:
- A. 100 and 1
- B. 100 and 2
- C. 85 and 5
- D. 85 and 10
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-07pnorm(-2.5, lower.tail = TRUE) # Option B: Mean = 100, SD = 2, P = 0.0062 ## [1] 0.006209665pnorm(2, lower.tail = FALSE) # Option C: Mean = 85, SD = 5, P = 0.0228 ## [1] 0.02275013pnorm(1, lower.tail = FALSE) # Option D: Mean = 85, SD = 10, P = 0.1587 ## [1] 0.1586553Empirical 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
- \(f(x, \mu, \sigma^2)\): This is
the density at a given point \(x\) for
a normal distribution with mean \(\mu\)
and variance \(\sigma^2\).
- \(\frac{1}{\sqrt{2\pi\sigma^2}}\):
The normalizing constant ensures the total area under the curve equals
1, making it a valid probability density function.
- \(\exp\left(-\frac{1}{2\sigma^2} (x -
\mu)^2\right)\): The exponential term shapes the bell curve,
measuring how far \(x\) is from the
mean \(\mu\) in terms of the standard
deviation \(\sigma\) (where \(\sigma = \sqrt{\sigma^2}\)).
- This formula is the basis for calculating densities in your earlier
R code, such as
dnorm(x, mean, sd)or the customdnormal()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- Two variables were observed as 95 and 97 from a normal distribution. The mean and SD of the distribution are most likely to be:
- A. 100 and 1
- B. 100 and 2
- C. 85 and 5
- D. 85 and 10
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-04Density function of multi-variate normal distribution
- 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:
- A. 95 and 5
- B. 5 and 95
- C. 50 and 50
- D. I do not know
\[ 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
- \(f(x; \mu, H)\): This is the
density at a given point \(x\) (a
vector of \(n\) variables) for a
multivariate normal distribution with mean vector \(\mu\) (an \(n
\times 1\) vector) and covariance matrix \(H\) (an \(n
\times n\) symmetric positive-definite matrix).
- \(\frac{1}{(2\pi)^{n/2}
|H|^{1/2}}\): The normalizing constant ensures the total volume
under the distribution equals 1, where \(n\) is the number of variables, and \(|H|\) is the determinant of the covariance
matrix \(H\). This accounts for the
dimensionality and spread of the distribution.
- \(\exp\left(-\frac{1}{2} (x - \mu)^T
H^{-1} (x - \mu)\right)\): The exponential term measures the
Mahalanobis distance between \(x\) and
\(\mu\), weighted by the inverse of the
covariance matrix \(H^{-1}\). This
captures the multivariate spread and correlation structure, generalizing
the univariate normal’s squared distance \((x
- \mu)^2 / \sigma^2\).
- This formula extends the univariate normal PDF (from your previous question) to multiple dimensions, accounting for correlations between variables via \(H\).
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-06va=5
ve=95
V=2*K*va+ve
dmnormal(x,mean,V)## [,1]
## [1,] 6.302263e-17va=50
ve=50
V=2*K*va+ve
dmnormal(x,mean,V)## [,1]
## [1,] 3.255017e-06Simulation 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")