script_slides0_consit.R

#- EST0116 INFERENCIA II
#- Slides 0

# slide 23 ----------------------------------------------------------------

estu_sim1 <- function(eps, theta, K = 100) {
  ns <- seq(5, 10000, 1000)
  prob <- numeric()
  for (n in ns) {
    cont <- numeric()
    for (k in 1:K) {
      x <- rnorm(n, theta)
      cont <- append(cont, ifelse(abs(mean(x) - theta) < eps, 1, 0))
    }
    prob <- append(prob, mean(cont))
  }
  plot(ns, prob, type = "o", ylim = c(0, 1))
  abline(h = 1, col = "red")
}

par(mfrow = c(2, 2))
estu_sim1(eps = 0.5, theta = 0)
estu_sim1(eps = 0.1, theta = 0)
estu_sim1(eps = 0.05, theta = 0)
estu_sim1(eps = 0.01, theta = 0)

par(mfrow = c(2, 2))
estu_sim1(eps = 0.5, theta = -10)
estu_sim1(eps = 0.1, theta = -10)
estu_sim1(eps = 0.05, theta = -10)
estu_sim1(eps = 0.01, theta = -10)

par(mfrow = c(2, 2))
estu_sim1(eps = 0.5, theta = 10)
estu_sim1(eps = 0.1, theta = 10)
estu_sim1(eps = 0.05, theta = 10)
estu_sim1(eps = 0.01, theta = 10)

# slide 24 ----------------------------------------------------------------

estu_sim2 <- function(theta, K = 100, ...) {
  ns <- seq(5, 1000, 50)
  eqm <- numeric()
  for (n in ns) {
    esp <- numeric()
    for (k in 1:K) {
      x <- rnorm(n, theta)
      esp <- append(esp, (mean(x) - theta)^2)
    }
    eqm <- append(eqm, mean(esp))
  }
  plot(ns, eqm, type = "o", ...)
  abline(h = 0, col = "red")
}

par(mfrow = c(1, 1))
estu_sim2(theta = 0)
estu_sim2(theta = 10)
estu_sim2(theta = -10)


# slide 26 ----------------------------------------------------------------

estu_sim3 <- function(theta, K = 100, ...) {
  ns <- seq(5, 1000, 50)
  eam <- numeric()
  for (n in ns) {
    esp <- numeric()
    for (k in 1:K) {
      x <- rnorm(n, theta)
      esp <- append(esp, abs(mean(x) - theta))
    }
    eam <- append(eam, mean(esp))
  }
  plot(ns, eam, type = "o", ...)
  abline(h = 0, col = "red")
}

par(mfrow = c(2, 1))
estu_sim3(theta = 0, ylim = c(0, 0.5))
estu_sim3(theta = 10, ylim = c(0, 0.5))
estu_sim3(theta = -10, ylim = c(0, 0.5))


# slide 28 ----------------------------------------------------------------

estu_sim4 <- function(theta, K = 100, ...) {
  ns <- seq(5, 1000, 50)
  eqm <- numeric()
  for (n in ns) {
    esp <- numeric()
    for (k in 1:K) {
      x <- runif(n, 0, theta)
      esp <- append(esp, (max(x) - theta)^2)
    }
    eqm <- append(eqm, mean(esp))
  }
  plot(ns, eqm, type = "o", ...)
  abline(h = 0, col = "red")
}

par(mfrow = c(2, 2))
estu_sim4(theta = 1)
estu_sim4(theta = 5)
estu_sim4(theta = 10)
estu_sim4(theta = 50)


# slide 30 ----------------------------------------------------------------

estu_sim5 <- function(alpha, beta, K = 100, ...) {
  ns <- seq(5, 1000, 50)
  eqm.alpha <- numeric()
  eqm.beta <- numeric()
  for (n in ns) {
    esp.alpha <- numeric()
    esp.beta <- numeric()
    for (k in 1:K) {
      x <- rgamma(n, shape = alpha, rate = beta)
      alpha.est <- (n/(n - 1))*(mean(x))^2/var(x)
      beta.est <- (n/(n - 1))*mean(x)/var(x)
      esp.alpha <- append(esp.alpha, (alpha.est - alpha)^2)
      esp.beta <- append(esp.beta, (beta.est - beta)^2)
    }
    eqm.alpha <- append(eqm.alpha, mean(esp.alpha))
    eqm.beta <- append(eqm.beta, mean(esp.beta))
  }
  par(mfrow = c(1, 2))
  plot(ns, eqm.alpha, type = "o")
  abline(h = 0, col = "red")
  plot(ns, eqm.beta, type = "o")
  abline(h = 0, col = "red")
  par(mfrow = c(1, 1))
}

estu_sim5(alpha = 10, beta = 2)
estu_sim5(alpha = 10, beta = 10)
estu_sim5(alpha = 10, beta = 20)



# slide 31 ----------------------------------------------------------------

library(MASS)
estu_sim4 <- function(alpha, K = 100, ...) {
  ns <- seq(5, 1000, 50)
  eqm <- numeric()
  for (n in ns) {
    esp <- numeric()
    for (k in 1:K) {
      x <- rgamma(n, shape = alpha)
      est <- fitdistr(x, "gamma", rate = 1)
      emv <- as.numeric(est$estimate[1])
      esp <- append(esp, (emv - alpha)^2)
    }
    eqm <- append(eqm, mean(esp))
  }
  plot(ns, eqm, type = "o", ...)
  abline(h = 0, col = "red")
}

estu_sim4(alpha = 1, ylim = c(0, 1))
estu_sim4(alpha = 5, ylim = c(0, 1))
estu_sim4(alpha = 10, ylim = c(0, 1))