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

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

# Práctica 7

# 1) 12 piezas de pan blanco. X = "Porcentaje de carbohidratos"
carbohidratos <- c(76.93, 76.88, 77.07, 76.68, 76.39, 75.09, 77.67, 76.88, 78.15, 76.5, 77.16, 76.42)
n <- length(carbohidratos)
n

#    a) Estimar la esperanza
a <- mean(carbohidratos)
a

#    b) Estimar la mediana
b <- median(carbohidratos)
b

#    c) Estimar la prob que X sea menor o igual a 76.5
cond <- carbohidratos <= 76.5
c <- length(carbohidratos[cond])/n
c


# 2) m.a de dist N(u, sigma^2). Obtener estimadores de máxima verosimilitud
# func de verosimilitud, L(u) = productoria de las p(xi, tita)
# 1/sg*sqrt(2*pi) * exp(-(x-u)^2/(2*sg^2))
#    a) u, siendo sg^2 = sg^2 obs, conocida
# si sg^2 es conocida, luego (1/sg*sqrt(2*pi))^n es constante
# luego => L(u, sg^2 ; x) = cte * exp(-sumatoria((x-u)^2)/(2*sg^2))
# => l(u, sg^2 ; x) = cte - (sumatoria((x-u)^2)/(2*sg^2))
# maximizamos derivando e igualando a 0
# minimizar h(u) = sumatoria((x-u)^2)
# => u = promedio de X

#    b) sg^2, siendo u = u obs, conocida
# ...
#    c) (u,sg^2)
# ...


# 3) a) Máquina envasa caramelos, peso en gr de cada bolsa es v.a ~ N(u, sg^2)
bolsas <- c(210, 197, 187, 217, 194, 208, 220, 199, 193, 203, 181, 212, 188, 196, 185)
n <- length(bolsas)
n

u <- mean(bolsas)
u

v <- mean(bolsas^2) - mean(bolsas)^2
v
v2 <- sum((bolsas - mean(bolsas))^2) / n
v2
v3 <- var(bolsas) * (n-1) / n
v3

#    b) Se realizan 20 mediciones de u. Cada observación X = u + e, donde e es el error aleatorio
# Suponiendo e ~ N(0, 0.01), estimar u
observaciones <- c(25.11, 25.02, 25.16, 24.98, 24.83, 25.05, 24.94, 25.04, 24.99, 24.96, 25.03, 24.97, 24.93, 25.12, 25.01, 25.12, 24.90, 24.98, 25.10, 24.96)
n <- length(observaciones)
n

u <- mean(observaciones)
u

sum((observaciones - u)^2)/n
var(observaciones)*(n-1)/n

#    c)
obs_c <- c(12.51, 11.66, 11.91, 12.25, 11.54, 11.36, 12.40, 12.19, 12.88, 12.16, 12.69, 12.91, 12.12, 11.02, 12.53, 11.77, 12.72, 10.56, 11.52, 11.66, 12.25, 12.09, 11.48, 12.36)
n_c <- length(obs_c)
n_c
sum((obs_c - 12)^2) / n_c
var(obs_c)*(n_c-1) / n_c


# 4) X = "Duracion lampara" ~ Exp(tita). n lamparas al azar, X1, ..., Xn
#    a.i) Estimador de momentos y de maxima verosimilitud de tita
# Momentos => momento poblacional = momento muestral
# => E(X) = sum(Xi)/n
# E(X) = 1/tita
# 1/tita = sum(Xi)/n => tita = n/sum(Xi)

# MV
# => L(tita) = product(tita * exp(-tita*xi)) = tita^n * exp(-tita*sum(xi))
# => log(L(tita)) = log(tita^n * exp(-n*tita*xi)) = n*log(tita) - tita*sum(xi)
# => (log(L(tita)))' = 0 para encontrar el maximo (respecto de tita)
# n/tita - sum(xi) = 0 => n/tita = sum(xi) => tita = n/sum(xi)

#     ii) Del percentil 90 (cuantil 0.9)

# 90 percentil
# P(X <= t) = 0.9 = 1 - exp(-tita*t)
# exp(-tita*t) = 0.1
# -tita*t = ln(0.1) = -ln(10)
# t = -ln(10)/-tita = ln(10)/tita
# t = ln(10)*sum(xi)/n

