###########################
# MS-C1620
# Statistical inference
# Lecture 5


# We use the significance level 0.05 throughout the script


################################################
# What is a chi-squared distribution?
#
# If X~N(0,1), then X^2 is NOT normal.
# It has "chi-squared distribution" (with parameter 1).
X = rnorm(100000)
C = X^2
hist(X,50)
hist(C,50)
X[1:5]
C[1:5]
var(C)

# If you add TWO squared standard normals, you get
# again a new distribution.  It is chi2 with parameter 2.
X = rnorm(100000)
Y = rnorm(100000)
C = X^2 + Y^2
C[1:5]
hist(C,50)
var(C)

# If you add FIVE squared standard normals, you get
# chi2 with parameter five.
X1 = rnorm(100000)
X2 = rnorm(100000)
X3 = rnorm(100000)
X4 = rnorm(100000)
X5 = rnorm(100000)
C = X1^2 + X2^2 + X3^2 + X4^2 + X5^2
C[1:5]
hist(C,50)
var(C)



################################################
# Normality testing

# Consider three data sets:
# 1. The sepal lengths of 100 irises of the species Versicolor and Virginica
# 2. The January sunspot numbers for the years 1749-1848
# 3. A random sample of size 100 from the standard normal distribution

x1 <- iris[51:150, 1]
x2 <- sunspots[seq(1, by = 12, length.out = 100)]
x3 <- rnorm(100)

# Do the samples look normal?
hist(x1, breaks = 10)
hist(x2, breaks = 10)
hist(x3, breaks = 30)


# Multiple histograms of normal samples
par(mfrow = c(4, 2))
invisible(replicate(8, hist(rnorm(100), breaks = 10)))
par(mfrow = c(1, 1))


# Bowman-Shenton/Jarque-Bera test
library(tseries)
jarque.bera.test(x1)
# No evidence against normality

jarque.bera.test(x2)
# Strong evidence aganist normality

x3 <- rnorm(1000)
jarque.bera.test(x3)


# QQ-plot
qqnorm(x1)
qqline(x1)
# The sample is slightly right-skew?

qqnorm(x2)
qqline(x2)
# The sample is heavily right-skew

qqnorm(x3)
qqline(x3)


# Shapiro-Wilk test
shapiro.test(x1)
# No evidence against normality

shapiro.test(x2)
# Strong evidence against normality

shapiro.test(x3)







#####################################################################
# Multinomial distribution


# Probability of having eye color
# BROWN | GREEN | BLUE
# 0.43  | 0.32  | 0.25

# We observe 100 people

# A sample from Multinomial(100, 0.43, 0.32, 0.25)
rmultinom(1, 100, c(0.43, 0.32, 0.25))
# Counts will typically NOT be exactly in the
# expected proportions. (It is a random sample!)


#####################################################################
# Chi-squared tests



##########################
# X^2 goodness-of-fit for discrete data

#Distribution of hair and eye color and sex in 592 statistics students.
HairEyeColor

# H0: the distribution (probabilities) of hair colors is (0.20, 0.50, 0.10, 0.20) for
# (Black, Brown, Red, Blond)?
hair <- margin.table(HairEyeColor, 1)
hair
hair/sum(hair)

# Observed and expected
hair
sum(hair)*c(0.2, 0.5, 0.1, 0.2)

chisq.test(hair, p = c(0.20, 0.50, 0.10, 0.20))
# No evidence against H0





##########################
# X^2 goodness-of-fit for discretized data

# Is the sunspot data we used before from an exponential distribution
x2 <- sunspots[seq(1, by = 12, length.out = 100)]
hist(x2, breaks = 15, probability = TRUE, ylim = c(0, 0.022))


# Estimate the parameters
lambdahat <- 1/mean(x2)
lambdahat

expfit <- function(x){ dexp(x, rate = lambdahat) }
curve(expfit, add = TRUE, col = 3, lwd = 2)


# Break the distribution into 8 categories, (0, 25, 50,...) + leftover category
hist(x2, breaks = 15, xlim = c(-1, 201), ylim = c(0, 23))
abline(v = seq(0, 200, by = 25), col = "red", lwd = 2)

# Observed frequencies
cum_x2_obs <- sapply(seq(0, 200, by = 25), function(i) sum(x2 <= i))[-1]
x2_obs <- c(cum_x2_obs[1], diff(cum_x2_obs), 0)
x2_obs

# Expected frequencies
curve(expfit, col = 3, lwd = 2, from = 0, to = 201)
abline(h = 0, lwd = 1)
segments(seq(0, 200, by = 25), rep(0, 9), seq(0, 200, by = 25), expfit(seq(0, 200, by = 25)), col = 2, lwd = 2)

p_exp <- diff(pexp(seq(0, 200, by = 25), rate = lambdahat))
x2_exp <- 100*p_exp
c(x2_exp, 100 - sum(x2_exp))

# The probabilities don't sum to one so we add one extra category
# Degrees of freedom is "categories" - 1 - "est. parameters" = 9 - 1 - 1 = 7
pchisq(chisq.test(x2_obs, p = c(p_exp, 1 - sum(p_exp)))$statistic, df = 7, lower.tail = FALSE)
# No evidence against the null that the data comes from an exponential distribution






##########################
# X^2 test of homogeneity

# This data set provides information on the fate of passengers on the fatal maiden voyage
# of the ocean liner 'Titanic', summarized according to economic status (class), sex, age and survival.
Titanic

# Let's assume that the 1st, 2nd and 3rd class ticket amounts (and the crew size) were fixed and
# each category had a specific probability to survive
Titanic_2 <- margin.table(Titanic, c(1, 4))
Titanic_2

# Is the probability to survive the same in each category?
chisq.test(Titanic_2)
# Strong evidence against H0




##########################
# X^2 test of independence

# Cross-classification of a sample of British males according to each subject's
# occupational status and his father's occupational status.
occupationalStatus
heatmap(occupationalStatus, Colv = NA, Rowv = NA, revC = TRUE)
# Think: How would the heat map look
# if the two variables are independent?

# Are the occupational statuses of a father and his son independent?
chisq.test(occupationalStatus)
# Strong evidence against H0
# (the warning has to do with some cell counts being very low,
# and the extreme classes could be combined to give a more reliable result)


##########################
# X^2 overfitting example

# Two samples, each representing counts from
# multinomial with n=1000 and p=(.25,.25,.25,.25).
p = c(.25, .25, .25, .25)

# First sample, really a random sample.
X = rmultinom(1, 1000, p)
X

# Second sample, not random but counts are "set"
# to almost equal the expected counts
Y = c(252, 248, 250, 250)
Y

chisq.test(X, p=p)
chisq.test(Y, p=p)
# Second test has statistic very close to zero,
# and very HIGH p-value. Suspicious? Was the data
# really random, and by miracle the counts were
# almost exactly the expected counts?