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)) } f_empirica_t <- function(cx, t) { # mean(cx <= t) length(cx[cx <= t]) / length(cx) } # Intervalos de confianza # Longitud long_n <- function (zalpha, sgc, n) { zalpha*sqrt(sgc/n)*2 } obtener_repeticiones <- function(zalpha, sgc, longitud) { sgc/(longitud/(2*zalpha))^2 } IC_varianza_conocida <- function(alpha, sgc, n, un) { zalpha <- qnorm(alpha/2) raiz <- sqrt(sgc/n) valor <- (zalpha*raiz) cat("zalpha: ", zalpha, ", raiz: ", raiz, "\n") cat("IC: (", un+valor, ", ", un-valor, ")\n") } IC_varianza_desconocida <- function(alpha, Sc, n, un) { zalpha <- qt(alpha/2, n-1) raiz <- sqrt(Sc/n) valor <- (zalpha*raiz) cat("zalpha: ", zalpha, ", raiz: ", raiz, "\n") cat("IC: (", un+valor, ", ", un-valor, ")\n") } IC_media_conocida <- function(alpha, datos, u) { n <- length(datos) suma <- sum((datos - u)^2) za <- qchisq(alpha/2, n, lower.tail = FALSE) zb <- qchisq(1-(alpha/2), n, lower.tail = FALSE) cat("n: ", n, "suma: ", suma, "za: ", za, ", zb: ", zb, "\n") cat("IC: (", suma/za, ", ", suma/zb, ")\n") } IC_media_desconocida <- function(alpha, datos) { n <- length(datos)-1 sc <- var(datos) print(sum((datos-mean(datos))^2)/24) za <- qchisq(alpha/2, n, lower.tail = FALSE) zb <- qchisq(1-(alpha/2), n, lower.tail = FALSE) cat("n: ", n, "sc: ", sc, "za: ", za, ", zb: ", zb, "\n") cat("IC: (", sc*n/za, ", ", sc*n/zb, ")\n") } IC_exponencial <- function(alpha, datos) { n <- length(datos) suma <- sum(datos) za <- qchisq(1-(alpha/2), 2*n, lower.tail = FALSE) zb <- qchisq(alpha/2, 2*n, lower.tail = FALSE) cat("n: ", n, "suma: ", suma, "za: ", za, ", zb: ", zb, "\n") cat("IC: (", za/(2*suma), ", ", zb/(2*suma), ")\n") } IC_asintotico_exponencial <- function(alpha, datos) { n <- length(datos) promedio <- mean(datos) za <- qnorm(alpha/2) cat("n: ", n, "promedio: ", promedio, "za: ", za, ", zb: ", -za, "\n") cat("IC: (", 1/(promedio + (-za*promedio/sqrt(n))), ", ", 1/(promedio - (-za*promedio/sqrt(n))), ")\n") } IC_asintotico_binomial <- function(alpha, promedio, n) { za <- qnorm(alpha/2) cat("n: ", n, "promedio: ", promedio, "za: ", za, ", zb: ", -za, "\n") cat("IC: (", promedio + (za*sqrt(promedio*(1-promedio)/n)), ", ", promedio - (za*sqrt(promedio*(1-promedio)/n)), ")\n") } # Tests region_rechazo <- function(uh0, vh0, promedio, n, alpha, comparar) { zalpha <- qnorm(1-alpha) valor <- (promedio - uh0) / sqrt(vh0/n) cat("valor: ", valor, ", zalpha: ", zalpha, "\n") if (comparar(valor, zalpha)) { print("Rechaza H0!") } else { print("No rechaza H0") } } region_rechazo_t <- function(uh0, sc, promedio, n, alpha, comparar) { zalpha <- qt(1-alpha, n-1) valor <- (promedio - uh0) / sqrt(sc/n) cat("valor: ", valor, ", zalpha: ", zalpha, "\n") if (comparar(valor, zalpha)) { print("Rechaza H0!") } else { print("No rechaza H0") } } region_rechazo_chi <- function(sc, vh0, n, alpha, comparar) { zalpha <- qchisq(1-alpha, n-1) valor <- (n-1)*sc/vh0 cat("valor: ", valor, ", zalpha: ", zalpha, "\n") if (comparar(valor, zalpha)) { print("Rechaza H0!") } else { print("No rechaza H0") } } # R = {(promedio - u) / sqrt(vh0/n) >= (uh0 - u) / sqrt(vh0/n) + zalpha} # (promedio - u) / sqrt(vh0/n) ~ N(0,1) # => pi(u) = 1 - pnorm((uh0 - u)/sqrt(vh0/n) + zalpha) # P(EI) = calcular_f_potencia(...) # P(EII) = 1 - calcular_f_potencia(...) calcular_f_potencia_mayor <- function(uh0, vh0, uh1, n, alpha) { 1 - pnorm((uh0-uh1)/sqrt(vh0/n) + qnorm(1-alpha)) } # R = {(n-1)*s^2 / sg0^2 <= qchisq(alpha, n-1)}, multiplica por sg0^2 / sg1^2 # R = {(n-1)*s^2 / sg1^2 <= sg0^2*qchisq(alpha, n-1) / sg1^2} # (n-1)*s^2 / sg1^2 ~ chisq(n-1) # => pi(u) = pchisq(sg0^2*qchisq(alpha, n-1) / sg1^2) calcular_f_potencia_chisq_menor <- function(sg0c, sg1c, n, alpha) { pchisq(sg0c*qchisq(alpha, n-1) / sg1c, n-1) } # pnorm((uh0-uh1)/sqrt(vh0/n) + qnorm(1-alpha)) <= valor # (uh0-uh1) / sqrt(vh0/n) + qnorm(1-alpha) <= qnorm(valor) # (uh0-uh1) / sqrt(vh0/n) <= qnorm(valor) - qnorm(1-alpha) # (uh0-uh1) <= (qnorm(valor) - qnorm(1-alpha)) * sqrt(vh0/n) # (uh0-uh1)/(qnorm(valor) - qnorm(1-alpha)) <= sqrt(vh0/n) # ((uh0-uh1) / (qnorm(valor) - qnorm(1-alpha)))^2 <= vh0/n # |n| >= vh0 / ((uh0-uh1) / (qnorm(valor) - qnorm(1-alpha)))^2 calcular_n_error_tipo_2 <- function(valor, uh0, vh0, uh1, alpha) { vh0 / ((uh0-uh1)/(qnorm(valor)-qnorm(1-alpha)))^2 } #################################################################################################### #################################################################################################### # Práctica 9 # 1) Rendimiento por ha de soja ~ N(37, 25). Un vendedor promete mayor rendimiento. # Cultiva 10 parcelas de 1 ha, aplica T.H. de nivel 0.05, suponiendo que el nuevo rendimiento por # ha ~ N(u, 25) # a) H0: u <= 37, H1: u > 37 (mayor rendimiento) # Xi ~ N(u, 25) => promedio bajo H0 ~ N(37, 25/n) # (promedio - 37) / sqrt(25/n) ~ N(0,1) => R={(promedio - 37) / sqrt(25/n) >= qnorm(1-0.05)} # => R={promedio >= qnorm(1-0.05)*sqrt(25/n) + 37} # b) rendimientos <- c(37, 39.50, 41.70, 42, 40, 41.25, 43, 44.05, 38, 38.50) mean(rendimientos) region_rechazo(37, 25, mean(rendimientos), length(rendimientos), 0.05, function(s, t) { s > t }) # Con los valores observados, rechazaríamos H0, lo que nos indica que la nueva soja efecftivamente # tiene mayor rendimiento # Calculamos P(EII) = P(Rc) (cuando u = 40) = 1 - P(R) = 1 - pi(40) # pi(40) = P((promedio - 40) / sqrt(25/n) >= (37-40) / sqrt(25/n) + qnorm(1-0.05)) 1 - calcular_f_potencia_mayor(37, 25, 40, 10, 0.05) # P(EII) = 0.4003223 # d) Calcular la cantidad n de parcelas para P(EII) <= 0.05, u = 40 calcular_n_error_tipo_2(0.05, 37, 25, 40, 0.05) # n ~= 30.06159 => n >= 31 1 - calcular_f_potencia_mayor(37, 25, 40, 30, 0.05) 1 - calcular_f_potencia_mayor(37, 25, 40, 31, 0.05) # < 0.05 # 2) X = tension concreto, u > 300. Se realizan 15 mediciones i.i.d. u0 = 304, s = 10, alpha = 0.05 # a) Construir test # H0: u <= 300, H1: u > 300, # (promedio - 300)/sqrt(10^2/15) ~ t(n-1), R = {(promedio - 300) / sqrt(10^2/15) >= qt(1-0.05, 14)} region_rechazo_t(300, 10^2, 304, 15, 0.05, function(s,t) { s > t }) # No rechazaríamos H0 # b) Acotar el p-valor # Sabemos que el p-valor es > alpha dado que no rechazamos H0 # Buscamos el valor máximo que puede tener calculando la probabilidad P(Z >= zobs) # p-valor = P(Z >= zobs) = P((promedio - 300) / sqrt(10^2/15) >= (304-300)/sqrt(100/15)) 1 - pt((304-300)/sqrt(100/15), 14) # 0.07182 # 3) X = tiempo de activacion, X ~ N(u, sg^2). # H0: sg >= 6 , H1: sg < 6 , n=11 datos_3 <- c(27, 41, 22, 27, 23, 35, 30, 24, 27, 28, 22) n <- length(datos_3) # a) region_rechazo_chi(var(datos_3), 6^2, n, 0.95, function(s, t) { s < t }) # Con estas mediciones, vemos que no convendría rechazar H0 en favor de H1 (el desvío del tiempo # de activación no sería menor que 6 segundos según lo observado) # b) p-valor, debe ser mayor a 0.05 pues no podemos rechazar H0, y vale: pchisq(9.378788, 10) # 0.5034284 # c) pi(sg^2), es creciente? cuanto vale en sg = 6 ? calcular_f_potencia_chisq_menor(6^2, 6^2, 11, 0.05) # pi(u=u0) = alpha datos_3c <- seq(1, 7, by=0.5) # plot(datos_3c, sapply(datos_3c, function(t) {calcular_f_potencia_chisq_menor(6^2,t^2,11,0.05)} )) # abline(v=6) # Es decreciente # 4) 1/10 prefiere A. 200 fumadores, muestra que 26 prefieren la marca A # H0: p = 1/10, H1: p > 1/10 # X = "cantidad de fumadores que prefieren A", X ~ Bi(200, p) # Por TCL, sabemos que: (promedio - p) / sqrt(p*(1-p) / n) ~ N(0,1) # Bajo H0, el estadistico: (promedio - 1/10) / sqrt(1/10*(1-1/10) / n) ~ N(0,1) # Luego, R = {(promedio - 1/10) / sqrt(1/10*(1-1/10) / n) >= zalpha} region_rechazo(1/10, (1/10*9/10), 26/200, 200, 0.05, function(s, t) { s > t }) # No hay suficientes datos que indiquen que debamos rechazar H0 en preferencia por H1 # Es decir, no hay suficiente información para decir que con la publicidad se obutiveron seguidores # b) p-valor = P(Z >= zobs) = P(Z >= 1.414214) = 1 - P(Z < 1.414214) = 1 - F(1.414214) 1 - pnorm(1.414214) # 0.07864954 # c) "decir que no fue efectiva (aceptar H0) cuando la preferencia subió a 0.15 (H1 es valida)" # Es un error de tipo II => P(EII) = 1 - pi(0.15) # R = {promedio >= zalpha * sqrt(p0*(1-p0) / n) + p0} # R = {(promedio-p1) / sqrt(p1*(1-p1)/n) >= (zalpha*sqrt(p0*(1-p0)/n)+(p0-p1)) / sqrt(p1*(1-p1)/n)} zalpha <- 1.644854 p0 <- 0.1 n <- 200 p1 <- 0.15 pnorm( (zalpha*sqrt(p0*(1-p0)/n)+(p0 - p1)) / sqrt(p1*(1-p1)/n) ) # 0.2748061 # e) # 0.05 >= pnorm( (zalpha*sqrt(p0*(1-p0)/n)+(p0 - p1)) / sqrt(p1*(1-p1)/n) ) (((zalpha*sqrt(p0*(1-p0))) / sqrt(p1*(1-p1)) - qnorm(0.05)) * sqrt(p1*(1-p1)) / -(p0 - p1))^2 # n ~= 467.2397, n >= 468 n <- 467 pnorm( (zalpha*sqrt(p0*(1-p0)/n)+(p0 - p1)) / sqrt(p1*(1-p1)/n) ) # 0.05008014 n <- 468 pnorm( (zalpha*sqrt(p0*(1-p0)/n)+(p0 - p1)) / sqrt(p1*(1-p1)/n) ) # 0.04974664 # f) # Verdadero # Verdadero # Verdadero # Falso, es decreciente # OJO! Este ejercicio puede estar mal resuelto, nunca lo pude validar # 5) Di = "duracion de lamparita", Di ~ E(lambda), 1/lambda > 50 # n = 40, promedio = 53 # a) # H0: 1/lambda <= 50, lambda >= 1/50 , H1: 1/lambda > 50, lambda < 1/50 , alpha = 0.05 # 2*lambda*promedio ~ chisq(2n) # R = {2*lambda*promedio >= qchisq(1-alpha, 2n)} lambda <- 1/50 promedio <- 53 alpha <- 0.05 n <- 40 2*lambda*promedio >= qchisq(1-0.05, 2*n)/40 # No se rechaza # c) Podemos aproximar por una normal usando el TCL # (sum(Di)/n - 1/lambda)*sqrt(n)*lambda ~ N(0,1) region_rechazo(50, 50^2, 53, 40, 0.05, function(s, t) { s > t })