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

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

# Integrar
integrar <- function(fx, lower, upper) {
  return(integrate(fx, lower, upper)$value)
}

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


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

integrand <- function(x) {
  ifelse(x < -1 | x > 1, 0, 0.75*(1-(x^2)))
}
integrar(integrand, -Inf, 1)
integrar(integrand, 0, 1)
integrar(integrand, -0.5, 0.5)
integrar(integrand, -1, -0.25) + integrar(integrand, 0.25, 1)

x <- seq(-2, 2, by = 0.1)
# plot(x, sapply(x, integrand), type = "l", main = "Grafico fx densidad")
obtener_acumulada(integrand)(0)
Fx <- obtener_acumulada(integrand)
# plot(-2:2, sapply(-2:2,Fx), type = "l", main = "Grafico Fx acumulada")

# 2) Sea X una v.a continua con Fx
#    a) Valor de theta?
#       Como la Fx es continua en todo el rango -Inf, Inf, debe serlo en X=2, y como 
#       lim x-> 2 debe ser igual a 1 para ser continua, entonces 2^3/theta = 1
#       theta = 2^3 = 8
Fx <- function(x) {
  ifelse(x < 0, 0, ifelse(x>=2, 1, x^3/8))
}
x <- seq(-1, 3, by = 0.1)
# plot(x, sapply(x,Fx), type = "l", main = "Grafico Fx acumulada")

#    b) Calcular usando Fx(x),
#     i) P(X <= 1) = Fx(1)
p_x_i <- Fx(1)
p_x_i

#    ii) P(0.5 <= X <= 1) = Fx(1) - Fx(0.5), pues al ser Fx abs. continua, no importan los bordes
p_x_ii <- Fx(1) - Fx(0.5)
p_x_ii

#   iii) P(0.5 <= X <= 1 | X < 1) = P(int) / P(X < 1) = P(0.5 <= X < 1) / P(X < 1)
p_x_iii <- p_x_ii / p_x_i
p_x_iii

#    c) Mediana u de esta distribución. u = ?
# Al ser Fx continua, podemos pensar que 0.5 = Fx(u) => 0.5 = u^3/8
u <- (0.5*8)^(1/3)
u

#    d) Encontrar fx(x)
# fx <- deriv(~ ((x^3)/8), "x")
fx <- obtener_densidad_de_acumulada(expression((x^3)/8))
x <- seq(0,2, by=0.1)
# plot(x, eval(fx), type = "l")
# abline(v=u, col="red")

# 3) Sea U una v.a. ~ U(0, theta), para theta > 0
#    a) Hallar Fu(x)
# Fu(x) = {
#   - 0,        si x < 0
#   - x/theta,  si 0 <= x <= theta
#   - 1,        si theta < x
# }

#    b) P(1 <= U <= 3) = 0.5, qué valores puede tomar theta?
# P(1 <= U <= 3) = Fu(3) - Fu(1) = 3/theta - 1/theta = 2/theta = 0.5 => theta = 4

#    c) P(U^2 < 2) con theta = 4, P(U < sqrt(2)) = F(sqrt(2)) = sqrt(2) / 4
 
# 4) Y v.a con fy(t) = {
#   - t/25,         si 0 <= y < 5
#   - 2/5 - t/25,   si 5 <= y < 10
#   - 0,            otro caso
# }

#    a) Fy
fy_0 <- function(y) {
  ifelse(y >= 0 & y <= 5, y/25, 0)
}
Fy_0 <- integrar(fy_0, 0, 5)

fy_1 <- function(y) {
  ifelse(y >= 5 & y < 10, (2/5) - (y/25), 0)
}
Fy_1 <- integrar(fy_1, 5, 10)

Fy_3 <- 0

x0 <- seq(0, 5, by = 0.5)
x1 <- seq(5, 10, by = 0.5)
x <- c(x0, x1)
fy <- c(sapply(x0, fy_0), sapply(x1, fy_1))
plot(x, fy, type = "l")

