Week 5 Discussion
2026-10-01
Part 1
The (student) t distribution converges to normal distribution as the degrees of freedom increase (beyond 120). Please plot a normal distribution, and a few t distributions on the same chart with 2, 5, 15, 30, 120 degrees of freedom.
# Normal distribution curve
curve(dnorm(x, mean = 0, sd = 1),
from = -4,
to = 4,
main = strwrap("Convergence of T Distribution to Normal Distribution as Degrees of Freedom Increases", width = 60),
xlab = "",
ylab = "Density",
col = "black",
lwd = 2,)
# T distribution curve with 2 degrees of freedom
curve(dt(x, df = 2),
from = -4,
to = 4,
col = "red",
lwd = 2,
add = TRUE)
# T distribution curve with 5 degrees of freedom
curve(dt(x, df = 5),
from = -4,
to = 4,
col = "orange",
lwd = 2,
add = TRUE)
# T distribution curve with 15 degrees of freedom
curve(dt(x, df = 15),
from = -4,
to = 4,
col = "green",
lwd = 2,
add = TRUE)
# T distribution curve with 30 degrees of freedom
curve(dt(x, df = 30),
from = -4,
to = 4,
col = "blue",
lwd = 2,
add = TRUE)
# T distribution curve with 120 degrees of freedom
curve(dt(x, df = 120),
from = -4,
to = 4,
col = "purple",
lwd = 2,
add = TRUE)
legend("topright", legend = c("Normal", "df = 2", "df = 5", "df = 15", "df = 30", "df = 120"), col = c("black", "red", "orange", "green", "blue", "purple"), lty = 1, lwd = 2)
Part 2
Lets work with normal data below (1000 observations with mean of 108 and sd of 7.2).
Plot two charts - the normally distributed data (above) and the Z score distribution of the same data. Do they have the same distributional shape ? Why or why not ?
set.seed(123) # Set seed for reproducibility
# Defined variables
mu <- 108 # Mean
sigma <- 7.2 # Standard deviation
# Randomly generated observations
data_values <- rnorm(n = 1000, mean = mu, sd = sigma)
# Observed mean and standard deviation
observed_mean <- mean(data_values)
paste("The observed mean of the randomly generated sample is", round(observed_mean, digits = 2))## [1] "The observed mean of the randomly generated sample is 108.12"observed_sd <- sd(data_values)
paste("The observed standard deviation of the randomly generated sample is", round(observed_sd, digits = 2))## [1] "The observed standard deviation of the randomly generated sample is 7.14"# Create histogram plotting observations
hist(data_values, freq = FALSE, prob = TRUE,
breaks = "FD",
main = bquote("Distribution of 1000 Random Observations (" * bar(x) == 108.12 * "," ~ s[x] == 7.14 * ")"),
xlab = "",
xaxt = "n",
xlim = c(84, 134),
ylab = "Density",
col = "lightblue",
border = "gray")
axis(side = 1, at = seq(84, 134, by = 2))
abline(v = mu, col = "red", lwd = 1)
# Normal distribution curve
curve(dnorm(x, mu, sigma),
col = "darkgreen",
lwd = 2,
add = TRUE)
legend("topright",
legend = expression(
"Random Observations",
"Theoretical mean (" * mu == 108 * ")",
"Normal Distribution Curve"),
pch = c(15, NA, NA),
pt.cex = 2.0,
box.col = "gray",
lty = c(0, 1, 1),
lwd = c(NA, 2, 2),
col = c("lightblue", "red", "darkgreen"),
merge = FALSE,
seg.len = 1.5,
x.intersp = 0.8)
# Create histogram for z-score distribution
# Calculate z-scores for observations
z_scores <- (data_values - mu) / sigma
# Mean and standard deviation for observation z-scores
z_score_mean <- mean(z_scores)
paste("The mean z-score calculated from the randomly generated sample is", round(z_score_mean, digits = 2))## [1] "The mean z-score calculated from the randomly generated sample is 0.02"z_score_sd <- sd(z_scores)
paste("The standard deviation calculated from the randomly generated sample is", round(z_score_sd, digits = 2))## [1] "The standard deviation calculated from the randomly generated sample is 0.99"hist(z_scores, freq = FALSE, prob = TRUE,
breaks = "FD",
main = bquote("Z-score Distribution (" * bar(x)[z-score] == 0.02 * "," ~ s[x[z-score]] == 0.99 * ")"),
xlab = "",
xaxt = "n",
xlim = c(-3.5, 3.5),
ylab = "Density",
col = "salmon",
border = "brown")
axis(side = 1, at = seq(-3.5, 3.5, by = 0.5))
abline(v = 0, col = "cyan", lwd = 1.5)
# Z-score normal distribution curve
curve(dnorm(x, 0, 1),
col = "navy",
lwd = 2,
add = TRUE)
legend("topright",
legend = expression(
"Z-scores",
"Theoretical Z-score",
"mean (" * mu == 0 * ")",
"Normal Z-score",
"Distribution Curve"),
pch = c(15, NA, NA, NA, NA),
pt.cex = 2.0,
box.col = "gray",
lty = c(0, 1, 0, 1, 0),
lwd = c(NA, 2, NA, 2, NA),
col = c("salmon", "cyan", NA, "navy", NA),
merge = FALSE,
seg.len = 1.5,
x.intersp = 0.8)
Part 3.
In your own words, please explain what is p-value?
The p-value is the probability of observing a test statistic at least as extreme as the one obtained, assuming the null hypothesis is true. It is calculated as the area under the distribution curve from the test statistic outward, in one tail for a one-sided test or in both tails for a two-sided test. The significance level (\(\alpha\)) is chosen in advance and determines the critical value(s), the cutoff(s) beyond which the area equals \(\alpha\) and the null hypothesis is rejected. In a two-sided test, this area is split evenly, so each tail has area \(\alpha / 2\).
If P > \(\alpha\), the data are consistent with what we’d expect under the null hypothesis, so we fail to reject it. If P ≤ \(\alpha\), the result is statistically significant since data this extreme would be unlikely if the null were true, so we reject the null hypothesis.
# Shade normal distribution function
shadenorm = function(below=NULL, above=NULL, pcts = c(0.025,0.975), mu=0, sig=1, numpts = 500, color = "gray", dens = 40, justabove= FALSE, justbelow = FALSE, lines=FALSE,between=NULL,outside=NULL){
if(is.null(between)){
below = ifelse(is.null(below), qnorm(pcts[1],mu,sig), below)
above = ifelse(is.null(above), qnorm(pcts[2],mu,sig), above)
}
if(is.null(outside)==FALSE){
below = min(outside)
above = max(outside)
}
lowlim = mu - 4*sig # min point plotted on x axis
uplim = mu + 4*sig # max point plotted on x axis
x.grid = seq(lowlim,uplim, length= numpts)
dens.all = dnorm(x.grid,mean=mu, sd = sig)
if(lines==FALSE){
plot(x.grid, dens.all, type="l", xlab="X", ylab="Density") # label y and x axis
}
if(lines==TRUE){
lines(x.grid,dens.all)
}
if(justabove==FALSE){
x.below = x.grid[x.grid<below]
dens.below = dens.all[x.grid<below]
polygon(c(x.below,rev(x.below)),c(rep(0,length(x.below)),rev(dens.below)),col=color,density=dens)
}
if(justbelow==FALSE){
x.above = x.grid[x.grid>above]
dens.above = dens.all[x.grid>above]
polygon(c(x.above,rev(x.above)),c(rep(0,length(x.above)),rev(dens.above)),col=color,density=dens)
}
if(is.null(between)==FALSE){
from = min(between)
to = max(between)
x.between = x.grid[x.grid>from&x.grid<to]
dens.between = dens.all[x.grid>from&x.grid<to]
polygon(c(x.between,rev(x.between)),c(rep(0,length(x.between)),rev(dens.between)),col=color,density=dens)
}
}
# Example normal distribution for two-sided z-test with level of significance = 0.05
shadenorm(pcts = c(0.025, 0.975),
mu = 0,
sig = 1,
col = "red",)
title(main = bquote("Normal Distribution (" * alpha == 0.05 * ", Two-sided Z-test)"))
# P-value > 0.05
abline(v = 1.5, col = "green", lwd = 2)
# P-value < 0.05
abline(v = 2.2, col = "blue", lwd = 2)
legend("topleft",
expression(
"Rejection region",
paste("Test statistic (" * p > alpha * "),"),
paste("fail to reject" ~ H[0] * "",),
paste("Test statistic (" * p < alpha * "),"),
paste("reject" ~ H[0] * "")),
pch = c(22, NA, NA, NA, NA),
pt.bg = c(adjustcolor("red", alpha.f = 0.3), NA, NA, NA, NA),
pt.cex = 2.0,
box.col = "gray",
lty = c(0, 1, 0, 1, 0),
lwd = c(NA, 2, NA, 2, NA),
col = c("red", "green", NA, "blue", NA),
seg.len = 1.5,
merge = FALSE,
x.intersp = 0.8)