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

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

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

# Práctica 4

# 1) (X,Y) vector aleatorio con fx conjunta
# fxy(x,y) = { 2 <= x <= 3 && 2 <= y <= 3 ? k(x^2+y^2) : 0 }
#     a) k = 3/38
fun <- function(x, y) {  
  x^2 + y^2
}
valor_integral <- integrar_doble(fun, 2, 3, 2, 3)
k <- 1/valor_integral
k

#     b) P(2 < X < 2.6 , 2 < Y < 2.6)
fx <- function(x, y) {  
  k*(x^2 + y^2)
}
integrar_doble(fx, 2, 2.6, 2, 2.6)

#     c) Probabilidad de max(X,Y) < 2.6 ??
# => P(max(X,Y) < 2.6) = P(X < 2.6, Y < 2.6) = b)

#     d) fx = integrar(fxy(x,y) dy), fy se integra dx
xy_marginal <- function(x) {
  k*(x^2)+1/2
}

# 10)

cov_continua_generica(fx, xy_marginal, xy_marginal, 2, 3, 2, 3)
corr_continua_generica(fx, xy_marginal, xy_marginal, 2, 3, 2, 3)
esperanza_de_suma(xy_marginal, xy_marginal, 2, 3, 2, 3)
varianza_de_suma(fx, xy_marginal, xy_marginal, 2, 3, 2, 3)

# 2) 3 profesores, 2 graduados, 1 alumno. X = "numero de profes", Y = "numero de graduados"
#     a) Hallar pxy(x,y) y las marginales de X e Y

# x/y |   0   |   1   |   2   |   Px
#  0  |   0   |  2/15 |  1/15 |  1/5
#  1  |  1/5  |  2/5  |   0   |  3/5
#  2  |  1/5  |   0   |   0   |  1/5
#  3  |   0   |   0   |   0   |   0
#  Py |  2/5  |  8/15 |  1/15 |   1

# P(X=0, Y=1) =  1/6 * 2/5 + 1/3 * 1/5 = 2/15 = choose(2, 1) * choose(1,1) / choose(6,2)
# P(X=0, Y=2) =  2/6 * 1/5 = 1/15 = choose(2, 2) / choose(6,2)

# P(X=1, Y=0) =  3/6 * 1/5 + 1/6 * 3/5 = 1/5 = choose(3, 1) * choose(1,1) / choose(6,2)
# P(X=1, Y=1) =  3/6 * 2/5 + 2/6 * 3/5 = 2/5 = choose(3, 1) * choose(2, 1) / choose(6,2)

# P(X=2, Y=0) =  3/6 * 2/5 = 1/5 = choose(3, 2) / choose(6,2)


#     b) Probabilidad de que el alumno no forme parte
# P(X+Y = 2) = P(X=2, Y=0) + P(X=1, Y=1) + P(X=0, Y=2) = 1/5 + 2/5 + 1/15 = 2/3

# 10)
pxy <- matrix(c(0,2/15, 1/15, 1/5, 2/5, 0, 1/5, 0, 0, 0, 0, 0), 4, 3, byrow = TRUE)
px <- c(1/5, 3/5, 1/5, 0)
py <- c(2/5, 8/15, 1/15)

cov_puntual_generica(pxy, px, py, c(0,1,2,3), c(0,1,2))
corr_puntual_generica(pxy, px, py, c(0,1,2,3), c(0,1,2))

# 3) (X, Y) v.a. bidimensional continua, ~ U en (-1,0), (0,1), (1,1) y (2,0)
# Escrito

# 4) X = "tiempo de vida lamparita 1", Y = "lamparita 2", en 10^3 horas
# Independientes, ambas ~ Exp(1)

#     a) fxy(x,y) = fx(x)*fy(y) al ser independientes. fx(x) = lambda * exp(-lambda*x) * I(0, Inf)
# fxy(x,y) = exp(-x)*exp(-y)*I(0,Inf)(x)*I(0,Inf)(y) = exp(-x-y)*I(0,Inf)(x)*I(0,Inf)(y)

#     b) P(X+Y >= 2) = 1 - P(X+Y < 2) = 1 - P(0 < Y < 2 - X, 0 < X < 2) =>
# P(Y < 2 - X, 0 < X < 2) = integrar_doble(exp(-x-y), 0, 2, 0, 2-x)
# integrar(e^(-x)*e^(-y) dy, 0, 2-x) => e^(-x)*-e^(-y)(2-x, 0) = e^(-x)*(-e^(-2+x) + 1)
# integrar(e^(-x)*(-e^(-2+x) + 1) dx, 0, 2) => 1 - 3*e^(-2)

# => P(X+Y >= 2) = 1 - (1 - 3*e^(-2)) = 3*e^(-2)

