library(pracma)

# Acumulada
fx_puntual_generica <- function(x, p) {
  return(function(t) {
    x_menor_igual_t <- x[x <= t]
    suma <- 0
    if (length(x_menor_igual_t) == 0) {
      return(0)
    }
    for (i in 1:length(x_menor_igual_t)) {
      suma <- suma + p[i]
    }
    return(suma)
  })
}

# Binomial
px_binomial_generica <- function(n, p) {
  return(function(k) {
    return(choose(n, k)*(p^k)*((1-p)^(n-k)))
  })
}

ex_binomial <- function(n, p) {
  return(n*p)
}

vx_binomial <- function(n, p) {
  return(n*p*(1-p))
}

# Poisson
px_poisson_generica <- function(n, p) {
  lambda <- n*p
  return(function(k) {
    return(exp(-lambda)*(lambda^k)/factorial(k))
  })
}

ex_poisson <- function(n, p) {
  return(n*p)
}

vx_poisson <- function(n, p) {
  return(n*p)
}

# Hipergeometrica
# n = cant_muestra
# N = cant_total
# D = cant_exitos
px_hiper_generica <- function(n, N, D) {
  return(function(k) {
    return(choose(D, k)*choose(N-D, n-k)/choose(N, n))
  })
}

ex_hiper <- function(n, N, D) {
  return(n*D/N)
}

vx_hiper <- function(n, N, D) {
  return(((N-n)/(N-1))*n*(D/N)*(1-(D/N)))
}

# Geométrica
px_geom_generica <- function(p) {
  return(function(k) {
    return(p*(1-p)^(k-1))
  })
}

px_geom_mayor <- function(k, p) {
  return((1-p)^k)
}

ex_geom <- function(p) {
  return(1/p)
}

vx_geom <- function(p) {
  return((1-p)/(p^2))
}

# Binom negativa
px_nbin_generica <- function(r,p) {
  return(function(k) {
    return(choose(k-1, r-1)*(p^r)*(1-p)^(k-r))
  })
}

ex_nbin <- function(r, p) {
  return(r/p)
}

vx_nbin <- function(r,p) {
  return(r*(1-p)/(p^2))
}

# Esperanza y Varianza puntual
ex_puntual_generica <- function(x,p) {
  return(sum(x*p))
}

vx_puntual_generica <- function(x,p) {
  ex <- ex_puntual_generica(x,p)
  return(ex_puntual_generica((x-ex)^2, p))
}

# Multinomial

px_multinomial_generica <- function(n, array_p) {
  return(function (array_x) {
    dividendo <- prod(factorial(array_x))
    multiplicando <- 1
    for (i in 1:length(array_x)) {
      multiplicando <- multiplicando * (array_p[i]^array_x[i])
    }
    return((factorial(n)/dividendo)*multiplicando)
  })
}

####################################################################################################

# Integrar
integrar <- function(fx, lower, upper) {
  return(integrate(fx, lower, upper)$value)
}

integrar_doble <- function(fx, lower, upper, lower2, upper2) {
  return(integral2(fx, lower, upper, lower2, upper2)$Q)
}

# Ex generica
ex_generica <- function(fx, lower, upper) {
  return(integrar(function(x) {x*fx(x)}, lower, upper))
}

# Vx generica
vx_generica <- function(fx, ex, lower, upper) {
  return(integrar(function(x) {fx(x)*((x-ex)^2)}, lower, upper))
}

obtener_acumulada <- function(fx) {
  return(function(t) {
    return(integrar(fx, -Inf, t))
  })
}

# Uniforme
fx_uniforme_generica <- function(a, b) {
  return(function(x) {
    ifelse(x < a | x > b, 0, 1/(b-a))
  })
}

ex_uniforme_generica <- function(a, b) {
  return((a+b)/2)
}

vx_uniforme_generica <- function(a,b) {
  return(((b-a)^2)/12)
}

# Normal
fx_normal_generica <- function(u, sg) {
  return(function(x) {
    exponent <- -((x-u)^2/(2*(sg^2)))
    return(exp(exponent)/(sg*sqrt(2*pi)))
  })
}

Fx_normal_std_generica <- function(x) {
  integrand <- function(t) {exp(-(t^2)/2)/(sqrt(2*pi))}
  return(integrar(integrand, -Inf, x))
}

