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

fx_inversa_generica <- function(fx) {
  suma <- 0
  px <- c()
  for (acum in fx) {
    px <- append(px, acum-suma)
    suma <- acum
  }
  return(px)
}

ex_generica <- function(x,p) {
  return(sum(x*p))
}

vx_generica <- function(x,p) {
  ex <- ex_generica(x,p)
  return(ex_generica((x-ex)^2, p))
}

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

# 1) 15 artículos, 4 son defectuosos. Se eligen 3 con reposición. Si X es nº defectuosos,
#    a) Hallar la función de prob. puntual asociada a X y graficar
#       px(t) = P(X = t) = F(t) - F(t-)
#       Rango de X => {0, 1, 2, 3}
#       p(def) = 4/15    p(def(c)) = 11/15
#       px(0) = P(X = 0) = "no salio ninguno defectuoso" = 11/15^3
#       px(1) = P(X = 1) = "salio uno defectuoso" = (11/15^2 * 4/15) * 3
#       px(2) = P(X = 2) = "salieron dos defectuosos" = (11/15 * 4/15^2) * 3
#       px(3) = P(X = 3) = "salieron todos defectuosos" = 4/15^3
#       t={0,1,2,3}

p_def <- 4/15
p_no_def <- 11/15
x <- c(0,1,2,3)
p_x <- c(p_no_def^3, (p_no_def^2*p_def)*3, (p_no_def*p_def^2)*3, p_def^3)
p_x

p_x_2 <- c()
for (k in x) {
  comb_n_k = factorial(3) / (factorial(3-k)*factorial(k))
  p_x_2 <- append(p_x_2, comb_n_k*(p_def^k)*(1-p_def)^(3-k))
}
p_x_2

plot(x, p_x)

#    b) Cuál es la prob. de que al menos 2 artículos sean defectuosos?
#       P(X >= 2) = P(DDF) + P(DFD) + P(FDD) + P(DDD) = P(X=2) + P(X = 3)
#       o también, por complemento, donde Y es nº de no defectuosos
#       P(Y < 2) = P(Y <= 1) = P(X=0) + P(X=1)
#       P(X >= 2) = 1-P(Y < 2) = 1-P(X=0) + P(X=1)
p_min_defectuosos <-  ((p_no_def*p_def^2)*3) + p_def^3
p_min_defectuosos
#    c) Hallar la función de dist. acumulada de X y graficar
#       Fx(t) = {
#         0,      si t < 0
#         0.394,  si 0 < t <= 1
#         0.824,  si 1 < t <= 2
#         0.98,   si 2 < t <= 3
#         1,      si t > 3
#       }

acumulada <- function(t) {
    if (t < 0) {
      return(0)
    }
    if (t <= 1) {
      return(0.394)
    }
    if (t <= 2) {
      return(0.824)
    }
    if (t <= 3) {
      return(0.98)
    }
    return(1)
}

x <- -2:6
plot(x, sapply(x, acumulada))

#    d) Estimar simulando las probs del item a)
# defectuoso = 0, no defectuoso = 1
def <- 0
no_def <- 1
articulos <- c(rep(def, 4), rep(no_def, 11))

extraer_articulos = function(n) {
  sample(articulos, n, replace = TRUE)
}

matriz <- matrix(replicate(1000, extraer_articulos(3)), 3)
total <- ncol(matriz)
sin_defectuosos <- matriz[,matriz[1,] != 0 & matriz[2,] != 0 & matriz[3,] != 0]
condicion_un_def <- (matriz[1,] == 0 & matriz[2,] != 0 & matriz[3,] != 0) | (matriz[1,] != 0 & matriz[2,] == 0 & matriz[3,] != 0) | (matriz[1,] != 0 & matriz[2,] != 0 & matriz[3,] == 0)
un_defectuoso <- matriz[,condicion_un_def]
condicion_dos_def <- (matriz[1,] == 0 & matriz[2,] == 0 & matriz[3,] != 0) | (matriz[1,] == 0 & matriz[2,] != 0 & matriz[3,] == 0) | (matriz[1,] != 0 & matriz[2,] == 0 & matriz[3,] == 0)
dos_defectuosos <- matriz[,condicion_dos_def]
tres_defectuosos <- matriz[,matriz[1,] == 0 & matriz[2,] == 0 & matriz[3,] == 0]