# 5) A y B procesan a medida que van llegando. X = "tiempo que tarda A", Y = "tiempo que tarda B"
# X ~ Exp(lambda1) , Y ~ Exp(lambda2) ,  independientes
# fxy(x,y) = lambda1*lambda2*exp(-(lambda1*x)-(lambda2*y))*I(0,Inf)(x)*I(0,Inf)(y)
# P(0 < X < Inf, X < Y < Inf) = integrar_doble(lambda1*lambda2*exp(-(lambda1*x)-(lambda2*y)))

# 8) X = "número de caras obtenidas", Rx = {0,1,2,3}, X ~ Bi(3, 0.5)
# Se extraen sin reposición a+1 bolillas, Y = "número de rojas extraídas", Ry = {0,1}
# px(x) = choose(3, x) * (0.5)^3
# Y ~ Hiper(5, 1, a+1)

# Py | X = x
#  Y  |   0   |   1   |
# X=0 |  0.8  |  0.2  |
# X=1 |  0.6  |  0.4  |
# X=2 |  0.4  |  0.6  |
# X=3 |  0.2  |  0.8  |

N <- 5
D <- 1
a <- 0
py_x0 <- px_hiper_generica(a+1,N,D)
py_x0(0)
py_x0(1)

a <- 1
py_x1 <- px_hiper_generica(a+1,N,D)
py_x1(0)
py_x1(1)

a <- 2
py_x2 <- px_hiper_generica(a+1,N,D)
py_x2(0)
py_x2(1)

a <- 3
py_x3 <- px_hiper_generica(a+1,N,D)
py_x3(0)
py_x3(1)


# b)

# x/y |   0   |   1   |   Px
#  0  |  1/10 |  1/40 |  1/8
#  1  |  9/40 |  3/20 |  3/8
#  2  |  3/20 |  9/40 |  3/8
#  3  |  1/10 |  1/40 |  1/8
#  Py |  1/2  |  1/2  |   1

# c) No son independientes, pues px(x)*py(y) != pxy(x,y), como por ejemplo, 1/2*1/8 != 1/10

# d) P(X = 2 | 2B) = P(2B | X=2)*P(X=2) / P(2B) = P(Y = 1 | X = 2) * P(X=2) / P(2B)
# P(2B) = P(Y = 0 | X = 1) + P(Y = 1 | X = 2) = (9/40 * 2) * P(X=2)
p_2b <- py_x2(1) + py_x1(0)
p_x2_2b <- py_x2(1) / p_2b
p_x2_2b


# 14) Longitud 100 caracteres, pueden ser A (0.7), B (0.12), C (0.1) o D (0.08).
#     Suponiendo independencia,
#     a) Distribución de (A,B,C,D)
# Cada experimento mide el "número de veces que ocurre", (A,B,C,D) ~ Multi(100, pA, pB, pC, pD)

#     b) distrib. marginales de A, B, C y D
# ~ Bi(100, p_exito)


# 15) L = "longitud de varilla" ~ N(5, 0.25), 40 varillas al azar
# A = "mide menos de 4cm", B = "entre 4 Y 4.8cm", C = "entre 4,8 y 5,5cm", D = "más de 5,5cm"

# X ~ Multi(pA, pB, pC, pD)
# pA = P(L < 4) = Fx(4)
Fx <- Fx_normal_generica_a_std(5, 0.25)
pA <- Fx(4)
# pB = P(4 < L < 4.8) = Fx(4.8) - Fx(4)
pB <- Fx(4.8) - Fx(4)
# pC = P(4.8 < L < 5.5) = Fx(5.5) - Fx(4.8)
pC <- Fx(5.5) - Fx(4.8)
# pD = P(L > 5.5) = 1 - Fx(5.5)
pD <- 1 - Fx(5.5)

px_multi <- px_multinomial_generica(40,c(pA,pB,pC,pD))
px_multi(c(1,13,18,8))

# 19) 14 ratas especie A, 11 especie B. Proba de ue muera cualquiera es 0.1, independiente
#     a) P(A+B > 4), A = "cantidad de ratas que murieron de A", B = "... de B"
# A ~ Bi(14, 0.1) , B ~ Bi(11, 0.1)
# P(A+B > 4) = 1 - P(A+B <= 4) = 1 - Fa+b(4), A+B ~ Bi(14+11, 0.1)

px <- px_binomial_generica(14+11,0.1)
rx <- 0:(14+11)
fx <- fx_puntual_generica(rx, sapply(rx, px))
1-fx(4)

#     b) P(A | A + B = 2) = pxy(x,y) / px(x)
px_a <- px_binomial_generica(14,0.1)
px_b <- px_binomial_generica(11,0.1)
px_a(2)*px_b(0) / px(2)