Fx_normal_generica_a_std <- function(u, sgc) {
  sg <- sqrt(sgc)
  return(function(t) {
    x <- ((t-u)/sg)
    return(Fx_normal_std_generica(x))
  })
}

arg_Fx_normal_a_std <- function(u, sgc) {
  sg <- sqrt(sgc)
  return(function(p) {
    return(qnorm(p)*sg + u)
  })
}

# Gamma
ex_gamma_generica <- function(alpha, lambda) {
  return(alpha/lambda)
}

vx_gamma_generica <- function(alpha, lambda) {
  return(alpha/(lambda^2))
}

# Covarianza
cov_puntual_generica <- function(pxy, px, py, x, y) {
  suma_pxy <- 0
  suma_px <- 0
  suma_py <- 0
  index_px <- 1
  for (i in x) {
    suma_px <- suma_px + i*px[index_px]
    index_py <- 1
    for (j in y) {
      suma_pxy <- suma_pxy + i*j*pxy[index_px,index_py]
      if (i == x[1]) {
        suma_py <- suma_py + j*py[index_py]
      }
      index_py <- index_py + 1
    }
    index_px <- index_px + 1
  }
  cat("E(XY) =", suma_pxy, "\n")
  cat("E(X) =", suma_px, "\n")
  cat("E(Y) =", suma_py, "\n")
  return(suma_pxy - (suma_px * suma_py))
}

cov_continua_generica <- function(fxy, fx, fy, x1, x2, y1, y2) {
  cov_xy <- integrar_doble(function(x,y) {return(x*y*fxy(x,y))}, x1, x2, y1, y2)
  ex_x <- ex_generica(fx, x1, x2)
  ey_y <- ex_generica(fy, y1, y2)
  cat("E(XY) =", cov_xy, "\n")
  cat("E(X) =", ex_x, "\n")
  cat("E(Y) =", ey_y, "\n")
  return(cov_xy - (ex_x*ey_y))
}

# Correlacion
corr_puntual_generica <- function(pxy, px, py, x, y) {
  cov <- cov_puntual_generica(pxy, px, py, x, y)
  cat("Covarianza", cov, "\n")
  desvio_x <- sqrt(vx_puntual_generica(x, px))
  cat("Desvio X", desvio_x, "\n")
  desvio_y <- sqrt(vx_puntual_generica(y, py))
  cat("Desvio Y", desvio_y, "\n")
  return(cov / (desvio_x * desvio_y))
}

corr_continua_generica <- function(fxy, fx, fy, x1, x2, y1, y2) {
  cov <- cov_continua_generica(fxy, fx, fy, x1, x2, y1, y2)
  cat("Covarianza", cov, "\n")
  ex_x <- ex_generica(fx, x1, x2)
  desvio_x <- sqrt(vx_generica(fx, ex_x, x1, x2))
  cat("Desvio X", desvio_x, "\n")
  ey_y <- ex_generica(fy, y1, y2)
  desvio_y <- sqrt(vx_generica(fy, ey_y, y1, y2))
  cat("Desvio Y", desvio_y, "\n")
  return(cov / (desvio_x * desvio_y))
}

# Esperanza de suma
esperanza_de_suma <- function(fx, fy, x1, x2, y1, y2) {
  ex_x <- ex_generica(fx, x1, x2)
  ey_y <- ex_generica(fy, y1, y2)
  return (ex_x + ey_y)
}

varianza_de_suma <- function(fxy, fx, fy, x1, x2, y1, y2) {
  ex_x <- ex_generica(fx, x1, x2)
  vx_x <- vx_generica(fx, ex_x, x1, x2)
  cat("X: Esperanza =", ex_x, "Varianza =", vx_x, "\n")
  ey_y <- ex_generica(fy, y1, y2)
  vy_y <- vx_generica(fy, ey_y, y1, y2)
  cat("Y: Esperanza =", ey_y, "Varianza =", vy_y, "\n")
  cov_xy <- cov_continua_generica(fxy, fx, fy, x1, x2, y1, y2)
  cat("Covarianza = ", cov_xy, "\n")
  return (vx_x + vy_y + 2*cov_xy)
}

# Chebyshev
calcular_chebyshev <- function(vx, epsilon) {
  return(vx/(epsilon^2))
}

f_empirica_t <- function(cx, t) { # mean(cx <= t)
  length(cx[cx <= t]) / length(cx)
}

# Intervalos de confianza
# Longitud
long_n <- function (zalpha, sgc, n) {
  zalpha*sqrt(sgc/n)*2
}