#     d)

lamparas_4 <- 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(lamparas_4)
n

tita <- n/sum(lamparas_4)
tita

# P(X <= x_90) = 1-e^(-lambda*x_90)
x_90 <- log(10)/tita
x_90

# f_t <- function(cx, t) { # mean(cx <= t)
#   length(cx[cx <= t]) / nrow(cx)
# }
#
# => f_t(x_90) = 0.9 => 0.9*n = length(cx[cx <= t])
# 0.9 * 20 = 18
# ordenando de menor a mayor, vemos que x_90 = 104.19
mean(lamparas_4 <= 104.19)

#    e)

lamparas_4_e <- scan("lamparas.txt") # Descargar archivo desde el campus
lamparas_4_e

n <- length(lamparas_4_e)
n

tita_e <- n/sum(lamparas_4_e)
tita_e

# P(X <= x_90) = 1-e^(-lambda*x_90)
x_90_e <- log(10)/tita_e
x_90_e

x_90_e_empirico <- sort(lamparas_4_e)[0.9*length(lamparas_4_e)]
x_90_e_empirico
mean(lamparas_4_e <= x_90_e_empirico)

x_90_e / x_90_e_empirico


# 5) X = "llamadas que recibe una central en un día" ~ P(tita)
# Estimar P(X = 40), se registran llamadas durante 48 dias
# X1, X2, ..., X48

# Momentos => E(Xi) = sum(Xi)/48
# E(X) = tita
# tita = sum(Xi)/48

# MV => L(k, tita) = prod(tita^k*exp(-tita)/k!)
# tita^(sum(k)) * exp(-tita*n) * prod(1/k!)
# ln(tita^(sum(k))) + ln(exp(-tita*n)) + ln(prod(1/k!))
# derivando e igualando a 0 =>
# sum(k)/tita - n = 0 => tita = sum(k)/n => igual a momentos

# P(X = 40) = tita^40*exp(-tita)/40!
# Como no sabemos tita, usamos la estimacion tita = sum(k)/48
# P(X = 40) = (sum(k)/n)^40 * exp(-sum(k)/n)/40!

#     d)
llamadas_5_d <- c(40, 46, 46, 45, 42, 44, 50, 31, 41, 42, 50, 34, 62, 32, 46, 50, 39, 42, 44, 41, 47, 42, 41, 50, 34, 47, 38, 40, 44, 45, 35, 51, 38, 41, 39, 34, 48, 35, 40, 40, 43, 36, 40, 49, 45, 47, 34, 45)
n <- length(llamadas_5_d)
n

tita <- sum(llamadas_5_d)/n
tita

p_x_40 <- tita^40*exp(-tita)/factorial(40)
p_x_40

llamadas_5_d_200 <- scan("llamadas.txt") # Descargar archivo desde el campus
n_200 <- length(llamadas_5_d_200)
n_200

tita_200 <- sum(llamadas_5_d_200)/n_200
tita_200

p_x_40_200 <- tita_200^40*exp(-tita_200)/factorial(40)
p_x_40_200


# 6) f(x;tita) = 1/tita*x^(-1+1/tita)*I(0,1)(x)
# Momentos => integral de x*1/tita*x^(-1+1/tita) entre 0 y 1 = 1/(tita+1)
# 1/(tita+1) = sum(xi)/n
# tita = n/sum(xi) - 1


# MV => L(xi, tita) = prod(1/tita*xi^(-1+1/tita))
# (1/tita)^n * prod(xi^(-1+1/tita))
# log((1/tita)^n) + log(prod(xi^(-1+1/tita)))
# n*log(1/tita) + (-1+1/tita)*sum(log(xi))
# -n/tita = sum(log(xi))/
# => tita = -sum(log(xi)) / n

# es insesgado

longitudes <- c(0.5, 0.7, 0.8, 0.95, 0.9, 0.6, 0.2, 0.85, 0.3, 0.2, 0.76, 0.55, 0.48, 0.8, 0.76, 0.13, 0.15, 0.67, 0.9, 0.95)
n <- length(longitudes)
n

tita <- n / sum(longitudes) - 1
tita

tita_mv <- -(sum(log(longitudes))) / n
tita_mv