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

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

# Práctica 6

# 1) Dos sedes de fabrica de alfajores, Q y P. Cada sede empaqueta cajas de 4 unidades.
# X = "cantidad de alfajores defectuosos en una caja"
# Y = "sede de donde proviene la caja", Y=0 = Q, Y=1 = P

alfajores <- read.table("alfajores.txt", header = TRUE) # Descargar el archivo desde el campus
table(alfajores)

# n = 500 cajas al azar
n <- nrow(alfajores)

#    a) Probabilidad de que la caja sea de Quilmes (Y=0)
condicion_quilmes <- alfajores$fabrica == 0
n_quilmes <- nrow(alfajores[condicion_quilmes,])
n_quilmes
p_quilmes <- n_quilmes/n
p_quilmes

#    b) Probabilidad de que la caja sea de Quilmes (Y=0) y tenga 3 defectuosos (X=3)
condicion_3_defectuosos <- alfajores$defectuosos == 3
n_quilmes_3_defectuosos <- nrow(alfajores[condicion_quilmes & condicion_3_defectuosos,])
n_quilmes_3_defectuosos
p_quilmes_3_defectuosos <- n_quilmes_3_defectuosos/n
p_quilmes_3_defectuosos

#    c) f de probabilidad conjunta
pXY = t(table(alfajores)) / n
pXY

#    d) Esperanza y varianza de X (defectuosos)
ex <- mean(alfajores$defectuosos)
ex
ex <- var(alfajores$defectuosos)
ex

#    e) P(X=3 | Y=0) = pXY(3,0) / pY(0)
p_e <- p_quilmes_3_defectuosos/p_quilmes
p_e

#    f) P(X=3) = P(X=3 | Y=0) * P(Y=0) + P(X=3 | Y=1) * P(Y=1)
condicion_pilar <- alfajores$fabrica == 1
p_pilar <- nrow(alfajores[condicion_pilar,]) / n
p_pilar
p_pilar_3_defectuosos <- nrow(alfajores[condicion_pilar & condicion_3_defectuosos,]) / n
p_pilar_3_defectuosos
p_f <- p_quilmes_3_defectuosos + p_pilar_3_defectuosos
p_f
# es la marginal X = 3


# 2)

lamparas <- read.table("./lamparas.txt", header = FALSE) # Descargar el archivo desde el campus
table(lamparas)

n <- nrow(lamparas)
n

#    a) Probabilidad de que una lámpara dure más de 30 horas
condicion_mas_30_horas <- lamparas > 30
p_mas_30_horas <- length(lamparas[condicion_mas_30_horas]) / n
p_mas_30_horas

#    b) Fx empírica de este conjunto de datos
f_t <- function(cx, t) { # mean(cx <= t)
  length(cx[cx <= t]) / nrow(cx)
}
1-f_t(lamparas, 30)

# mean -> media
# median -> mediana (la mitad de ...)
# quantile -> quantiles
# mean(..., alpha) -> media podada

# dispersión
# mean(...^2)-mean(...)^2 = varianza
# var(...) = S^2 (divide por n-1), varianza muestral
# IQR(...) = distancia intercuartil, Q3 - Q1
# mad(..., constant = 1) = mediana del modulo de xi - mediana (~x)
# P(X > t) = 0.9 => P(X <= t) = 0.1 = F(t)
f_t(lamparas, 2.41)