# Fy(t) = {
#   - 0,                      si t < 0 o t >= 10
#   - t^2/50,                 si 0 <= t < 5
#   - (2t/5 - t^2/50) - 1,    si 5 <= t < 10
#   - 1,                      si t >= 10
# }

#   4.b) Calcular E(Y) y V(Y)
# Fy(u) = 0.5 = (2u/5 - u^2/50) - 1 => u = 5
# V(Y) = E(X^2) - E(X)^2 = E((X - E(X))^2) = integrate(f(x) * (x-ux)^2 dx, -Inf, Inf)

fy <- function(y) {ifelse(y >= 0 & y <= 5, y/25, 0)}
fy_2 <- function(y) {ifelse(y >= 5 & y < 10, (2/5) - (y/25), 0)}

ex1 <- ex_generica(fy, 0, 5)
ex2 <- ex_generica(fy_2, 5, 10)
ex <- ex1 + ex2
ex

vx_generica(fy, ex, 0, 5) + vx_generica(fy_2, ex, 5, 10)


# Dist normal
# qnorm(p, u, sigma) => devuelve para qué valor fx calcula p

# 5) n puntos al azar en [0, 1], independientemente, con dist uniforme.
#    X = "cantidad de puntos que caen en el intervalo [0, p] (0 < p < 1)"
#    Qué distribución tiene X?
# fx(x) = y(x) * I[0,p](x)
#   También es uniforme, pues y(x) es uniforme y su fy(y) = 1/y

frst <- Fx_normal_std_generica(1.25)
scnd <- Fx_normal_std_generica(-0.38)

# 6) Sea Z una v.a. ~ N(0,1)
#    a) P(0 <= Z <= 2) = Fx(2) - Fx(0)
a <- Fx_normal_std_generica(2) - Fx_normal_std_generica(0)
a

#    b) P(|Z| <= 2.5) = P(-2.5 <= Z <= 2.5) = Fx(2.5) - Fx(-2.5)
b <- Fx_normal_std_generica(2.5) - Fx_normal_std_generica(-2.5)
b

#    c) P(Z >= -1.37) = 1 - P(Z < -1.37) = 1 - Fx(-1.37) 
c <- 1 - Fx_normal_std_generica(-1.37)
c

#    d) c tal que P(Z < c) = 0.98
qnorm(0.98, 0, 1)

#    e) c tal que P(|Z| <= c) = 0.90 = P(-c <= Z <= c) = P(Z <= c) - (1-P(Z <= c))
# 1.90 = 2 * P(Z <= c) => P(Z <= c) = 0.95
qnorm(0.95, 0, 1)

#    f)
alphas <- c(0.1, 0.05, 0.025, 0.01)
sapply(alphas, function(x) { qnorm(1-x, 0, 1)})

# 7) X ~ N(5, 0.25)
Fx <- Fx_normal_generica_a_std(5, 0.25)
#   a) P(4.75 <= X <= 5.5)
Fx(5.5) - Fx(4.75)

#   b) P(|X| > 5.25) = 1 - P(|X| <= 5.25) = 1 - (P(-5.25 <= X <= 5.25))
#  P(-5.25 <= X <= 5.25) = Fx(5.25) - Fx(-5.25)
1 - Fx(5.25) - Fx(-5.25)

#   c) c tal que P(|X - 5| <= c) = 0.9 = P(-c + 5 <= X <= c + 5) = Fx(c+5) - Fx(-c + 5)
#      Fx(c+5) - (1 - Fx(c+5)) = 2*Fx(c+5) - 1 = 0.9
#      Fx(c+5) = 0.95
#      Z = (X - u) / sg => Z = ((c+5) - u) / sg => c = qnorm(Fx(...)) * sg + u - 5
#      Fz((c+5 - u) / sg) = 0.95
(qnorm(0.95)*0.5) + 5 - 5