# P(X=0) = 0.394
p_x0 <- 0.394
p_experimento_0 <- ncol(sin_defectuosos) / total
p_experimento_0 / p_x0

# P(X=1) = 0.43
p_x1 <- 0.43
p_experimento_1 <- ncol(un_defectuoso) / total
p_experimento_1 / p_x1

# P(X=2) = 0.156
p_x2 <- 0.156
p_experimento_2 <- ncol(dos_defectuosos) / total
p_experimento_2 / p_x2

# P(X=3) = 0.019
p_x3 <- 0.019
p_experimento_3 <- ncol(tres_defectuosos) / total
p_experimento_3 / p_x3

fx_inversa_generica(c(0.394, 0.824, 0.980, 1))

p_experimento_0+p_experimento_1+p_experimento_2+p_experimento_3 / 1

# 5) Esperanza => E(x) = sum(xi*px(xi)) desde i = 1 hasta inf
e_x1 <- p_x0*0 + p_x1*1 + p_x2*2 + p_x3*3
e_x1

e_simulada_x1 <- p_experimento_0*0 + p_experimento_1*1 + p_experimento_2*2 + p_experimento_3*3
e_simulada_x1/e_x1


# 6)
ex_1 <- ex_generica(c(0,1,2,3), c(p_x0,p_x1,p_x2,p_x3))
ex_1

# 7) 
vx_generica(c(0,1,2,3), c(p_x0,p_x1,p_x2,p_x3))
p_x0*(0-ex_1)^2 + p_x1*(1-ex_1)^2 + p_x2*(2-ex_1)^2 + p_x3*(3-ex_1)^2

# 2)
# X VA tiene Fx(x) = {
#   0,     si x < 1
#   0.3,   si 1 <= x < 3
#   0.4,   si 3 <= x < 6
#   0.6,   si 6 <= x < 12
#   1,     si 12 <= x
# }
# 
# Rx = {1, 3, 6, 12}
# Usar que px(t) = Fx(t) - Fx(t-)
# Fx(t) = sum(px(xi)) para xi <= t
# Fx(0) = px(0) = 0
# Fx(1) = Fx(0) + px(1) = 0.3 => px(1) = 0.3
# Fx(3) = Fx(1) + px(3) = 0.4 => px(3) = 0.1
# Fx(6) = Fx(3) + px(6) = 0.6 => px(6) = 0.2
# Fx(12) = Fx(6) + px(12) = 1 => px(12) = 0.4

# b)
# P(3 < X <= 6) = {
#   - con px => P(3 < X <= 6) = P(X<=6)-P(X<=3) = px(6)+px(3)+px(1) - px(3)+px(1) = 0.2 = px(6)
#   - con Fx => Fx(6) - Fx(3) = 0.6 - 0.4 = 0.2
# }

# P(3<= X <= 6) = {
#   - con px => px(6) + px(3) = 0.3
#   - con Fx => Fx(6) - Fx(3-) = 0.6 - 0.3 = 0.3
# }

# P(X >= 4) = P(X > 3) {
#   - con px => px(6) + px(12) = 0.2 + 0.4 = 0.6
#   - con Fx => 1 - P(X <= 3) = 1-Fx(3) = 1-0.4 = 0.6
# }

# P(X >= 6) = 1-P(X < 6) {
#   - con px => px(6) + px(12) = 0.2 + 0.4 = 0.6
#   - con Fx => 1 - P(X < 6) = 1-Fx(6-) = 1-0.4 = 0.6
# }

# 3)

x <- c(1,3,6,12)
px <- c(0.3, 0.1, 0.2, 0.4)

fx_2 <- fx_generica(x, px)
lapply(x, fx_2)

simular_2 <- function(n) {
  return(sample(x, n, replace = TRUE, prob = px))
}

# c)
simular_2(5)

# d)
data <- simular_2(1000)
total <- length(data)
# i)
i <- data[data > 3 & data <= 6]
prob_i <- length(i)/total

prob_i / 0.2

# ii)
ii <- data[data >= 3 & data <= 6]
prob_ii <- length(ii)/total

prob_ii / 0.3

