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")
}

# Tests

region_rechazo <- function(uh0, vh0, promedio, n, alpha, comparar) {
  zalpha <- qnorm(1-alpha)
  valor <- (promedio - uh0) / sqrt(vh0/n)
  cat("valor: ", valor, ", zalpha: ", zalpha, "\n")
  if (comparar(valor, zalpha)) {
    print("Rechaza H0!")
  } else {
    print("No rechaza H0")
  }
}

region_rechazo_t <- function(uh0, sc, promedio, n, alpha, comparar) {
  zalpha <- qt(1-alpha, n-1)
  valor <- (promedio - uh0) / sqrt(sc/n)
  cat("valor: ", valor, ", zalpha: ", zalpha, "\n")
  if (comparar(valor, zalpha)) {
    print("Rechaza H0!")
  } else {
    print("No rechaza H0")
  }
}

region_rechazo_chi <- function(sc, vh0, n, alpha, comparar) {
  zalpha <- qchisq(1-alpha, n-1)
  valor <- (n-1)*sc/vh0
  cat("valor: ", valor, ", zalpha: ", zalpha, "\n")
  if (comparar(valor, zalpha)) {
    print("Rechaza H0!")
  } else {
    print("No rechaza H0")
  }
}

# R = {(promedio - u) / sqrt(vh0/n) >= (uh0 - u) / sqrt(vh0/n) + zalpha}
# (promedio - u) / sqrt(vh0/n) ~ N(0,1)
# => pi(u) = 1 - pnorm((uh0 - u)/sqrt(vh0/n) + zalpha)
# P(EI) = calcular_f_potencia(...)
# P(EII) = 1 - calcular_f_potencia(...)
calcular_f_potencia_mayor <- function(uh0, vh0, uh1, n, alpha) {
  1 - pnorm((uh0-uh1)/sqrt(vh0/n) + qnorm(1-alpha))
}

# R = {(n-1)*s^2 / sg0^2 <= qchisq(alpha, n-1)}, multiplica por sg0^2 / sg1^2
# R = {(n-1)*s^2 / sg1^2 <= sg0^2*qchisq(alpha, n-1) / sg1^2}
# (n-1)*s^2 / sg1^2 ~ chisq(n-1)
# => pi(u) = pchisq(sg0^2*qchisq(alpha, n-1) / sg1^2)
calcular_f_potencia_chisq_menor <- function(sg0c, sg1c, n, alpha) {
  pchisq(sg0c*qchisq(alpha, n-1) / sg1c, n-1)
}

# pnorm((uh0-uh1)/sqrt(vh0/n) + qnorm(1-alpha)) <= valor
# (uh0-uh1) / sqrt(vh0/n) + qnorm(1-alpha) <= qnorm(valor)
# (uh0-uh1) / sqrt(vh0/n) <= qnorm(valor) - qnorm(1-alpha)
# (uh0-uh1) <= (qnorm(valor) - qnorm(1-alpha)) * sqrt(vh0/n)
# (uh0-uh1)/(qnorm(valor) - qnorm(1-alpha)) <= sqrt(vh0/n)
# ((uh0-uh1) / (qnorm(valor) - qnorm(1-alpha)))^2 <= vh0/n
# |n| >= vh0 / ((uh0-uh1) / (qnorm(valor) - qnorm(1-alpha)))^2
calcular_n_error_tipo_2 <- function(valor, uh0, vh0, uh1, alpha) {
  vh0 / ((uh0-uh1)/(qnorm(valor)-qnorm(1-alpha)))^2
}