##### Testing Alpha
#####

# This checks that the proposed test is appropriate
# according to the Type I error rates. Note that the
# test is Type-I appropriate for a given value of
# alpha if the observed rejection rate is alpha. Also
# it is appropriate for all possible values of alpha
# if the p-values follow a standard Uniform distribution.


### Preamble
source("http://rfs.kvasaheim.com/math322.R")



### The Forsberg test for kurtosis

# Generate data under the null
x = rnorm(15)

# Perform the test (and retain the p-value)
forsberg.test(x)


# Do this many, many, many times

pval = numeric()
for(i in 1:1e4) {
  x = rnorm(15)
  pval[i] = forsberg.test(x)
}


# Look at the distribution of p-values
hist(pval)

# Check that the test is fine for alpha = 0.05
alpha = 0.05
binom.test(sum(pval<alpha),length(pval), p=alpha)

# Check that it is appropriate for all alpha
ks.test(pval,"punif")






### two-sample t-test, unequal variances

pval = numeric()
for(i in 1:1e4) {
  x1 = rnorm(15, m=0,s=1)
  x2 = rnorm(15, m=0,s=5)
  pval[i] = t.test(x1,x2)$p.value
}


# Look at the distribution of p-values
hist(pval)

# Check that the test is fine for alpha = 0.05
alpha = 0.05
binom.test(sum(pval<alpha),length(pval), p=alpha)

# Check that it is appropriate for all alpha
ks.test(pval,"punif")







### two-sample t-test, unequal variances, wrong test, equal n

pval = numeric()
for(i in 1:1e4) {
  x1 = rnorm(15, m=0,s=1)
  x2 = rnorm(15, m=0,s=5)
  pval[i] = t.test(x1,x2, var.equal=TRUE)$p.value
}


# Look at the distribution of p-values
hist(pval)

# Check that the test is fine for alpha = 0.05
alpha = 0.05
binom.test(sum(pval<alpha),length(pval), p=alpha)

# Check that it is appropriate for all alpha
ks.test(pval,"punif")






### two-sample t-test, unequal variances, wrong test, UNequal n

pval = numeric()
for(i in 1:1e4) {
  x1 = rnorm(25, m=0,s=1)
  x2 = rnorm(15, m=0,s=5)
  pval[i] = t.test(x1,x2, var.equal=TRUE)$p.value
}


# Look at the distribution of p-values
hist(pval)

# Check that the test is fine for alpha = 0.05
alpha = 0.05
binom.test(sum(pval<alpha),length(pval), p=alpha)

# Check that it is appropriate for all alpha
ks.test(pval,"punif")