prob_within_epsilon <- function(estimator_values, theta, eps){ mean(abs(estimator_values - theta) <= eps) }
ECM <- function(estimator_values, theta){ mean( (estimator_values - theta)^2 ) }
set.seed(123) mu <- 10 sigma <- 3
n_grid <- c(5,10,20,50,100,200,500,1000) # tamaños de muestra R <- 5000 # réplicas Monte Carlo
var_xbar <- var_x1 <- prob_eps_xbar <- prob_eps_x1 <- ecm_xbar <- ecm_x1 <- numeric(length(n_grid))
epsilon <- 0.3 # tolerancia
for(k in seq_along(n_grid)){ n <- n_grid[k] xbar <- x1 <- numeric(R) for(r in 1:R){ x <- rnorm(n, mean=mu, sd=sigma) xbar[r] <- mean(x) x1[r] <- x[1] } var_xbar[k] <- var(xbar) var_x1[k] <- var(x1) prob_eps_xbar[k] <- prob_within_epsilon(xbar, mu, epsilon) prob_eps_x1[k] <- prob_within_epsilon(x1, mu, epsilon) ecm_xbar[k] <- ECM(xbar, mu) ecm_x1[k] <- ECM(x1, mu) }
plot(n_grid, var_xbar, type=“b”, pch=19, xlab=“n”, ylab=“Varianza”, main=“Normal: Varianzas de Xbar vs X1”) lines(n_grid, var_x1, type=“b”, pch=17) legend(“topright”, legend=c(“Var(Xbar)”,“Var(X1)”), pch=c(19,17))
plot(n_grid, prob_eps_xbar, type=“b”, pch=19, ylim=c(0,1), xlab=“n”, ylab=paste0(“P(|Tn - mu| <=”, epsilon, “)”), main=“Normal: Probabilidad de estar cerca de mu”) lines(n_grid, prob_eps_x1, type=“b”, pch=17) abline(h=1, lty=2, col=“gray”) legend(“bottomright”, legend=c(“Xbar”,“X1”), pch=c(19,17))
plot(n_grid, ecm_xbar, type=“b”, pch=19, xlab=“n”, ylab=“ECM”, main=“Normal: ECM de Xbar vs X1”) lines(n_grid, ecm_x1, type=“b”, pch=17) legend(“topright”, legend=c(“ECM(Xbar)”,“ECM(X1)”), pch=c(19,17))
set.seed(123) p <- 0.3 n_grid <- c(5,10,20,50,100,200,500,1000) R <- 5000 epsilon_b <- 0.05
var_xbar_b <- var_x1_b <- prob_eps_xbar_b <- prob_eps_x1_b <- ecm_xbar_b <- ecm_x1_b <- numeric(length(n_grid))
for(k in seq_along(n_grid)){ n <- n_grid[k] xbar <- x1 <- numeric(R) for(r in 1:R){ x <- rbinom(n, size=1, prob=p) xbar[r] <- mean(x) x1[r] <- x[1] } var_xbar_b[k] <- var(xbar) var_x1_b[k] <- var(x1) prob_eps_xbar_b[k] <- prob_within_epsilon(xbar, p, epsilon_b) prob_eps_x1_b[k] <- prob_within_epsilon(x1, p, epsilon_b) ecm_xbar_b[k] <- ECM(xbar, p) ecm_x1_b[k] <- ECM(x1, p) }
plot(n_grid, var_xbar_b, type=“b”, pch=19, xlab=“n”, ylab=“Varianza”, main=“Bernoulli: Varianzas de Xbar vs X1”) lines(n_grid, var_x1_b, type=“b”, pch=17) legend(“topright”, legend=c(“Var(Xbar)”,“Var(X1)”), pch=c(19,17))
plot(n_grid, prob_eps_xbar_b, type=“b”, pch=19, ylim=c(0,1), xlab=“n”, ylab=paste0(“P(|Tn - p| <=”, epsilon_b, “)”), main=“Bernoulli: Probabilidad de estar cerca de p”) lines(n_grid, prob_eps_x1_b, type=“b”, pch=17) abline(h=1, lty=2, col=“gray”) legend(“bottomright”, legend=c(“Xbar”,“X1”), pch=c(19,17))
plot(n_grid, ecm_xbar_b, type=“b”, pch=19, xlab=“n”, ylab=“ECM”, main=“Bernoulli: ECM de Xbar vs X1”) lines(n_grid, ecm_x1_b, type=“b”, pch=17) legend(“topright”, legend=c(“ECM(Xbar)”,“ECM(X1)”), pch=c(19,17))
Ejercicio 1
set.seed(1)
mu <- 10 sigmas <- c(1,3,5) epsilons <- c(0.2, 0.5) n_grid <- c(5,10,20,50,100,200,500,1000) R <- 4000
prob_eps <- array(NA, dim = c(length(n_grid), length(sigmas), length(epsilons)), dimnames = list(paste0(“n=”, n_grid), paste0(“sigma=”, sigmas), paste0(“eps=”, epsilons))) var_xbar <- matrix(NA, nrow=length(n_grid), ncol=length(sigmas), dimnames=list(paste0(“n=”, n_grid), paste0(“sigma=”, sigmas))) ecm_xbar <- var_xbar
for(s in seq_along(sigmas)){ sigma <- sigmas[s] for(k in seq_along(n_grid)){ n <- n_grid[k] xbar <- replicate(R, mean(rnorm(n, mu, sigma))) var_xbar[k,s] <- var(xbar) ecm_xbar[k,s] <- mean((xbar - mu)^2) for(e in seq_along(epsilons)){ eps <- epsilons[e] prob_eps[k,s,e] <- mean(abs(xbar - mu) <= eps) } } }
par(mfrow=c(1,3)) matplot(n_grid, var_xbar, type=“b”, pch=1:3, lty=1, xlab=“n”, ylab=“Var(Xbar)”, main=“Normal: Var(Xbar) vs n”); legend(“topright”, col=1:3, pch=1:3, legend=colnames(var_xbar)) matplot(n_grid, ecm_xbar, type=“b”, pch=1:3, lty=1, xlab=“n”, ylab=“ECM(Xbar)”, main=“Normal: ECM(Xbar) vs n”); legend(“topright”, col=1:3, pch=1:3, legend=colnames(ecm_xbar))
plot(n_grid, prob_eps[, “sigma=3”, “eps=0.2”], type=“b”, ylim=c(0,1), main=“P(|Xbar-mu|<=eps) con n”, xlab=“n”, ylab=“prob”) lines(n_grid, prob_eps[, “sigma=1”, “eps=0.2”], type=“b”, pch=2, col=2) lines(n_grid, prob_eps[, “sigma=5”, “eps=0.2”], type=“b”, pch=3, col=3) lines(n_grid, prob_eps[, “sigma=3”, “eps=0.5”], type=“b”, pch=4, col=4) abline(h=1, lty=2) legend(“bottomright”, legend=c(“σ=3, ε=0.2”,“σ=1, ε=0.2”,“σ=5, ε=0.2”,“σ=3, ε=0.5”), col=c(1,2,3,4), pch=c(1,2,3,4)) par(mfrow=c(1,1))
Ejercicio 2
set.seed(2)
ps <- c(0.1, 0.5, 0.9) n_grid <- c(5,10,20,50,100,200,500,1000) R <- 4000 eps <- 0.05
var_b <- matrix(NA, nrow=length(n_grid), ncol=length(ps), dimnames=list(paste0(“n=”,n_grid), paste0(“p=”,ps))) prob_b <- var_b
for(j in seq_along(ps)){ p <- ps[j] for(k in seq_along(n_grid)){ n <- n_grid[k] xbar <- replicate(R, mean(rbinom(n, 1, p))) var_b[k,j] <- var(xbar) prob_b[k,j] <- mean(abs(xbar - p) <= eps) } }
par(mfrow=c(1,2)) matplot(n_grid, var_b, type=“b”, pch=1:3, lty=1, xlab=“n”, ylab=“Var(Xbar)”, main=“Bernoulli: Var(Xbar)=p(1-p)/n”) legend(“topright”, legend=colnames(var_b), pch=1:3, col=1:3)
matplot(n_grid, prob_b, type=“b”, pch=1:3, lty=1, ylim=c(0,1), xlab=“n”, ylab=paste0(“P(|Xbar-p|<=”, eps, “)”), main=“Bernoulli: Convergencia en prob.”) abline(h=1, lty=2) legend(“bottomright”, legend=colnames(prob_b), pch=1:3, col=1:3) par(mfrow=c(1,1))
Ejercicio 3
set.seed(3)
mu <- 5; sigma <- 2 n_grid <- c(5,10,20,50,100,200,500,1000) R <- 5000 eps <- 0.3
var_x1 <- numeric(length(n_grid)) prob_x1 <- numeric(length(n_grid))
for(k in seq_along(n_grid)){ n <- n_grid[k] # Ojo: Tn=X1, por lo que su distribución no depende de n x1_vals <- rnorm(R, mu, sigma) # basta 1 muestra por réplica var_x1[k] <- var(x1_vals) prob_x1[k] <- mean(abs(x1_vals - mu) <= eps) }
par(mfrow=c(1,2)) plot(n_grid, var_x1, type=“b”, pch=19, xlab=“n”, ylab=“Var(X1)”, main=“Var(X1) vs n (no depende de n)”) abline(h=sigma^2, lty=2, col=“gray”)
plot(n_grid, prob_x1, type=“b”, pch=19, ylim=c(0,1), xlab=“n”, ylab=paste0(“P(|X1-mu|<=”, eps, “)”), main=“Probabilidad de estar ‘cerca’ no ↑ a 1”) abline(h=1, lty=2, col=“gray”) par(mfrow=c(1,1))