arg_Fx_normal_a_std(5, 0.25)(0.95)-5

#   d) El 90-percentil de X => P(X <= c) = 0.9
arg_Fx_normal_a_std(5, 0.25)(0.9)

# 8) Índice cefálico I (anchura del cráneo como % de la longitud) es una v.a ~ N(u, sgc).
#    58% con I <= 75        => P(I <= 75) = Fi(75) = 0.58
#    38% con 75 < I <= 80   => P(75 < I <= 80) = Fi(80) - Fi(75) = 0.38 => Fi(80) = 0.96
#    4%  con I > 80         => P(I > 80) = 1 - P(I <= 80) = 1 - Fi(80) = 0.04
#    Fx(75) = Fz(75 - u / sg) = 0.58
a <- qnorm(0.58)
a
#    Fx(80) = Fz(80 - u / sg) = 0.96
b <- qnorm(0.96)
b

# {
#   (75 - u) / sg = 0.2018935 => u = 75 - 0.2018935*sg
#   (80 - u) / sg = 1.750686 => (5 + 0.2018935*sg) / sg => 5/sg + 0.2018935 = 1.750686
# }
sg <- 5/ (b-a)
u <- 75 - (a*sg)
sg
u

#    a) Hallar fx de I
fx <- fx_normal_generica(u, sg)
x <- seq(64.3, 84.3, by=0.1)
# plot(x, sapply(x, fx), type = "l", main = "Gráfico de fx densidad")
Fz <- Fx_normal_generica_a_std(u, sg^2)
# plot(x, sapply(x, Fz), type = "l", main = "Gráfico de Fx acumulada")

#    b) Calcular P(78 <= I <= 82)
Fz(82) - Fz(78)


# 9) IQR es la distancia entre Q3 y Q1, siendo Q3 el valor que obtiene el 75-percentil,
#    y Q1, el 25-percentil
#    IQR = Q3 - Q1 = qnorm(0.75) - qnorm(0.25)
dist_iqr <- qnorm(0.75) - qnorm(0.25)
dist_iqr

x <- seq(-5, 5, by = 0.1)
# plot(x, dnorm(x), type="l", main = "Distancia interquartil de N(0,1)")
# abline(v=-dist_iqr, col="red")
# abline(v=dist_iqr, col="red")

# 11) Red de computadoras, tiempo destintado a búsqueda tiene una fx,
#     fx(t) = c(100 - t)*I[0,100](t)
#     a) Hallar la constante c
# integrar(fx, -Inf, Inf) = 1 => la indicadora vale 1 entre 0 y 100, y 0 para el resto
# integrar(fx, 0, 100) => c*(100 - t) dt => c*integrar((100-t dt)) = 1
# 100t - t^2 / 2 = 1/c => 100*100 - 100^2 / 2 = 1/c
dividendo <- 100^2 - (100^2)/2
dividendo
c <- 1/dividendo
c

integrando <- function(t) {
  0.0002*(100 - t)
}

# plot(0:100, sapply(0:100, integrando), type="l")

# Comprobamos que la integral de fx de 0 a 100 es 1
integrar(integrando, 0, 100)

#     b) Clasificación de usuario segun porcentaje de tiempo: {
#          - tipo 1,    si T < 25%
#          - tipo 2,    si 25% <= T < 50%
#          - tipo 3,    si 50% <= T < 75%
#          - tipo 4,    si T > 75%
#        }
# U = "tipo de usuario", Ru={1,2,3,4}
# 


# 12) Diametro D (en dm) es una v.a con fd(x) = k*x*I(0,10)(x)
#     a) Hallar el valor de k
# integrar(k*x*I(0,10)(x), -Inf, Inf) = 1 => I vale 0 si no está en (0,10) =>
# integrar(k*x, 0, 10) => k*x dx => k(integrar(x, 0, 10)) => k*(x^2/2 | (10,0))
k <- 2/(10^2)
k