# iii)
iii <- data[data >= 4]
prob_iii <- length(iii)/total

prob_iii / 0.6

# iv)
iv <- data[data >= 6]
prob_iv <- length(iv)/total

prob_iv / 0.6


# ver 5), 6) y 7) más arriba

# 8) X ~ B(p)
# Bernoulli => X toma los valores 1 (éxito) o 0 (fracaso), el parámetro p es la P(éxito)

# b)
# V(x) = E((x-E(x))^2) = (0-p)^2*(1-p) + (1-p)^2*(p) = p^2*(1-p) + (1-p)^2*p = (1-p)*(p^2 + p(1-p))
# = (1-p)*(p*((1))) = (1-p)*p

# 9) 70% de consultas requiere acceso a db. Un sistema recibe 25 consultas independientes
#    a) prob. de que exactamente 20 requieran acceso a db
# X = "cantidad de consultas que requieren acceso a db"
# Rx = {0,..., 25}
# X ~ Bi(n = 25, p = 0.7)
# P(X = 20) = px(20) = (25 20)*p^20*(1-p)^(5)
px_ej9 <- px_binomial_generica(25, 0.7)
p_x_20 <- 53130*((0.7)^20)*(0.3^5)
p_x_20
px_ej9(20)

p_x_21 <- 12650*((0.7)^21)*(0.3^4)
p_x_21
px_ej9(21)

p_x_22 <- 2300*((0.7)^22)*(0.3^3)
p_x_22
px_ej9(22)

p_x_23 <- 300*((0.7)^23)*(0.3^2)
p_x_23
px_ej9(23)

p_x_24 <- 25*((0.7)^24)*(0.3^1)
p_x_24
px_ej9(24)

#    ii) el nro de consultas que requieran acceso a db esté entre 20 y 24 inclusive
# P(20 <= X <= 24) = Fx(24) - Fx(20-) = px(20) + px(21) + px(22) + px(23) + px(24)
p_20_x_24 <- p_x_20 + p_x_21 + p_x_22 + p_x_23 + p_x_24
p_20_x_24

# con Fx
x <- 0:25
probs_ej9 <- sapply(x, px_ej9)
fx_ej9 <- fx_generica(x, probs_ej9)
p_20_x_24_2 <- fx_ej9(24) - fx_ej9(19)
p_20_x_24_2

#     b) Calcular esperanza y varianza
e_x <- ex_generica(0:25, probs_ej9)
e_x

e_x2 <- ex_binomial(25, 0.7)
e_x2

v_x <- vx_generica(0:25, probs_ej9)
v_x


plot(x, probs_ej9)
points(e_x,0, col="red", pch=16)


# 10) Dos dados, uno equilibrado y otro cargado (1 y 2 con prob 1/3, resto 1/12)
#     Un dado al azar 3 veces (independientemente). X = "numero de veces que sale 1 ó 2"
#     a) Si eligió el cargado, distrib de X
#        Rx = {0, 1, 2, 3}
n <- 3
r_x <- 0:n
p_exito <- 2/3
px <- px_binomial_generica(n, p_exito)
fx <- fx_generica(r_x, sapply(r_x,px))
plot(r_x, sapply(r_x, fx), type="s")

## Acumulada para el dado normal
# p_exito <- 1/3
# px <- px_binomial_generica(n, p_exito)
# fx <- fx_generica(r_x, sapply(r_x,px))
# plot(r_x, sapply(r_x, fx), type="s")


# 11) Se inspeccionan 10 de 50 industrias. X="cantidad de ind. visitadas en infracción"
#   a) i) 35 no cumplen, cuál es la dist de visitadas en infracción? 
# X ~ Bi(10, 0.7)
# Rx = {0, 1, ..., 10}
n <- 10
r_x <- 0:n
p_exito <- 35/50
px <- dhyper(r_x, 35, 15, n)
fx <- fx_generica(r_x, px)
plot(r_x, px, type="s")

ex <- ex_generica(r_x, px)
points(ex,0, col="red", pch=16)

vx <- vx_generica(r_x, px)


#       ii) P(X=6)?
px[7]
px_hiper_generica(10, 50, 35)(6)

