# Parámetros
mu <- 5
sigma <- 10
n <- 1000
# Generar datos
set.seed(123)
datos <- rnorm(n, mu, sigma)
# Crear grid para densidad teórica
x <- seq(-20, 60, length.out = 1000)
f_teorica <- dnorm(x, mu, sigma)
# Calcular diferentes estimaciones kernel
f_default <- density(datos)
f_bw_small <- density(datos, bw = 0.2)
f_bw_large <- density(datos, bw = 0.8)
f_bw_optimal <- density(datos, bw = "SJ") # Método óptimo
# Graficar
par(mfrow = c(1, 2))
# Gráfico 1: Teórica vs Kernel por defecto
plot(x, f_teorica, type = "l", lwd = 2, col = "blue",
main = "Densidad Teórica vs Kernel",
xlab = "x", ylab = "Densidad")
lines(f_default, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Teórica", "Kernel (default)"),
col = c("blue", "red"), lwd = 2, lty = 1:2)
# Gráfico 2: Comparación de anchos de banda
plot(x, f_teorica, type = "l", lwd = 2, col = "blue",
main = "Efecto del Ancho de Banda",
xlab = "x", ylab = "Densidad")
lines(f_bw_small, col = "green", lwd = 2)
lines(f_bw_large, col = "purple", lwd = 2)
lines(f_bw_optimal, col = "orange", lwd = 2)
legend("topright",
legend = c("Teórica", "h=0.2", "h=0.8", "Óptimo (SJ)"),
col = c("blue", "green", "purple", "orange"), lwd = 2)
par(mfrow = c(1, 1))
# segundo
# Parámetros de la distribución normal
n <- 100 # Tamaño de muestra
# 1. Generar muestra aleatoria
set.seed(123) # Para reproducibilidad
muestra <- rnorm(n, mean = mu, sd = sigma)
# 2. Crear secuencia para la densidad teórica
x <- seq(-30, 60, length.out = 1000)
densidad_teorica <- dnorm(x, mean = mu, sd = sigma)
# 3. Estimación kernel con ancho de banda por defecto
densidad_kernel <- density(muestra)
# 4. Gráfico comparativo
plot(x, densidad_teorica, type = "l", lwd = 2, col = "blue",
main = paste("Densidad Teórica vs Estimación Kernel (n =", n, ")"),
xlab = "x", ylab = "Densidad",
ylim = c(0, max(densidad_teorica, densidad_kernel$y)))
# Superponer la estimación kernel
lines(densidad_kernel, col = "red", lwd = 2, lty = 2)
# Añadir puntos de la muestra
rug(muestra, col = "darkgray")
# Añadir leyenda
legend("topright",
legend = c("Densidad teórica N(0,1)", "Estimación kernel", "Muestra"),
col = c("blue", "red", "darkgray"),
lwd = c(2, 2, NA), lty = c(1, 2, NA), pch = c(NA, NA, "|"))
# TERCERO
# Parámetros de la distribución normal
mu <- 5 # Media
sigma <- 10 # Desviación estándar
n <- 100 # Tamaño de muestra
# 1. Generar muestra aleatoria
set.seed(123) # Para reproducibilidad
muestra <- rnorm(n, mean = mu, sd = sigma)
# 2. Crear secuencia para la densidad teórica
x <- seq(-20, 40, length.out = 1000)
densidad_teorica <- dnorm(x, mean = mu, sd = sigma)
# 3. Estimación kernel con kernel Gaussiano (σ² = 1)
# Esto se logra usando kernel = "gaussian" (opción por defecto)
densidad_kernel <- density(muestra, kernel = "gaussian", bw = "nrd0")
# 4. Gráfico comparativo
plot(x, densidad_teorica, type = "l", lwd = 2, col = "blue",
main = "Estimación Kernel Gaussiana (σ² = 1)",
xlab = "x", ylab = "Densidad",
ylim = c(0, max(densidad_teorica, densidad_kernel$y)))
# Superponer la estimación kernel
lines(densidad_kernel, col = "red", lwd = 2, lty = 2)
# Añadir puntos de la muestra
rug(muestra, col = "darkgray")
# Añadir leyenda
legend("topright",
legend = c("Densidad teórica N(0,1)", "Estimación kernel Gaussiano", "Muestra"),
col = c("blue", "red", "darkgray"),
lwd = c(2, 2, NA), lty = c(1, 2, NA), pch = c(NA, NA, "|"))
# Mostrar información sobre la estimación
cat("Parámetros de la estimación kernel:\n")
## Parámetros de la estimación kernel:
cat("Kernel usado:", densidad_kernel$kernel, "\n")
## Kernel usado:
cat("Ancho de banda (h) usado:", densidad_kernel$bw, "\n")
## Ancho de banda (h) usado: 3.170318
cat("Número de puntos de evaluación:", length(densidad_kernel$x), "\n")
## Número de puntos de evaluación: 512
# CUARTO PUNTO
# Generar datos de ejemplo: 1000 observaciones con media 5 y desviación 10
set.seed(123)
datos <- rnorm(1000, mean = 5, sd = 10)
# 1. Estimación de densidad por defecto
dens_default <- density(datos)
plot(dens_default, main = "Densidad por defecto", col = "blue", lwd = 2)
abline(v = 5, col = "red", lty = 2) # Línea en la media real
# 2. Modificar el ancho de banda (h) - más pequeño (menos suavizado)
dens_h_small <- density(datos, bw = 1) # bw = h
plot(dens_h_small, main = "Ancho de banda pequeño (h=1)", col = "darkgreen", lwd = 2)
abline(v = 5, col = "red", lty = 2)
# 3. Modificar el ancho de banda - más grande (más suavizado)
dens_h_large <- density(datos, bw = 5)
plot(dens_h_large, main = "Ancho de banda grande (h=5)", col = "purple", lwd = 2)
abline(v = 5, col = "red", lty = 2)
# 4. Cambiar el kernel (lo que afecta σₖ² indirectamente)
# Kernels disponibles: "gaussian", "epanechnikov", "rectangular",
# "triangular", "biweight", "cosine", "optcosine"
dens_epan <- density(datos, kernel = "epanechnikov")
plot(dens_epan, main = "Kernel Epanechnikov", col = "orange", lwd = 2)
abline(v = 5, col = "red", lty = 2)
# 5. Selección automática de ancho de banda (métodos alternativos)
dens_bw_nrd0 <- density(datos, bw = "nrd0") # Método por defecto
dens_bw_SJ <- density(datos, bw = "SJ") # Método Sheather & Jones
par(mfrow = c(1, 2))
plot(dens_bw_nrd0, main = "Método NRD0", col = "blue", lwd = 2)
abline(v = 5, col = "red", lty = 2)
plot(dens_bw_SJ, main = "Método SJ", col = "darkred", lwd = 2)
abline(v = 5, col = "red", lty = 2)
# QUINTO PUNTO
# -------------------------------
# 1. Configuración inicial
# -------------------------------
rm(list = ls())
set.seed(123)
# Paquetes necesarios
if (!require("ggplot2")) install.packages("ggplot2")
## Cargando paquete requerido: ggplot2
library(ggplot2)
# -------------------------------
# 2. Generación de datos simulados
# -------------------------------
n <- 1000 # Tamaño de muestra
datos <- c(rnorm(n/2, mean = -3, sd = 1.5), rnorm(n/2, mean = 3, sd = 1))
# Densidad verdadera (mezcla de normales)
true_dens <- function(x) {
0.5 * dnorm(x, -3, 1.5) + 0.5 * dnorm(x, 3, 1)
}
# -------------------------------
# 3. Configuración de kernels y anchos de banda
# -------------------------------
kernels <- c("gaussian", "epanechnikov", "rectangular") # Kernels disponibles en R
bws <- c(0.2, 0.5, 1, 2) # Valores de ancho de banda (h) a evaluar
# Varianzas teóricas de los kernels (σ²_K)
kernel_vars <- c(
"gaussian" = 1, # σ²_K = 1
"epanechnikov" = 0.2, # σ²_K ≈ 0.2
"rectangular" = 1/3 # σ²_K ≈ 0.333
)
# -------------------------------
# 4. Simulación y gráficos de densidad estimada
# -------------------------------
par(mfrow = c(length(kernels), length(bws)), mar = c(3, 3, 2, 1))
for (k in kernels) {
for (h in bws) {
dens_est <- density(datos, bw = h, kernel = k)
plot(dens_est,
main = paste0("Kernel: ", k, "\nh = ", h, ", σ²_K = ", round(kernel_vars[k], 3)),
col = "blue", lwd = 2, xlim = c(-8, 8), ylim = c(0, 0.25),
xlab = "", ylab = "")
curve(true_dens(x), add = TRUE, col = "red", lty = 2, lwd = 1.5)
legend("topright",
legend = c("Estimado", "Verdadero"),
col = c("blue", "red"), lty = c(1, 2), cex = 0.8)
}
}
# -------------------------------
# 5. Evaluación cuantitativa (MISE aproximado)
# -------------------------------
# Función para calcular el MISE aproximado
compute_mise <- function(dens_est, true_dens) {
x_points <- dens_est$x
delta <- diff(x_points)[1]
est_dens <- dens_est$y
true_dens_values <- true_dens(x_points)
mise <- sum((est_dens - true_dens_values)^2) * delta
return(mise)
}
# Data frame para almacenar resultados
results <- expand.grid(kernel = kernels, bw = bws)
results$mise <- NA
# Calcular MISE para cada combinación
for (i in 1:nrow(results)) {
k <- as.character(results$kernel[i])
h <- results$bw[i]
dens_est <- density(datos, bw = h, kernel = k)
results$mise[i] <- compute_mise(dens_est, true_dens)
}
# Mostrar resultados en una tabla ordenada
results$kernel_var <- kernel_vars[results$kernel]
results <- results[order(results$kernel, results$bw), ]
print(results)
## kernel bw mise kernel_var
## 1 gaussian 0.2 0.0009318167 1.0000000
## 4 gaussian 0.5 0.0008555817 1.0000000
## 7 gaussian 1.0 0.0067483801 1.0000000
## 10 gaussian 2.0 0.0274451103 1.0000000
## 2 epanechnikov 0.2 0.0008764302 0.2000000
## 5 epanechnikov 0.5 0.0008795729 0.2000000
## 8 epanechnikov 1.0 0.0079821451 0.2000000
## 11 epanechnikov 2.0 0.0348822692 0.2000000
## 3 rectangular 0.2 0.0009749357 0.3333333
## 6 rectangular 0.5 0.0009233362 0.3333333
## 9 rectangular 1.0 0.0088359545 0.3333333
## 12 rectangular 2.0 0.0411965390 0.3333333
# Gráfico de MISE vs. ancho de banda por kernel
ggplot(results, aes(x = bw, y = mise, color = kernel)) +
geom_line(size = 1) +
geom_point(size = 2) +
labs(title = "Error Cuadrático Medio Integrado (MISE) Aproximado",
x = "Ancho de banda (h)", y = "MISE",
color = "Kernel") +
theme_minimal() +
scale_color_manual(values = c("blue", "red", "green")) +
facet_wrap(~kernel, scales = "free_y")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# -------------------------------
# 6. Análisis de sensibilidad a h y σ²_K
# -------------------------------
# Gráfico 3D opcional (requiere el paquete 'plotly')
if (require("plotly")) {
plot_ly(results, x = ~bw, y = ~kernel_var, z = ~mise,
type = "scatter3d", mode = "markers",
color = ~kernel, size = 5) %>%
layout(scene = list(xaxis = list(title = "Ancho de banda (h)"),
yaxis = list(title = "σ²_K (varianza del kernel)"),
zaxis = list(title = "MISE")),
title = "Sensibilidad del MISE a h y σ²_K")
}
## Cargando paquete requerido: plotly
## Warning: package 'plotly' was built under R version 4.4.3
##
## Adjuntando el paquete: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
# SEXTO PUNTO
# -------------------------------
# 1. Configuración inicial
# -------------------------------
rm(list = ls())
set.seed(123)
# Paquetes necesarios
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("dplyr")) install.packages("dplyr")
## Cargando paquete requerido: dplyr
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(dplyr)
# -------------------------------
# 2. Generación de datos simulados
# -------------------------------
n <- 1000 # Tamaño de muestra por simulación
true_mean <- 5
true_sd <- 10
# Función de densidad verdadera (normal)
true_dens <- function(x) {
dnorm(x, mean = true_mean, sd = true_sd)
}
# -------------------------------
# 3. Parámetros del kernel
# -------------------------------
h <- 1.5 # Ancho de banda fijo
kernel <- "gaussian" # Kernel gaussiano
sigma2_K <- 1 # σ²_K para kernel gaussiano
# -------------------------------
# 4. Generar 100 estimaciones kernel
# -------------------------------
n_simulations <- 100
x_range <- seq(true_mean - 3*true_sd, true_mean + 3*true_sd, length.out = 200)
# Matriz para almacenar todas las estimaciones
all_estimates <- matrix(NA, nrow = length(x_range), ncol = n_simulations)
for (i in 1:n_simulations) {
# Generar nuevos datos para cada simulación
sim_data <- rnorm(n, mean = true_mean, sd = true_sd)
# Calcular densidad kernel
dens <- density(sim_data, bw = h, kernel = kernel, from = min(x_range), to = max(x_range), n = length(x_range))
# Almacenar estimación
all_estimates[, i] <- dens$y
}
# -------------------------------
# 5. Calcular la media de las estimaciones
# -------------------------------
mean_estimate <- rowMeans(all_estimates)
# -------------------------------
# 6. Crear dataframe para ggplot
# -------------------------------
plot_data <- data.frame(
x = x_range,
True = true_dens(x_range),
Mean_Estimate = mean_estimate
)
# Preparar datos para las 100 estimaciones (para geom_line)
sim_data_long <- data.frame(
x = rep(x_range, n_simulations),
y = as.vector(all_estimates),
Simulation = rep(1:n_simulations, each = length(x_range))
)
# -------------------------------
# 7. Crear la figura final
# -------------------------------
ggplot() +
# Líneas grises para las 100 estimaciones individuales
geom_line(data = sim_data_long, aes(x = x, y = y, group = Simulation),
color = "gray70", alpha = 0.3, linewidth = 0.3) +
# Línea roja para la densidad verdadera
geom_line(data = plot_data, aes(x = x, y = True),
color = "red", linewidth = 1.2, linetype = "dashed") +
# Línea azul para la media de las estimaciones
geom_line(data = plot_data, aes(x = x, y = Mean_Estimate),
color = "blue", linewidth = 1.2) +
# Añadir leyenda y formato
labs(title = paste("100 Estimaciones Kernel Gaussiano (h =", h, ", σ²_K =", sigma2_K, ")"),
subtitle = "Línea roja: Densidad verdadera | Línea azul: Media de 100 estimaciones",
x = "x", y = "Densidad") +
# Ajustar límites del eje y para mejor visualización
ylim(0, max(plot_data$True, plot_data$Mean_Estimate) * 1.1) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
# -------------------------------
# 8. Opcional: Mostrar algunas estadísticas
# -------------------------------
cat("Parámetros usados:\n")
## Parámetros usados:
cat(" - Número de simulaciones:", n_simulations, "\n")
## - Número de simulaciones: 100
cat(" - Tamaño de muestra por simulación:", n, "\n")
## - Tamaño de muestra por simulación: 1000
cat(" - Ancho de banda (h):", h, "\n")
## - Ancho de banda (h): 1.5
cat(" - Kernel:", kernel, "\n")
## - Kernel: gaussian
cat(" - σ²_K:", sigma2_K, "\n\n")
## - σ²_K: 1
# Calcular MISE entre la media de estimaciones y la densidad verdadera
delta_x <- diff(x_range)[1]
mise <- sum((mean_estimate - true_dens(x_range))^2) * delta_x
cat("Error Cuadrático Medio Integrado (MISE) entre la media de estimaciones y la densidad verdadera:", mise, "\n")
## Error Cuadrático Medio Integrado (MISE) entre la media de estimaciones y la densidad verdadera: 3.926955e-06