# SEE HOW WELL A T TEST WORKS ON DATA FROM A T DISTRIBUTION.
#
# R. M. Neal, 2003.
# DO TESTS WHEN NULL HYPOTHESIS IS TRUE. Produces histograms of the
# p-values when the null hypothesis is true, for data coming from
# t distributions with various degrees of freedom, as given by the first
# argument, and for various sample sizes, as given by the second argument.
# The number of simulations to do for each test is given as the third
# argument (default 1000), and the number of bins in the histograms
# is given as the fourth argument (default 20).
do.null.tests <- function (df.list, ss.list, k=1000, bins=20)
{
# Set up to put six plots on one page.
par (mfrow=c(3,2))
# Produce histograms for each number of degrees of freedom, and each
# sample size.
for (df in df.list)
{
for (ss in ss.list)
{
hist (t.test.sim(k,ss,df), nclass=bins, prob=T,
main="", xlab=paste(k,"pvalues for n =",ss,"and df =",df))
}
}
}
# SIMULATE T TESTS ON DATA FROM A T DISTRIBUTION. Simulates k data
# sets, each consisting of n data points that are drawn independently
# from the t distribution with df degrees of freedom. For each data
# set, the p-value for a two-sided t test of the null hypothesis that
# the mean is mu (default 0) is computed. The value returned by this
# function is the vector of these k p-values.
t.test.sim <- function (k, n, df, mu=0)
{
pvalues <- numeric(k)
for (i in 1:k)
{ x <- rt (n, df)
pvalues[i] <- t.test.pvalue(x,mu)
}
pvalues
}
# FIND THE P-VALUE FOR A TWO-SIDED T TEST. The data is given by the first
# argument, x, which must be a numeric vector. The mean under the null
# hypothesis is given by the second argument, mu, which defaults to zero.
#
# Note: This function is just for illustrative purposes. The p-value
# can be obtained using the built-in t.test function with the expression:
#
# t.test(x,mu=mu)$p.value
t.test.pvalue <- function (x, mu=0)
{
if (!is.numeric(x) || !is.numeric(mu) || length(mu)!=1)
{ stop("Invalid argument")
}
n <- length(x)
if (n<2)
{ stop("Can't do a t test with less than two data points")
}
t <- (mean(x)-mu) / sqrt(var(x)/n)
2 * pt (-abs(t), n-1)
}