setwd("C:/Courses/STA4502/2009-Fall/R/")

data(swiss); swiss
x <- swiss$Infant.Mortality

# =============================
# Binomial Test for Quantiles  
# =============================


## Conover (1971), p. 97f.
## Under (the assumption of) simple Mendelian inheritance, a cross
##  between plants of two particular genotypes produces progeny 1/4 of
##  which are "dwarf" and 3/4 of which are "giant", respectively.
##  In an experiment to determine if this assumption is reasonable, a
##  cross results in progeny having 243 dwarf and 682 giant plants.
##  If "giant" is taken as success, the null hypothesis is that p =
##  3/4 and the alternative that p != 3/4.
binom.test(c(682, 243), p = 3/4)
binom.test(682, 682 + 243, p = 3/4)   



# ========================
# RUNS TEST FOR RANDOMNESS
# ========================

# YOU NEED INSTALL THE lawstat PACKAGE FIRST.
library(lawstat)
runs.test(x, plot.it=T, alternative="two.sided")

y = arima.sim(n = 20, list(ar = c(0.5)));
y; median(y)
cbind(y=round(y,2), g.l=sign(y<median(y)))
runs.test(y, plot.it=T, alternative="two.sided")

# Compute a combination number 
factorial(15)/(factorial(7)*factorial(8))


# =========================
# CHISQ TEST FOR NORMALITY
# =========================

chisq.norm <- function(x, p= c(20, 40, 60, 80)/100){
    mid.points <- quantile(x, p=p)
    cutoff <- c(min(x)-1, mid.points, max(x) +1);  
    obs.ct <- table(cut(x, breaks=cutoff, labels = NULL)); 
    label <- names(obs.ct); 
    sd.z <- (mid.points-mean(x))/sd(x); 
    p.temp <- pnorm(sd.z)
    exp.p <- c(pnorm(sd.z[1]),  p.temp[-1] - p.temp[-length(p.temp)], 
        1-pnorm(sd.z[length(sd.z)]))
    exp.ct <- length(x)*exp.p
    tabl <- data.frame(bin=label, obs.count=as.vector(obs.ct), z=c(sd.z, NA), 
        expected.prob=exp.p, expected.count=exp.ct, 
        chi.cmp=(obs.ct-exp.ct)^2/exp.ct)
    rownames(tabl) <- NULL
    print(tabl)
    chisq <- chisq.test(x=obs.ct, p = exp.p)     
    print(chisq)
}
chisq.norm(x) 

# ANOTHER VERSION AVAILABLE IN THE PACKAGE "nortest"
library(nortest)
pearson.test(x, n.classes=5, adjust = F)


# =============================================
# Q-Q Plot and Shapiro-Wilk Test for Normality
# =============================================

# Here is the normal Q-Q plot 
par(mfrow=c(1,1), mar=rep(6, 4))
qqnorm(x)
qqline(x, lwd=2, col="blue")

shapiro.test(x)



# ======================================
# Kolmogorov Test and Lilliefors Test
# ======================================

# Suppose that we want to test if the data were from Normal(20, 3)
mu0 <- 20; sigma0 <- 3
ks.test(x, y="pnorm", mean=mu0, sd=sigma0)


# A more detailed Look at the Kolmogorov test 
x.bar <- mean(x); x.bar; s <- sd(x); s

n <- length(x)
Fn <- ecdf(x)
order.stat <- knots(Fn)
Fn.ecdf <- unique(Fn(sort(x)))
Fn.1 <- c(0, Fn.ecdf[-length(Fn.ecdf)])
F.star <- pnorm(order.stat, mean=mu0, sd=sigma0)
D.plus <- F.star-Fn.ecdf; 
D.minus <- F.star-Fn.1;
work.sheet <- data.frame(order.stat=order.stat, Fn=Fn.ecdf, Fn.1=Fn.1, F.star, 
    D.plus= D.plus, D.minus= D.minus)
work.sheet


T <- c(max(abs(D.plus)), max(abs(D.minus))); T
t <- max(T); 
p.value <- 0
for (k in 1:10){
 tmp <-  2*(-1)^(k-1)*exp(-2* k^2*t^2*n)
 print(cbind(k, tmp))
 p.value <- p.value + tmp
}
p.value 


# Lilliefors Test for Normality
# -----------------------------

library(nortest)
help(package="nortest")
lillie.test(x)








#