Punto 1

Graficamos la función de densidad de una variable aleatoria con distribución normal de media \(\mu=50\) y varianza \(\sigma^2=100\). Ésta es la distribución que siguen los puntajes de la prueba de Matemáticas del Icfes, para la cual se centran los puntajes en 50 y se toma una desviación estándar de 10 puntos.

mu = 50
sigma = 10

Como los puntajes toman valores de 0 a 100, asignamos este intervalo dentro del código de la gráfica:

curve(dnorm(x, mu, sigma), from = 0, to = 100,
      main = "Función de densidad de N(50, 10)")

Punto 2

Tomamos una muestra aleatoria de tamaño 100 mediante la función rnorm(), habiendo fijado previamente una semilla para reproducibilidad de los resultados del código:

set.seed(1305)
data <- rnorm(100, mu, sigma)
head(data) #Los primeros 6 valores generados
## [1] 53.58056 50.97156 56.01596 56.88797 51.89886 38.56211

Punto 3

Ejecutamos la estimación kernel de la densidad a partir de la muestra aleatoria generada en el apartado anterior. Esto, mediante la función density() con kernel = "gaussian":

d <- density(data, kernel = "gaussian")
plot(d, lwd = 2, main = "Estimación kernel con kernel gaussiano")

Comparamos visualmente la estimación con la curva teórica:

plot(d, lwd = 2, col = 2, main = "Estimación kernel con kernel gaussiano")
curve(dnorm(x, mu, sigma), from = 0, to = 100, col = 1, add = TRUE)
legend(3, .0375, lty = c(1, 1), legend = c("Densidad teórica", "Estimación kernel"), col = c(1, 2))

Punto 4

El ancho de banda \(h\) se modifica mediante el parámetro bw dentro de la función density():

plot(d, lwd = 2, col = 2, main = "Estimación kernel con kernel gaussiano y varios anchos de banda", xlab = "")
curve(dnorm(x, mu, sigma), from = 0, to = 100, col = 1, add = TRUE)
d2 <- density(data, kernel = "gaussian", bw = 2) #Estimación con bw=2
lines(d2, lwd = 2, col = 3)
d3 <- density(data, kernel = "gaussian", bw = 3) #Estimación con bw=3
lines(d3, lwd = 2, col = 4)
d4 <- density(data, kernel = "gaussian", bw = 5) #Estimación con bw=5
lines(d4, lwd = 2, col = 5)
legend(3, .0375, lty = c(1, 1, 1, 1, 1),
       legend = c("Densidad teórica", "bw por defecto", "bw=2", "bw=3", "bw=5"),
       col = c(1, 2, 3, 4, 5))

\(\sigma^2_K=1\) está implícito en el kernel gaussiano. Para cambiarlo, tendríamos que hacer una implementación personalizada del kernel.

Punto 5

Trabajaremos con el ancho de banda óptimo:

\[ h=1.059\,\widehat{\sigma}\,n^{-1/5}=1.059\,\widehat{\sigma}\,\cdot\,100^{-1/5}=0.4216\,\widehat{\sigma}. \]

Hacemos cinco simulaciones distintas con muestras aleatorias de tamaño \(n=100\) de la distribución \(\mathrm{N}(50, 10)\):

set.seed(13051)
data1 <- rnorm(100, mu, sigma)
set.seed(13052)
data2 <- rnorm(100, mu, sigma)
set.seed(13053)
data3 <- rnorm(100, mu, sigma)
set.seed(13054)
data4 <- rnorm(100, mu, sigma)
set.seed(13055)
data5 <- rnorm(100, mu, sigma)

E introducimos cada muestra dentro de density(), tomando como ancho de banda bw = 0.4216*sd(data), donde sd(data) calcula la desviación estándar empírica de cada una de las cinco muestras. Las curvas estimadas se sintetizan en el siguiente gráfico:

dbw1 <- density(data1, kernel = "gaussian", bw = 0.4216*sd(data1)) #Estimación con data1
plot(dbw1, lwd = 2, col = 2, main = "Estimación kernel con kernel gaussiano y varias muestras", xlab = "", ylim = c(0, .045))
curve(dnorm(x, mu, sigma), from = 0, to = 100, col = 1, add = TRUE)
dbw2 <- density(data2, kernel = "gaussian", bw = 0.4216*sd(data2)) #Estimación con data2
lines(dbw2, lwd = 2, col = 3)
dbw3 <- density(data3, kernel = "gaussian", bw = 0.4216*sd(data3)) #Estimación con data3
lines(dbw3, lwd = 2, col = 4)
dbw4 <- density(data4, kernel = "gaussian", bw = 0.4216*sd(data4)) #Estimación con data4
lines(dbw4, lwd = 2, col = 5)
dbw5 <- density(data5, kernel = "gaussian", bw = 0.4216*sd(data5)) #Estimación con data5
lines(dbw5, lwd = 2, col = 6)
legend(13, .0455, lty = c(1, 1, 1, 1, 1, 1),
       legend = c("Densidad teórica", "con muestra data1", "con muestra data2", "con muestra data3", "con muestra data4", "con muestra data5"),
       col = c(1, 2, 3, 4, 5, 6))

Nótese que \(\widehat{f}_K(x)\) se ajusta considerablemente bien a la curva teórica con cada muestra aleatoria de tamaño 100 que se genera.

Punto 6

set.seed(13056)
secuencia <- seq(-4, 4, length.out = 100)
f_hat <- matrix(0, nrow = 100, ncol = length(secuencia))

for (b in 1:100) {
  muestrab <- rnorm(100, mean = mu, sd = sigma)
  densidadb <- density(muestrab, kernel = "gaussian", bw = 5, from = -4, to = 4, n = length(secuencia))
  f_hat[b, ] <- densidadb$y
}

media <- colMeans(f_hat)
plot(secuencia, media, type = "l", col = "blue", lwd = 2,
     main = "Estimaciones Kernel y su Media", xlab = "x", ylab = "Densidad")
for (b in 1:100) {
  lines(secuencia, f_hat[b, ], col = rgb(1, 0, 0, alpha = 0.1))
}
lines(secuencia, media, col = "blue", lwd = 2)
lines(secuencia, dnorm(secuencia, mu, sigma), col = "black", lwd = 2, lty = 2)
legend("topright", legend = c("Media estimaciones", "Verdadera f(x)"),
       col = c("blue", "black"), lwd = 2, lty = c(1, 2))

Punto 7

sesgado <- media - dnorm(secuencia, mu, sigma)

plot(secuencia, sesgado, type = "l", col = "purple", lwd = 2,
     main = "Sesgo Estimado de la Densidad Kernel", ylab = "Sesgo", xlab = "x")
abline(h = 0, lty = 2)