obtener_repeticiones <- function(zalpha, sgc, longitud) {
  sgc/(longitud/(2*zalpha))^2
}

IC_varianza_conocida <- function(alpha, sgc, n, un) {
  zalpha <- qnorm(alpha/2)
  raiz <- sqrt(sgc/n)
  valor <- (zalpha*raiz)
  cat("zalpha: ", zalpha, ", raiz: ", raiz, "\n")
  cat("IC: (", un+valor, ", ", un-valor, ")\n")
}

IC_varianza_desconocida <- function(alpha, Sc, n, un) {
  zalpha <- qt(alpha/2, n-1)
  raiz <- sqrt(Sc/n)
  valor <- (zalpha*raiz)
  cat("zalpha: ", zalpha, ", raiz: ", raiz, "\n")
  cat("IC: (", un+valor, ", ", un-valor, ")\n")
}

IC_media_conocida <- function(alpha, datos, u) {
  n <- length(datos)
  suma <- sum((datos - u)^2)
  za <- qchisq(alpha/2, n, lower.tail = FALSE)
  zb <- qchisq(1-(alpha/2), n, lower.tail = FALSE)
  cat("n: ", n, "suma: ", suma, "za: ", za, ", zb: ", zb, "\n")
  cat("IC: (", suma/za, ", ", suma/zb, ")\n")
}

IC_media_desconocida <- function(alpha, datos) {
  n <- length(datos)-1
  sc <- var(datos)
  print(sum((datos-mean(datos))^2)/24)
  za <- qchisq(alpha/2, n, lower.tail = FALSE)
  zb <- qchisq(1-(alpha/2), n, lower.tail = FALSE)
  cat("n: ", n, "sc: ", sc, "za: ", za, ", zb: ", zb, "\n")
  cat("IC: (", sc*n/za, ", ", sc*n/zb, ")\n")
}

IC_exponencial <- function(alpha, datos) {
  n <- length(datos)
  suma <- sum(datos)
  za <- qchisq(1-(alpha/2), 2*n, lower.tail = FALSE)
  zb <- qchisq(alpha/2, 2*n, lower.tail = FALSE)
  cat("n: ", n, "suma: ", suma, "za: ", za, ", zb: ", zb, "\n")
  cat("IC: (", za/(2*suma), ", ", zb/(2*suma), ")\n")
}

IC_asintotico_exponencial <- function(alpha, datos) {
  n <- length(datos)
  promedio <- mean(datos)
  za <- qnorm(alpha/2)
  cat("n: ", n, "promedio: ", promedio, "za: ", za, ", zb: ", -za, "\n")
  cat("IC: (", 1/(promedio + (-za*promedio/sqrt(n))), ", ", 1/(promedio - (-za*promedio/sqrt(n))), ")\n")
}

IC_asintotico_binomial <- function(alpha, promedio, n) {
  za <- qnorm(alpha/2)
  cat("n: ", n, "promedio: ", promedio, "za: ", za, ", zb: ", -za, "\n")
  cat("IC: (", promedio + (za*sqrt(promedio*(1-promedio)/n)), ", ", promedio - (za*sqrt(promedio*(1-promedio)/n)), ")\n")
}

####################################################################################################
####################################################################################################

# Práctica 8

# 1) Sea X1, ..., Xn una m.a. de una población normal
#    a) Encontrar un int. confianza 1-alpha para la media cuando la varianza es conocida
# Si la varianza es conocida podemos estandarizar los Xi
# (un - u) / sqrt(sg^2/n) ~ Z, Z ~ N(0,1)
# P(-qnorm(alpha/2) < Z < qnorm(alpha/2)) = 1-alpha
# P(-qnorm(alpha/2)*sqrt(sg^2/n) < un - u < qnorm(alpha/2)*sqrt(sg^2/n)) = 1 - alpha
# P(-qnorm(alpha/2)*sqrt(sg^2/n) < u - un < qnorm(alpha/2)*sqrt(sg^2/n)) = 1 - alpha
# P(-qnorm(alpha/2)*sqrt(sg^2/n) + un < u < qnorm(alpha/2)*sqrt(sg^2/n) + un) = 1 - alpha
# => IC = (-qnorm(alpha/2)*sqrt(sg^2/n) + un, qnorm(alpha/2)*sqrt(sg^2/n) + un)