#   b) 500 inds, 350 no en infracción.
# para aproximar una hiper usamos la V.A binomial
p_exito <- 350/500
px_2 <- px_binomial_generica(n, p_exito)
list_px_2 <- sapply(r_x, px_2)
fx_2 <- fx_generica(r_x, list_px_2)
plot(r_x, list_px_2, type="s")

ex_aprox <- ex_generica(r_x, list_px_2)
points(ex_aprox,0, col="red", pch=16)

vx_aprox <- vx_generica(r_x, list_px_2)

#       ii) P(X=6)?
px_2(6)

# c)

ex
ex_hiper(10, 50, 35)
ex_aprox
e_x2 <- ex_binomial(10, p_exito)
e_x2


vx
vx_hiper(10, 50, 35)
vx_aprox


# 12) Ruleta, 38 secciones. 18 rojas/negras y 2 verdes. X = "nro. necesario de juegos hasta verde"
#     X ~ G(p = 2/38 = 1/19)
#     Rx = {1, 2, ..., n}
#     a) Prob de necesitar al menos 4 jugadas => P(X >= 4) = P(X > 3) = (1-2/38)^3
p_exito <- 2/38
px <- px_geom_generica(p_exito)
list_px <- sapply(1:420, px)
fx <- fx_generica(1:420, list_px)
p_mayor_3 <- 1-fx(3)
p_mayor_3

#     b) Fx(x)
plot(1:420, sapply(1:420, fx))

#     c) Si necesité 7 o más jugadas, P(X >= 10 | X >= 7) ? Comparar con a)
# Falta de memoria => P(X >= 10 | X >= 7) = P(X > 9 | X > 6) = P(X > 3)
p_mayor_3
#     d) Prob de necesitar num. impar de jugadas
# px(k) con k impar = px(2k-1) = sum (1-2/38)^(2k-1)*2/38  k=1 to inf
# 18/37
#     e) Ex y Vx
ex <- ex_generica(1:420, list_px)
ex
ex2 <- ex_geom(p_exito)
ex2


vx <- vx_generica(1:420, list_px)
vx
vx2 <- vx_geom(p_exito)
vx2


# 13) Y = "Jugadas hasta obtener exactamente 3 verdes"
#     Y ~ Bn(3, 2/38)
#     b) P(X = 5)
px <- dnbinom(0:2, 3, 2/38)
px[length(px)]

px_2 <- px_nbin_generica(3, 2/38)
px_2(5)

ex <- ex_nbin(3, 2/38)
ex
ex_2 <- ex_generica(1:1000, sapply(1:1000, px_2))
ex_2
ex_3 <- ex_generica(1:1000, dnbinom(0:999, 3, 2/38))
ex_3

vx <- vx_nbin(3, 2/38)
vx
vx_2 <- vx_generica(1:1000, sapply(1:1000, px_2))
vx_2
vx_3 <- vx_generica(1:1000, dnbinom(0:999, 3, 2/38))
vx_3


# 14) Encontrar palabra clave, elige 1 de 2 secuencias al azar.
#     10% de la primer secuencia contiene la palabra clave
#     5% de la segunda secuencia contiene la palabra clave
#     X = "cantidad de sitios explorados hasta encontrar uno con la palabra clave, primera sec"
#     Y = "..., segunda sec"
#     Rx = {1, 2, ..., inf}
#     X ~ G(0.1)
#     Y ~ G(0.05)
#     a) Termina al encontrar la palabra clave, P(X > 5)?
# P(X > 5) = P(X > 5 | sec1)*0.5 + P(X > 5 | sec2)*0.5
p_exito_sec1 <- 0.1
p_exito_sec2 <- 0.05
px_sec1 <- px_geom_generica(p_exito_sec1)
fx <- fx_generica(1:100,sapply(1:100, px_sec1))
plot(1:100, sapply(1:100, fx))
ex <- ex_geom(p_exito_sec1)
ex
abline(v=ex)
1-fx(5)
# o tambien, con la formula
px_5_sec1 <- px_geom_mayor(5, p_exito_sec1)
px_5_sec1

px_sec2 <- px_geom_generica(p_exito_sec2)
fx_2 <- fx_generica(1:100,sapply(1:100, px_sec2))
1-fx_2(5)
px_5_sec2 <- px_geom_mayor(5, p_exito_sec2)
px_5_sec2