#     b) P(4 <= D <= 6) = integrate(k*x, 4, 6)
fd <- function(x) {
  k*x
}
integrar(fd, 4, 6)

#     c) P(4 <= D <= 6 | D > 5) = P(5 < D < 6) / P(D > 5)
# P(D > 5) = 1 - P(D <= 5) = 1 - integrar(k*x, 0, 5)
integrar(fd, 5, 6) / (1-integrar(fd, 0, 5))

#     d) A = "arboles que tienen diametro entre 4 y 6"
#        A ~ Bi(3, 0.11), Ra={0,1,2,3}
# P(A = 2) = 1 * (0.11)^2 * (1-0.11)
p_a_2 <- choose(3, 2) * (0.2)^2 * (1-0.2)
p_a_2


# 13) Variedades I y II de tabaco, P(I) = 0.35, P(II) = 0.65
# Para I, contenido de nicotina v.a ~ N(1.9, 0.16), y para II, ~ N(2.2, 0.09)
#     a) P(Nicotina >= 2.1) = P(Nicotina >= 2.1 | I)*P(I) + P(Nicotina >= 2.1 | II)*P(II)
# Para I: 
# P(Nic >= 2.1) = 1 - P(Nic < 2.1) => Z = (X - u) / sg => Fz(z) = P(Z <= z) = P(X-u/sg <= 2.1-u/sg)
# P(Nic >= 2.1) = 1 - P(Z <= 2.1-u /sg) de una N(0,1)

Fx <- Fx_normal_generica_a_std(1.9, 0.16)
a <- 1-Fx(2.1)

Fx2 <- Fx_normal_generica_a_std(2.2, 0.09)
b <- 1-Fx2(2.1)

p_nic_2.1 <- a*0.35 + b*0.65
p_nic_2.1

#     b) P(I | Nic >= 2.1) = P(Nic >= 2.1 | I) / P(nic2.1)
p_i_dado_nic_2.1 <- a*0.35 / p_nic_2.1 
p_i_dado_nic_2.1


# 14) X toma valores entre 0 y 1, fx(x) = 4x^3, si 0 < x < 1
#     a) Mediana de X
integrando <- function(x) { 4*(x^3) }
integrar(integrando, 0, 1) # Es funcion de densidad
# x^4 entre 0 y k me da 0.5 => k^4 = 0.5 => k = raiz_cuarta(0.5) ~= 0.8408964
0.8408964^4

#     b) Densidad de la variable Z = 1 - X, estrictamente decreciente
# z = 1-x = g(x) ,  g-1(z) = 1-z  ,  g-1(z)' = -1
# fz = fx(g-1(z))*g-1(z)' = -fx(g-1(z))*g-1(z)' = -fx(1-z)*-1 = fx(1-z)
# 0 < x < 1 => 0 < 1-z < 1 => 1 > z > 0
# fz(z) = 4*(1-z)^3 entre 0 < z < 1


# 15)
# a) Fx(x) = integrate(xe^-(x^2/2), 0, t) = 1-e^(x^2/2)
# b) P(T < 1) = Fx(1) = 1 - e^(-1/2) = 0.39346934
# c) T^2 = Z, en el mismo rango (0, Inf) es estrictamente creciente
# g(T) = T^2 => g-1(Z) = Z^(1/2)  => g-1(z)'=Z^(-1/2)*1/2 = 1/2*sqrt(Z)
# fz(z) = fx(g-1(z))*g-1(z)' = fx(sqrt(z)) * 1/2*sqrt(z)
# fz(z) = sqrt(z)*e^-(|z|/2) * 1/2*sqrt(z) = 1/2 * e^(-z/2)