#    b) 20 mediciones de u. Xi = u + ei, donde ei es el error aleatorio de la i-ésima medición
# Promedio = 25.01, ei ~ N(0, 0.36)
#       i) Intervalo de confianza de nivel 0.95 para X
# Xi ~ N(u, 0.36), usamos lo anterior para obtener el IC de nivel 0.95
# alpha = 0.05 => alpha/2 = 0.025
zalpha <- qnorm(0.025)
zalpha
# P(-1.959964*sqrt(0.36/20) + 25.01 < u < 1.959964*sqrt(0.36/20) + 25.01) = 1 - 0.05 = 0.95
a <- -1.959964*sqrt(0.36/20) + 25.01
a

b <- 1.959964*sqrt(0.36/20) + 25.01
b

IC_varianza_conocida(0.05, 0.36, 20, 25.01)
IC_varianza_desconocida(0.05, 0.36, 20, 25.01)
# => IC = (24.74704, 25.27296)
long <- b-a
long

#      ii) long = b-a = 1.959964*sqrt(0.36/n) + 25.01 + 1.959964*sqrt(0.36/n) - 25.01
# long = 1.959964*sqrt(0.36/n))* 2 <= 0.05  =>  0.36/(0.05/(1.959964*2))^2 = n ~= 2212
obtener_repeticiones(1.959964, 0.36, 0.05)
long_2212 <- long_n(1.959964, 0.36, 2212)
long_2212
long_2213 <- long_n(1.959964, 0.36, 2213)
long_2213

#     iii)
obtener_repeticiones(1.959964, 0.6, 0.05) # > 3687
long_3687 <- long_n(1.959964, 0.6, 3687)
long_3687
long_3688 <- long_n(1.959964, 0.6, 3688)
long_3688

#      iv)
# Mientras más varíen los datos, más repeticiones se necesitarán para obtener un intervalo
# de confianza de igual longitud
# O dicho de otra manera, mientras más varíen los datos el IC tendrá una longitud mayor (es peor)

# 3) Sea X1, ..., Xn una m.a. de una población normal
#    a) Encontrar un IC 1-alpha para la media cuando la varianza es desconocida
# Si la varianza es desconocida podemos usar la distribución t-student
# (un - u) / sqrt(sg^2/n) ~ Z, Z ~ N(0,1) => reemplazamos sg por S
# => (un - u) / sqrt(S^2/n) ~ Z, Z ~ tn-1
# P(-qt(alpha/2, n-1) < Z < qt(alpha/2, n-1)) = 1-alpha
# P(-qt(alpha/2, n-1)*sqrt(S^2/n) < un-u < qt(alpha/2, n-1)*sqrt(S^2/n)) = 1-alpha
# P(un - qt(alpha/2, n-1)*sqrt(S^2/n) < u < un + qt(alpha/2, n-1)*sqrt(S^2/n)) = 1-alpha
# => IC = (un - qt(alpha/2, n-1)*sqrt(S^2/n), un + qt(alpha/2, n-1)*sqrt(S^2/n))

#    b) varianza desconocida, s = 0.76
#       i) Intervalo de confianza de nivel 0.95 para X = u + ei, X ~ (u, sg^2)
# Podemos estimar sg con s => usamos t-student, con n-1 = 19
zalpha <- qt(0.05/2, 20-1)
zalpha
a <- 25.01 - 2.093024*sqrt(0.76^2/20)
a

b <- 25.01 + 2.093024*sqrt(0.76^2/20)
b


# => IC = (24.65431, 25.36569)
IC_varianza_desconocida(0.05, 0.76^2, 20, 25.01)

# 4) rendimiento de almendros anual en kg = X, X ~ N(u, sg^2)
# n = 5, u = 525, s = 10
#    a) Estimar para u por IC de nivel 0.95 
n <- 5
un <- 525
s <- 10
alpha <- 0.05
IC_varianza_desconocida(alpha, s^2, n, un)
# IC: ( 512.5834 ,  537.4166 )

# 5) IC de nivel exacto 1 - alpha para la sg^2 cuando u es desconocida
# Pivot chi cuadrado con n grados de libertad
# sum((Xi - u0)^2)/sg^2 ~ Chi cuadrado con n grados de libertad
# P(z(1-alpha) < sum((Xi - u0)^2)/sg^2 < z(alpha)) = 1 - alpha => 1/sg^2 invierte la desigualdadd
# P(sum((Xi - u0)^2)/z(alpha) < sg^2 < sum((Xi - u0)^2)/z(1-alpha))
# => IC = (sum((Xi - u0)^2)/z(alpha), sum((Xi - u0)^2)/z(1-alpha))