px_5 <- (px_5_sec1+px_5_sec2)*0.5
px_5

#     b) Encontró en X = 6, P(segunda secuencia)?
# P(sec2 | X = 6) = P(X=6 | sec2) / (P(X=6 | sec1)+P(X=6 | sec2))
px_6_sec2 <- px_sec2(6)
px_6_sec1 <- px_sec1(6)
p_sec2_dado_6 <- px_6_sec2 / (px_6_sec2+px_6_sec1)
p_sec2_dado_6

#     c) Termina cuando se encuentran 2 sitios con la palabra clave, P(X = 10)?
# W = "cantidad de sitios explorados hasta encontrar 2 con la palabra clave, sec1"
# W~Bn(2,0.1)
# Rw = {1, 2, ..., inf}
# Z = "cantidad de sitios explorados hasta encontrar 2 con la palabra clave, sec2"
# Z~Bn(2,0.05)
# Rz = {1, 2, ..., inf}
# P(X = 10) = P(X=10 | sec1) * 0.5 + P(X=10 | sec2) * 0.5
pw <- px_nbin_generica(2, 0.1)
pw_10 <- pw(10)

pz <- px_nbin_generica(2, 0.05)
pz_10 <- pz(10)

p_10 <- (pw_10 + pz_10) * 0.5
p_10

# 15) Demanda semanal de cajones ~ P(2), 4 cajones al principio de la semana
#     X = "cantidad de cajones pedidos en la semana"
#     X~P(2)
#     Rx={0,1,...,inf}
#     a) Probabilidad de vender todo su stock durante la semana => P(X >= 4) = 1-P(X <= 3)
rx <- 0:100
px <- px_poisson_generica(1, 2)
list_px <- sapply(rx, px)
fx <- fx_generica(rx, list_px)
p_4 <- 1-fx(3)
p_4

#     b) Probabilidad de no poder cumplir con al menos un pedido => P(X > 4) = 1-P(X <= 4)
1-fx(4)

#     c) Distribución de cajones vendidos en una semana
# Y = "cantidad de cajones vendidos por semana
# Y ~ P(2)
# Ry = {0, 1, 2, 3, 4}
rx <- 0:4
px <- px_poisson_generica(2,1)
list_px <- sapply(rx, px)
fx <- fx_generica(rx, c(list_px[1:4],p_4))
plot(rx, sapply(rx, fx))

#     d) Cuantos cajones necesita / P(X <= k) >= 0.99
# P(X <= k) = Fx(k) = 0.99 = sum(px(i)), i=0 hasta k
rx <- 0:10
px <- px_poisson_generica(1, 2)
list_px <- sapply(rx, px)
fx <- fx_generica(rx, list_px)
sapply(0:10, fx)

# fx(6) = 0.9954, mientras que fx(5) = 0.9834 => deberia tener al menos 6 productos


# 16) Ubica 1000 libros por día. Si mal ubica un libro con prob 0.001, independiente. Aproximar
#     X = "cantidad de libros mal ubicados por día"
#     X ~ Bi(1000, p=0.001), pero como n = 1000 > 100, p < 0.01 y np = 1 < 20 => X ~ P(1)
#     Rx = {0, 1, ..., 1000}
#     a) P(X >= 1) = 1 - P(X <= 0) = 1 - P(X = 0)
rx <- 0:1000
px <- px_poisson_generica(1000, 0.001)
p_a <- 1-px(0)
p_a

#     b) P(X = 3)
p_b <- px(3)
p_b

px_exacta <- px_binomial_generica(1000, 0.001)
p_b_exacto <- px_exacta(3)
p_b_exacto


# 17) Pescador paga 100$ por participar. Cantidad de peces por pescador ~ P(4). 50$ por pez, max x6
#     X = "peces obtenidos"
#     Rx = {0, 1, ..., inf}
#     Y = "peces cobrables"
#     Ry = {0, 1, ..., 6}
#     G = "ganancia neta de un pescador"
#     Rg = {-100, -50, 0, 50, 100, 150, 200}
#     a) pg(g)
ry <- 0:6
py <- px_poisson_generica(4, 1)
list_py <- sapply(ry, py)
# Pescar 6 o > 6 es lo mismo, pues se puede cobrar hasta 6 máximo
# P(X=6) = 1 - sum P(X = i) para i = 0 hasta 5
py_6 <- 1 - sum(list_py[1:6])
probs <- c(list_py[1:6], py_6)
sum(probs)
rg <- c(-100, -50, 0, 50, 100, 150, 200)
rg
probs

