# ZAD1

load(url("http://ls.home.amu.edu.pl/data_sets/Centrala.RData"))
library(EnvStats)

epois(Centrala$Liczba,
      method = "mle/mme/mvue",
      ci = TRUE, ci.type = "two-sided", conf.level = 0.95, 
      ci.method = "exact")$interval$limits

epois(Centrala$Liczba,
      method = "mle/mme/mvue",
      ci = TRUE, ci.type = "two-sided", conf.level = 0.95, 
      ci.method = "pearson.hartley.approx")$interval$limits

epois(Centrala$Liczba,
      method = "mle/mme/mvue",
      ci = TRUE, ci.type = "two-sided", conf.level = 0.95, 
      ci.method = "normal.approx")$interval$limits

# ZAD 2

awarie <- read.table("http://ls.home.amu.edu.pl/data_sets/awarie.txt")

lambda_limits <- eexp(awarie$V1,
                      method = "mle/mme",
                      ci = TRUE, ci.type = "two-sided", conf.level = 0.95, 
                      ci.method = "exact")$interval$limits
rev(1 / lambda_limits)
rev(1 / lambda_limits^2)


# ZAD 3

median_cint <- function(x, conf_level = 0.95) {
  m <- mean(x^2)
  a <- 1 - conf_level
  n <- length(x)
  L <- sqrt(log(2) * m * (1 - qnorm(1 - a / 2) / sqrt(n)))
  R <- sqrt(log(2) * m * (1 + qnorm(1 - a / 2) / sqrt(n)))
  w <- list(title = "mediana", est = sqrt(log(2) * m), l = L, r = R, conf_level = conf_level)
  class(w) <- "confint"
  return(w)
}


print.confint <- function(x) {
  cat(x$conf_level * 100, "percent confidence interval:", "\n")
  cat(x$l, " ", x$r, "\n")
}

summary.confint <- function(x) {
  cat("\n", "Confidence interval of", x$title, "\n", "\n")
  cat(x$conf_level * 100, "percent confidence interval:", "\n")
  cat(x$l, " ", x$r, "\n")
  cat("sample estimate", "\n")
  cat(x$est, "\n")
}



x <- c(0.9, 6.2, 2.1, 4.1, 7.3, 
       1.0, 4.6, 6.4, 3.8, 5.0, 
       2.7, 9.2, 5.9, 7.4, 3.0, 
       4.9, 8.2, 5.0, 1.2, 10.1,
       12.2, 2.8, 5.9, 8.2, 0.5)
median_cint(x)
summary(median_cint(x))


# ZAD 4

conf_int <- function(x, conf_level = 0.95) {
  c(mean(x) - sd(x) * qt(1 - (1 - conf_level) / 2, length(x) - 1) / sqrt(length(x)), 
    mean(x) + sd(x) * qt(1 - (1 - conf_level) / 2, length(x) - 1) / sqrt(length(x)))
}

mu <- 1
sigma <- 3
nr <- 1000
df <- 3
lambda <- 3

n <- 10
cat(paste("n = ", n, sep = ""), "\n")

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rnorm(n, mu, sigma)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < mu) & (mu < temp_int[2]))
}
mean(temp)

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rchisq(n, df)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < df) & (df < temp_int[2]))
}
mean(temp)

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rexp(n, rate = lambda)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < 1 / lambda) & (1 / lambda < temp_int[2]))
}
mean(temp)

n <- 50
cat(paste("n = ", n, sep = ""), "\n")

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rnorm(n, mu, sigma)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < mu) & (mu < temp_int[2]))
}
mean(temp)

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rchisq(n, df)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < df) & (df < temp_int[2]))
}
mean(temp)

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rexp(n, rate = lambda)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < 1 / lambda) & (1 / lambda < temp_int[2]))
}
mean(temp)

n <- 100
cat(paste("n = ", n, sep = ""), "\n")

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rnorm(n, mu, sigma)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < mu) & (mu < temp_int[2]))
}
mean(temp)

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rchisq(n, df)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < df) & (df < temp_int[2]))
}
mean(temp)

temp <- numeric(nr)
set.seed(1234)
for (i_nr in seq_len(nr)) {
  x <- rexp(n, rate = lambda)
  temp_int <- conf_int(x)
  temp[i_nr] <- ((temp_int[1] < 1 / lambda) & (1 / lambda < temp_int[2]))
}
mean(temp)