# ==============================================================================
# --------------------------------- SEMINAR 9 ----------------------------------
# ==============================================================================


# ----------------------------------- TASK 1 -----------------------------------


T2019 <- c(-1.7, 1.7, 5.6, 9.4, 10.7, 20.7, 18.8, 18.9, 13.3, 9.5, 5.6, 1.9)
Taverage <- c(-2.8, -1.1, 2.5, 7.3, 12.3, 15.5, 16.9, 16.4, 12.8, 8.0, 2.7, -1.0)

# test the following hypothesis about the temperatures (using the PAIRED t-TEST)
# H0: T2019 = Taverage,
# H1: T2019 > Taverage (right-tailed alternative)

# Under the normality assumption of T2019, Taverage (check yourself!) we can test:
# H0: T2019 - Taverage = 0,
# H1: T2019 - Taverage > 0

Temp <- T2019 - Taverage

# if T2019 and Taverage are NORMALLY DISTIRBUTED, also Temp will be
# NORMALLY DISTRIBUTED => we can use simple 1-sample t-test

mu0 <- 0
alpha <- 0.05
n <- length(Temp)

mu <- mean(Temp)
sigma <- sd(Temp)

# HINT: use formulas from slaids 56-57 from LECTURE 08


# ....................................... A ....................................


# realization of the TEST STATISTIC:

t <- sqrt(n) * mu / sigma

# CRITICAL REGION:

CR1 <- qt(1 - alpha, n - 1)

c(CR1, Inf)
t

# CONCLUSION: we REJECT H0

# for VOLUNTEERS (visualization):

xx <- seq(-5, 5, 0.1)
density.t <- dt(xx, n - 1)

plot(xx, density.t, type = 'l', main = "Critical region",
     xlab = "x", ylab = "f(x)", sub = paste("mu=", mu))
rect(CR_1, 0, 5, 1, col = rgb(red = 1, green = 0, blue = 0, alpha = 0.2), lty = 0)
abline(v = t, col = 'blue', lwd = 2)


# ....................................... B ....................................


# using P-VALUE:

p.value <- 1 - pt(t, n - 1)


# ....................................... C ....................................


# using CONFIDENCE INTERVAL (left-sided):

q.t <- qt(1 - alpha, n - 1)
lower <- mu - q.t * (sigma / sqrt(n))

c(lower, Inf)


# ....................................... D ....................................


# using BUILT-IN FUNCTION:

t.test(T2019, Taverage, paired = TRUE, alternative = "greater")

# conclusion: we REJECT H0 at the significance level 0.05 in behalf of H1


# ----------------------------------- TASK 2 -----------------------------------


mach.1 <- c(29, 27, 29, 35, 29, 32, 28, 34, 32, 33)
mach.2 <- c(31, 28, 30, 28, 37, 29, 27, 27, 39, 33,
            31, 32, 31, 29, 32, 28, 27, 28, 24, 34)

alpha <- 0.05

n1 <- length(mach.1)
mu1 <- mean(mach.1)
sigma1 <- sd(mach.1)

n2 <- length(mach.2)
mu2 <- mean(mach.2)
sigma2 <- sd(mach.2)


# ....................................... A ....................................


# Test if the variances of both measurements are the same against 2-side alternative.

# HINT: use F-test for equal variances (see slaids 54-55 from LECTURE 8):
# H0: sigma1 = sigma2,
# H1: sigma1 != sigma2

# use any approach You want (from 1-3)


### 1) using TEST STATISTIC:

f <- sigma1^2 / sigma2^2

# critical region:

CR1 <- qf(alpha/2, n1 - 1, n2 - 1)
CR2 <- qf(1 - alpha/2, n1 - 1, n2 - 1)

c(0, CR1)
c(CR2, Inf)
f


### 2) using P-VALUE:

dist.f <- pf(f, n1 - 1, n2 - 1)

p.value <- 2 * min(dist.f, 1 - dist.f)


### 3) using CONFIDENCE INTERVAL:

q.f.1 <- qf(1 - alpha / 2, n1 - 1, n2 - 1)
q.f.2 <- qf(alpha / 2, n1 - 1, n2 - 1)

c(sigma1^2 / (sigma2^2 * q.f.1),
  sigma1^2 / (sigma2^2 * q.f.2))


### 4) using BUILT-IN FUNCTION:

var.test(mach.1, mach.2)

# conclusion: we do NOT REJECT H0 => we can use 2-sample T-test in (B)


# ....................................... B ....................................


# Test the hypothesis that the usual time is the same for both machines against
# 2-side alternative, assuming equal variance of the random sample distribution.

# HINT: use 2-sample T-test (see slaid 50-51 from the LECTURE 8)

# use any approach You want (from 1-3)


### 1) using TEST STATISTIC:

sigma <- sqrt(((n1 - 1) * sigma1^2 + (n2 - 1) * sigma2^2) / (n1 + n2 - 2))

t <- (mu1 - mu2) / (sigma * sqrt((n1 + n2) / (n1 * n2)))

# critical region:

CR1 <- -qt(1 - alpha/2, n1 + n2 - 2)
CR2 <- qt(1 - alpha/2, n1 + n2 - 2)

c(-Inf, CR1)
c(CR2, Inf)
t


### 2) using P-VALUE:

p.value <- 2 * min(pt(t, n1 + n2 - 2), 1 - pt(t, n1 + n2 - 2))


### 3) using CONFIDENCE INTERVAL:

q.t <- qt(1 - alpha / 2, n1 + n2 - 2)

c((mu1 - mu2) - q.t * sigma / sqrt((n1 * n2) / (n1 + n2)),
  (mu1 - mu2) + q.t * sigma / sqrt((n1 * n2) / (n1 + n2)))


### 4) using BUILT-IN FUNCTION:

t.test(mach.1, mach.2, var.equal = TRUE)

# conslusion: we do NOT REJECT H0

# if we REJECTED the equality of variances, we should use WELCHS TEST:
# t.test(mach.1, mach.2, var.equal = FALSE)


# ----------------------------------- TASK 3 -----------------------------------


# Test a hypothesis that the expected number of bombs per sector is equal to 1
# against two side alternative.

n_bombs <- c(0, 1, 2, 3, 4, 7)
n_sectors <- c(229, 211, 93, 35, 7, 1)

n <- sum(n_sectors)
events <- sum(n_bombs * n_sectors)

lambda0 <- 1
alpha <- 0.05


### 1) using a BUILT-IN FUNCTION:

# install.packages("exactci")
library(exactci)
poisson.exact(events, n, alternative = "two.sided")


### 2) using asymptotical CONFIDENCE INTERVAL:

q.n <- qnorm(1 - alpha / 2)
mu <- events / n

c(mu - q.n * sqrt(mu / n), mu + q.n * sqrt(mu / n))

# conclusion: we do NOT REJECT H0