#     b) Cuánto dinero espera ganar cada participante
ey <- ex_generica(ry, probs)
ey

eg <- 50*ey - 100
eg


# 18)
n <- 5
p <- 1/3
x <- 0:5
px <- dbinom(x, n, p)
fx <- pbinom(x, n, p)

px_2 <- px_binomial_generica(n, p)
list_px_2 <- sapply(x, px_2)
fx_2 <- fx_generica(x, list_px_2)
list_fx_2 <- sapply(x, fx_2)
# plot(x, px)
# plot(x, list_px_2)
# plot(x, fx, type = "s")
# plot(x, list_fx_2, type = "s")

# 19) Xt = "Tareas que llegan a un servidor por minuto", Xt ~ P(5) 
#     a) Prob de P(Xt < 5) en 30 segundos
# Xt2 = "Tareas que llegan a un server por cada 30 segs"
# Xt2 ~ P(2.5)
# P(Xt2 < 5) = P(Xt2 <= 4) = Fxt2(4)
# Rxt2 = {0, 1, ..., inf}
rx <- 0:10
px <- px_poisson_generica(2.5, 1)
list_px <- sapply(rx, px)
fx <- fx_generica(rx, list_px)
p_5 <- fx(4)
p_5

# Xt30 = "Tareas que llegan a un server cada media hora"
# Xt30 ~ P(150)
ex <- ex_poisson(150,1)
ex

# 20) X = "numero de veces que se bloquea una red por semana"
#     X ~ P(2)
#     a) En 2 semanas no se bloquee
# X2 = "numero de veces que se bloquea cada dos semanas"
# X2 ~ P(4)
# Rx = {0, 1, ..., inf}
px <- px_poisson_generica(4, 1)
px_0 <- px(0)
px_0

#     b) Periodo de 4 semanas, una semana en la que no se bloqueó
# Y = "cantidad de semanas en las que hubo bloqueos" cada 4 semanas
# Ry = {0, 1, 2, 3, 4}
# P(Y = 3) = P(X>0)^3 * P(X=0)
# P(X > 0) = 1 - P(X <= 0) = 1 - P(X = 0)
px <- px_poisson_generica(2, 1)
px_0 <- px(0)
py_1 <- ((1-px_0)^3) * px_0 * 4
py_1

# o con binomial,
# Y2 = "cantidad de semanas en las que hubo bloqueos"
# Y2 ~ Bi(4, P(X > 0))
# Ry2 = {0, 1, 2, 3, 4}
# P(X = 3)
n <- 4
rx <- 0:n
# P(X > 0) = 1 - P(X = 0)
p_exito <- 1-px_0
py2 <- px_binomial_generica(n, p_exito)
py2_3 <- py2(3)
py2_3


# Ejercicios adicionales
# 1)
n <- 15
rx <- 0:n
p_exito <- 1/5
px <- px_binomial_generica(n, p_exito)
list_px <- sapply(rx, px)
fx <- fx_generica(rx, list_px)
1-fx(5)

ex <- ex_generica(rx, list_px)
ex

plot(rx, list_px)
points(ex,0, col="red", pch=16)
list_px



# 2)
n <- 10
r_x <- 0:2
px <- dhyper(r_x, 1000, 29000, n)
fx <- fx_generica(r_x, px)

#       ii) P(X=6)?
px[3]


ex <- ex_generica(0:10, sapply(0:10, px_hiper_generica(45, 30000, 3600)))
ex
#   b) 500 inds, 350 no en infracción.
# para aproximar una hiper usamos la V.A binomial
r_x <- 0:n
p_exito <- 1/30
px_2 <- px_binomial_generica(n, p_exito)
list_px_2 <- sapply(r_x, px_2)
fx_2 <- fx_generica(r_x, list_px_2)
sapply(r_x, fx_2)

px_2(2)