# Zadatak iz statistike: Neka X ima familiju Exp(lambda). Formiramo test
# kolicnika vjerodostojnosti za H0(lambda=1), H1(lambda=lambda_1), lambda_1>1.
# Za razlicite vrijednosti lambda_1, empirijski izracunati moc testa.

# Dobija se da je kriticna oblast W={sum_xk <= 83.6}


test <- function(n, lambda) {
  x <- rexp(n, lambda)
  sum(x) <= 83.6
  
}

# Moc testa za H1(lambda=1.1)

r <- replicate(10 ^ 5, test(100, 1.1))
mean(r)

# Moc testa za H1(lambda=1.5)

r <- replicate(10 ^ 5, test(100, 1.5))
mean(r)

# Moc testa za H1(lambda=2)
r <- replicate(10 ^ 5, test(100, 2))
mean(r)  # teorijski je dobijeno ~1


# Intervali povjerenja




# Neka je dat uzorak iz Exp(lambda) raspodjele. Znamo da je ocjena mmv 
# 1/sum(xk). Na osnovu asimptotske normalnosti konstruisemo intervale 
# povjerenja (izvodjenje na casu!)

# Generisacemo N uzoraka iz eksponencijalne raspodjele i naci intervale
# povjerenja, a potom izracunati procenat intervala koji su obuhvatili
# pravu vrijednost parametra lambda

ip1 <- function(n, lambda, beta) {
  x <- rexp(n, lambda)
  # mmv ocjena za lambda
  lambda.hat <- 1 / mean(x)
  
  c <- qnorm((1 + beta) / 2)
  a <- lambda.hat * 1 / (1 + c / sqrt(n))
  b <- lambda.hat * 1 / (1 - c / sqrt(n))
  # Provjeravamo da li je obuhvaceno lambda
  (a <= lambda) & (lambda <= b)
}

# izvlacimo N uzoraka i racunamo procenat

N <- 10 ^ 5

r <- replicate(N, ip1(100, 0.5, 0.95))
mean(r)
# Racunamo apsolutnu gresku u odnosu na pravu vrijednost 
err1 <- mean(r) - 0.95

# Drugi interval 

ip2 <- function(n, lambda, beta) {
  x <- rexp(n, lambda)
  # mmv ocjena za lambda
  lambda.hat <- 1 / mean(x)
  
  c <- qnorm((1 + beta) / 2)
  a <- lambda.hat * (1 - c / sqrt(n))
  b <- lambda.hat * (1 + c / sqrt(n))
  # Provjeravamo da li je obuhvaceno lambda
  (a <= lambda) & (lambda <= b)
}

# izvlacimo N uzoraka i racunamo procenat

N <- 10 ^ 5

r <- replicate(N, ip2(100, 0.5, 0.95))
mean(r)
# Racunamo apsolutnu gresku u odnosu na pravu vrijednost 
err2 <- mean(r) - 0.95

err1-err2


# lambda = 1


r1 <- replicate(N, ip2(100, 1, 0.95))
r2 <- replicate(N, ip2(100, 1, 0.95))
mean(r1)-0.95
mean(r2)-0.95

# lambda = 2

r1 <- replicate(N, ip2(100, 2, 0.95))
r2 <- replicate(N, ip2(100, 2, 0.95))
mean(r1)-0.95
mean(r2)-0.95