datos_5b <- c(176.50, 191.50, 186.90, 181.10, 195.70, 188.10, 187.40, 185.10, 176.90, 191.20, 193.80, 187.00, 179.00, 173.00, 184.40, 199.60, 190.40, 206.80, 193.00, 177.10, 201.00, 192.50, 176.60, 180.10, 186.40)
IC_media_conocida(0.1, datos_5b, 185)
# IC: ( 49.97784 ,  128.7891 )

# 6) 
IC_media_desconocida(0.1, datos_5b)
# IC: ( 48.21915 ,  126.7943 )

# 8) X1, ..., Xn m.a. de E(lambda).
#    b) Hallar intervalo de confianza para lambda, de nivel exacto 1 - alpha
# P(qchisq(1-(alpha/2), 2*n, lower.tail = FALSE) < 2*lambda*sum(xi) < qchisq(alpha/2))
# P(z(1-alpha)/2*sum(xi) < lambda < z(alpha)/2*sum(xi))
# => IC = (z(1-alpha)/2*sum(xi), z(alpha)/2*sum(xi))

#    d) alpha = 0.05, n = 20
datos_4d <- c(39.08, 45.27, 26.27, 14.77, 65.84, 49.64, 0.80, 66.58, 69.60, 32.42, 228.36, 64.79, 9.38, 3.86, 37.18, 104.75, 3.64, 104.19, 8.17, 8.36)
n <- length(datos_4d)
n

suma_datos <- sum(datos_4d)
suma_datos

IC_exponencial(0.05, datos_4d)
# IC: ( 0.01242842 ,  0.03018552 )

#    e) 
# P(-z(alpha/2) <= sqrt(n)*(promedio - u) / sqrt(sg^2) <= z(alpha/2)) ~= 1 - alpha
# E(Xi) = 1/lambda, V(Xi) = 1/lambda^2 < infinito
# Por TCL => (sum(Xi) - n*u) / sqrt(sg^2*n) ~ N(0,1)
# luego, (promedio - 1/lambda)*sqrt(n) / (1/lambda)
# Como X ~ N(0,1), luego podemos usarlo como pivot
# P(-z(alpha/2) <= (promedio - 1/lambda)*sqrt(n)*lambda <= z(alpha/2)) ~= 1 - alpha
# P(-z(alpha/2) <= (promedio - 1/lambda)*sqrt(n) / promedio <= z(alpha/2))
# P(-z(alpha/2)*promedio/sqrt(n) + promedio <= 1/lambda <= z(alpha/2)*promedio + promedio)
# P(1/(z(alpha/2)*promedio/sqrt(n) + promedio) <= lambda <= 1/(-z(alpha/2)*promedio + promedio))

IC_asintotico_exponencial(0.05, datos_4d)
# IC: ( 0.01414688 ,  0.03622131 )

# 9) X1, X2, X3 => sum(Xi) ~ Bi(n, p)
# P(-zalpha <= (promedio - p) / sqrt(p*(1-p)/n) <= zalpha) => reemplazamos p por su promedioo
# P(-zalpha <= (promedio - p) / sqrt(promedio*(1-promedio)/n) <= zalpha)
# P(-zalpha*sqrt(promedio*(1-promedio)/n) + promedio <= p <= promedio + zalpha*sqrt(promedio*(1-promedio)/n))
# => IC = (-zalpha*sqrt(promedio*(1-promedio)/n) + promedio, zalpha*sqrt(promedio*(1-promedio)/n) + promedio)

IC_asintotico_binomial(0.10, 3/5, 1000)
# IC: ( 0.574518 ,  0.625482 )

long_n(qnorm(0.95), 3/5*2/5, 1000)
0.625482 - 0.574518
# Como promedio depende de n, podemos usar el peor valor de la varianza para estimar n
obtener_repeticiones(qnorm(0.95), 1/4, 0.02)
# n > 6763.859

long_n(qnorm(0.95), 1/4, 6763)
long_n(qnorm(0.95), 1/4, 6764)
IC_asintotico_binomial(0.10, 1/2, 6764)exactas20
# IC: ( 0.4900001 ,  0.5099999 )
0.5099999 - 0.4900001