# 16) V = "vida util en meses de un componente electrónico"
#     V ~ E(lambda) => ~ Gamma(1, lambda)
# P(V > 20) = 0.449 = 1 - P(V <= 20) => Fv(20) = 1 - 0.449 = 0.551 = 1 - e^(-lambda * 20)
# 0.449 = e^(-lambda * 20) => ln(0.449) / -20 = lambda = 0.04003661
#     a) E(V) = 1/lambda = 1/0.04003661 = 24.9771
#        Var(V) = 1/(lambda^2) = 623.8575062
x= seq(0,100,length.out = 100)
lambda = 0.04003661 # Desviacion estandar
dexp(x,rate=lambda)
plot(x,dexp(x,rate=lambda),type='l')
abline(v=24.9771)

plot(x, pexp(x, 0.04003661),type='l')
abline(h=0.5)

#     b) P(V > 10) = 1 - P(V <= 10) = 1 - Fv(10)
b <- 1-pexp(10, 0.04003661)
b

#     c) P(V > 30 | V > 20) = P(V > 10) = b
b


# 17) Xt ~ P(4*t). T="tiempo de espera hasta que llegue la primer tarea, en minutos"
# P(T <= 15) = F(t) = 1 - P(Xt2 = 0), donde Xt2 ~ P(1*t), "cantidad de tareas por cada 15 seg"
# P(T <= 15) = 1 - (e^(-1)*1 /1) 
a <- 1 - (exp(-1)*1 /1) 
a

pexp(1)

# P(T <= 15) = F(t) = 1 - P(Xt3 = 0), donde Xt3 ~ P(1/60*t), "cantidad de tareas por seg"
a <- 1 - exp(-4/60*15) 
a

pexp(15, 4/60)
pexp(0.25, 4)


# 18) Sistema de 5 * comps como el 14), conectados en serie
# Ai = "el i-ésimo componente dura por lo menos hasta t", con i={1,2,...,5}
# X = "momento en el cual el sistema se desconecta"
#     a) {X >= t} = A1t n A2t n ... n A5t = n(Ai) con i={1,2,...,5}
#     b) Como son independientes, seria una productoria, A1t * A2t * ...
# P(X >= t) = P(A1) * P(A2) * ... * P(A5)
# P(Ai) = P(Vi >= t) = 1 - P(Vi <= t) = 1 - (1 - e^(-0.04*t)) = e^(-0.04*t)
# => P(X >= t) = P(V1 >= t) * P(V2 >= t) ... = P(V1 >= t)^5, pues son iguales
# => P(X >= t) = (e^(-0.04*t))^5 = e^(-0.04*5*t) = e^(-0.2*t)
#     c) Fx(t) = {
#   0,               si t <= 0
#   1 - e^(-0.2*t), si t > 0
# }
#
# fx(t) = Fx(t)' = (1-e^(-0.2*t))' = e^(-0.2*t)*-0.2 => X ~ Exp(0.2)

# 19) 3 bloques, tiempo de compilacion en segundos T ~ Exp(1), independiente
#     a) X = "cantidad de bloques que tienen tiempo de comp. superior a 2 segundos"
#        X ~ Bi(3, P(B > 2)), donde P(B > 2) = exp(-2)
# P(X >= 1) = 1 - P(X < 1) = 1 - Fx(0) = 1 - P(X = 0)
p_exito <- exp(-2)
a <- 1 - (choose(3,0)*1*(1-p_exito)^3)
a

#     b) P(X = 2)
b <- (choose(3,2)*(p_exito^2)*(1-p_exito))
b



# 22) X ~ U(0, 1)
#     a) cX + d = Z  =>  fz(z) = fx(g-1(z))*g-1(z)'
# g(x) = cX + d  ,  g-1(z) = (z - d) / c  ,  g-1(z)' = 1/c
# fz(z) = fx((z - d) / c)) * z/c = 1/c * I[d, c+d](z), con d <= z <= c+d

# Por otro lado, si c < 0, la inversa de X da vuelta la desigualdad, g-1(z)' da negativa

# 25) Generar a partir de X ~ U(0,1), una variable aleatoria con distribución U(3,8)

# px(k) = 5^k * e^(-5) / k!