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 5 # 1) X1 ~ U(3,5), X2 ~ N(3,1/4), X3 ~ Gamma(6,3), X4 ~ U(2,4) en minutos, independientes. # Y = "tiempo que tarda un producto en pasar por toda la línea de producción" # Hallar E(Y) y V(Y) e_y <- ex_uniforme_generica(3, 5) + 3 + ex_gamma_generica(6,3) + ex_uniforme_generica(2, 4) e_y v_y <- vx_uniforme_generica(3, 5) + 1/4 + vx_gamma_generica(6,3) + vx_uniforme_generica(2, 4) v_y # 2) # a) Sea X una v.a., distribución desconocida. E(X) = 5, V(X) = 0.1 # P(4.5 < X < 5.5) ~= P(|X - 5| > 0.5) # P(|X - 5| > 0.5) = 1 - P(|X - 5| <= 0.5) = 1 - P(-0.5 <= X - 5 <= 0.5) = P(4.5 <= X <= 4.5) # Además, P(|X - 5| > 0.5) <= 0.1 / (0.5)^2 = 0.4 a <- 1 - calcular_chebyshev(0.1, 0.5) a # b) Sea X1, ..., Xn muestra aleatoria (iid) con E(X1) = 5, V(X1) = 0.1, y x̅ su promedio # Cota inferior cuando n = 10, P(4.5 < x̅ < 5.5 ~= P(x̅ - 5| > 0.5) # P(|x̅ - 5| > 0.5) <= 0.1 / 10*(0.5)^2 = 0.04 b <- 1 - calcular_chebyshev(0.1/10, 0.5) b # c) La cota tiende a 1 cuando n tiende a infinito # 3) X = "cantidad de aviones que aterrizan", X ~ P(10*t), t=20 minutos # Y ~ P(10*3*t) t = 60 min => 30 aviones cada hora # Determinar una cota inferior para que P(20 < Y < 40) # E(Y) = 30, y P(20 <= Y <= 40) = P(-10 <= Y - 30 <= 10) = P(|Y - 30| <= 10) = 1 - P(|Y - 30| > 10) # Por Chebyshev, P(|Y - 30| > 10) <= sg^2 / 10^2 , sg^2 = 30, luego P(|Y - 30| > 10) <= 0.3 # P(20 <= Y <= 40) = 1 - P(|Y - 30| > 10) >= 0.7 # 4) Por ley de los grandes numeros: # P(|fn - p| > 0.1) <= p*(1-p) / n*(0.1^2), y como p*(1-p) como mucho vale 0.25 cuando p=0.5, luego # P(|fn - p| > 0.1) <= 0.25 / n*0.01 = 25/n # Si tomamos el lim cuando n -> infinito, fn se acerca a p, luego mejoraríamos con más encuestas # b) 25/n = 0.1 luego n = 250 # 7) X = "peso de queso en kg", E(X) = 2, V(X) = 0.04 # a) n = 60, calcular de forma aproximada que 60 quesos pesen más de 122kg # P(sum(Qi) i=1 hasta 60 > 122) => Por TCL, (T60 - 60*2) / (sqrt(60*0.04)) -> N(0,1) # Luego P(sum(Qi) i=1 hasta 60 > 122) ~= 1 - P((T60 - 60*2) / (sqrt(60*0.04)) <= sqrt(15)/3) # => 1 - pnorm(sqrt(15)/3) 1-pnorm(sqrt(15)/3) # b) P(sum(Qi) i=1 hasta n >= 5000) >= 0.95 => Por TCL, (Tn - n*2) / (sqrt(n*0.04)) -> N(0,1) # P((Tn - n*2) / (sqrt(n*0.04)) >= 5000-(n*2) / (sqrt(n*0.04))) >= 0.95 => # 1 - P((Tn - n*2) / (sqrt(n*0.04)) < 5000-(n*2) / (sqrt(n*0.04))) >= 0.95 # 1 - FQ(...) = 0.95 => FQ(...) = 0.05 argumento <- qnorm(0.05) argumento # 5000 - n*2 / sqrt(n)*0.2 = argumento => argumento*sqrt(n) = 25000 - n*2/0.2 # => -1.644854*sqrt(n) + 10*n = 25000 # => n ~= 2508, luego con n = 2509 aproximadamente se debería poder satisfacer el pedido de 5000kg # 8) Sea X ~ Bi(100,0.8), E(X) = 80, V(X) = 16 # a) P(75 <= X <= 85) => por TCL, X ~ N(u, sg^2) = N(80, 16) # P(75 <= X <= 85) = P(-5 / 4 <= X - 80 / 4 <= 5 / 4) ~= F(5/4) - F(-5/4) aprox_a <- pnorm(5/4) - pnorm(-5/4) aprox_a real_a <- pbinom(85, 100, 0.8) - pbinom(75, 100, 0.8) aprox_a/real_a # b) P(X >= 80) = 1 - P(X <= 79) => por TCL, X ~ N(80, 16) # P(X <= 79) = pnorm() 1 - pnorm(79, 80, 4) 1-Fx_normal_generica_a_std(80, 16)(79) 1 - pbinom(79, 100, 0.8) # c) P(0.7 <= X/100 <= 0.8) = P(70 <= X <= 80) # 9) Diariamente A, B y C llenos de piedras. Xa, Xb, Xc v.a. independientes, donde # tienen esperanza 1.8, 3.8 y 4.1, y varianza 0.1, 0.18, 0.25, respectivamente # a) D = A+B+C => P(sum(Di) desde 1 hasta 360 > 3500) # E(D) = E(A)+E(B)+E(C) = 1.8+3.8+4.1 = 9.7 # V(D) = V(A)+V(B)+V(C) = 0.53 (pues son independientes) # P(sum(Di) > 3500) => Por TCL, Tn - 360*9.7 / sqrt(360*0.53) ~ N(0,1) # P((Tn - 3492)/(sqrt(360*0.53)) > (3500-3492)/(sqrt(360*0.53))) = # 1 - P((Tn - 3492)/(sqrt(360*0.53)) <= (3500-3492)/(sqrt(360*0.53))) # 1 - F((3500-3492)/(sqrt(360*0.53))) 1 - pnorm((3500-3492)/(sqrt(360*0.53))) # b) P(sum(Di) desde 1 hasta n > 3500) >= 0.9 # P((Tn - (n*9.7)) / (sqrt(n*0.53)) > (3500-(n*9.7)) / (sqrt(n*0.53))) >= 0.9 # 1 - P((Tn - (n*9.7)) / (sqrt(n*0.53)) <= (3500-(n*9.7)) / (sqrt(n*0.53))) >= 0.9 # P((Tn - (n*9.7)) / (sqrt(n*0.53)) <= (3500-(n*9.7)) / (sqrt(n*0.53))) <= 0.1 qnorm(0.1) # -1.281552 # (3500-(n*9.7)) / (sqrt(n*0.53)) <= -1.281552 # x ~= 362.656, luego como minimo x >= 363 1 - pnorm((3500-363*9.7)/(sqrt(363*0.53))) # P((Tn - (n*9.7)) > (3500-(n*9.7))) = P(|Tn - (n*9.7)| <= (3500-(n*9.7))) >= 0.1 # 0.1 = n*0.53/(3500-(n*9.7))^2 # Con chebyshev, x ~= 365.361, luego x >= 366 # 11) Probabilidad de ganar = 0.3 , se paga 1$ para jugar, si gana se reciben 5$ # a) Proba aproximada de que en 100 juegos, el jugador gane más de 90 pesos # P(sum(Ji) desde 1 hasta 100 > 90) # Ji = "cantidad de dinero intercambiado", J = 5*X - 1, X ~ Be(0.3) # E(J) = 5*E(X) - 1 = 5*0.3 - 1 = 0.5 # V(J) = 25*V(X) = 25*0.3*0.7 = 5.25 # P(sum(Ji) desde 1 hasta 100 > 90) => por TCL, (T100 - 100*0.5)/(sqrt(100*5.25)) ~ N(0,1) # P((T100 - 100*0.5)/(sqrt(100*5.25)) > (90-(100*0.5))/(sqrt(100*5.25))) = # 1 - P((T100 - 100*0.5)/(sqrt(100*5.25)) <= (90-(100*0.5))/(sqrt(100*5.25))) ~= # 1 - F((90-(100*0.5))/(sqrt(100*5.25))) 1 - pnorm((90-(100*0.5))/(sqrt(100*5.25))) # b) n = ? si P(sum(Ji) desde 1 hasta n > 90) >= 0.8 # P(sum(Ji) desde 1 hasta n > 90) => por TCL, (Tn - n*0.5)/(sqrt(n*5.25)) ~ N(0,1) # 1 - P((Tn - n*0.5)/(sqrt(n*5.25)) <= (90-(n*0.5))/(sqrt(n*5.25))) >= 0.8 # 0.2 >= P((Tn - n*0.5)/(sqrt(n*5.25)) <= (90-(n*0.5))/(sqrt(n*5.25))) qnorm(0.2) # -0.8416212 # (90-(n*0.5))/(sqrt(n*5.25)) >= -0.8416212 # n ~= 239.714 => al menos n = 240 1 - pnorm((90-(239.714*0.5))/(sqrt(239.714*5.25)))