A one-sample t-test is conducted over a range of mean values. The null hypothesis is that test mean is the true mean of the population.
The sample is of length 500 and is a random sample based on a standard normal distribution, which has a mean of 0. Thus, we expect values far from 0 will be rejected.
If you uncomment the cat() lines, you will get intermediate results. You can use the t.test function or a calculation, based on the t-distribution. We have a two sided alternative hypotheses, that is, we sum over both tails. This is the default for a t.test.
# ex22.R
N <- 500
mu <- seq(-.2,+.2,.001) # test for the population mean
pr <- numeric(length(mu))
rej <- numeric(length(mu))
samp <- rnorm(N)
mean.samp <- mean(samp)
sd.samp <- sd(samp)
for (i in 1:length(mu)) {
m <- mu[i]
t.crit <- (mean.samp-m)/(sd.samp/sqrt(N))
if (t.crit>=0) {
p <- 2*(1-pt(t.crit, df = N-1))
} else {
p <- 2*pt(t.crit, df = N-1 )
}
pr[i] <- p
res <- t.test(samp, mu = m)
#cat('\n-- t.crit = ',t.crit,'\tres$statistic =',res$statistic)
#cat('\n-- p =',p,'\tres$p.value =', res$p.value)
if (p<0.05) {
#cat('\nmu =',m,' rejected')
rej[i] <- 1
} else {
#cat('\nmu =',m,' not rejected')
rej[i] <- 0
}
}
plot(mu, rej, type = 'l',
main = "Means rejected (p-value: blue), (threshold: green)",
col = 'red', xlab = 'mean', ylab = "Reject = 1 (red)")
lines(mu, pr, col = 'blue')
abline(h = .05, col = 'green')
Output:
No comments:
Post a Comment