Una muestra aleatoria de seis autos de un determinado modelo evidencia que cada uno de ellos consume las siguientes cantidades en kms por litro:
\[ 18.6, \quad 18.4, \quad 19.4, \quad 20.4, \quad 19.4, \quad 20.5 \]
Queremos determinar la probabilidad de que el consumo medio muestral sea menor que 17.6 km/l, suponiendo que la distribución de la población es normal con media 17.
Definimos los valores del problema:
# Datos de la muestra
datos <- c(18.6, 18.4, 19.4, 20.4, 19.4, 20.5)
n <- length(datos) # Tamaño de la muestra
# Parámetros poblacionales
mu <- 17 # Media poblacional
# Valor de referencia para la media muestral
x_ref <- 17.6
La media muestral \(\bar{X}\) se calcula como: \[ \bar{X} = \frac{\sum X_i}{n} \]
La desviación estándar muestral \(s\) es: \[ s = \sqrt{\frac{\sum (X_i - \bar{X})^2}{n-1}} \]
El error estándar de la media es: \[ \sigma_{\bar{X}} = \frac{s}{\sqrt{n}} \]
# Media muestral
x_barra <- mean(datos)
# Desviación estándar muestral
s <- sd(datos)
# Error estándar de la media
sigma_x <- s / sqrt(n)
# Mostrar resultados
x_barra; s; sigma_x
## [1] 19.45
## [1] 0.8757854
## [1] 0.3575379
Para convertir el valor de interés \(x_{ref} = 17.6\) a la distribución normal estándar, usamos la transformación: \[ Z = \frac{X - \mu}{\sigma_{\bar{X}}} \]
# Cálculo del valor Z
z <- (x_ref - mu) / sigma_x
# Mostrar valor Z
z
## [1] 1.678144
La probabilidad acumulada se obtiene con la función
pnorm
:
\[ P(\bar{X} < 17.6) = P(Z < z) \]
# Probabilidad acumulada en la distribución normal estándar
probabilidad <- pnorm(z)
# Mostrar resultado
probabilidad
## [1] 0.9533405
La probabilidad de que el consumo medio muestral de gasolina sea menor que 17.6 km/l es 95.33%.
Esto significa que, en 95.33% de las muestras de seis autos, la media del consumo estará por debajo de 17.6 km/l. Si esta probabilidad es baja, entonces un consumo tan bajo sería poco común en la población.
Para utilizar la aproximación normal en proporciones, debemos
verificar que:
\[
np > 5 \quad \text{y} \quad n(1 - p) > 5
\]
# Parámetros
p <- 0.40 # Proporción poblacional
n <- 20 # Tamaño de la muestra
# Verificación de supuestos
np <- n * p
nq <- n * (1 - p)
# Mostrar resultados
np
## [1] 8
nq
## [1] 12
Se puede decir que si cumplen el supuesto de normalidad.
# Definir los valores
p <- 0.40 # Proporción poblacional
n <- 20 # Tamaño de la muestra
# Calcular la desviación estándar de la proporción muestral
sigma_p <- sqrt((p * (1 - p)) / n)
# Definir el valor de la proporción muestral a evaluar
p_muestra <- 0.50
# Calcular el puntaje Z
z <- (p_muestra - p) / sigma_p
# Calcular la probabilidad acumulada P(p_muestra < 0.50)
probabilidad <- pnorm(z)
# Mostrar resultado en porcentaje
paste("Probabilidad:", round(probabilidad * 100, 2), "%")
## [1] "Probabilidad: 81.93 %"
Este resultado confirma que el 40% de la población mayor de 40 años es consistente con la muestra, ya que no es raro obtener una proporción muestral menor al 50%. En consecuencia, no hay razones para dudar del valor poblacional asumido.
Una empresa ha recibido 120 solicitudes de trabajo de recién graduados en administración de empresas. Se sabe que el 40% de los administradores recién graduados son mujeres. Queremos calcular la probabilidad de que entre el 35% y el 45% de las solicitudes sean de mujeres.
# Parámetros del problema
n <- 120 # Tamaño de la muestra
p <- 0.40 # Proporción poblacional de mujeres
p1 <- 0.35 # Límite inferior del intervalo
p2 <- 0.45 # Límite superior del intervalo
La proporción muestral \(\hat{p}\) sigue aproximadamente una distribución normal con:
Media: \[ \mu_{\hat{p}} = p = 0.40 \]
Desviación estándar: \[ \sigma_{\hat{p}} = \sqrt{\frac{p(1 - p)}{n}} \]
# Cálculo de la desviación estándar de la distribución de la proporción muestral
sigma_p <- sqrt((p * (1 - p)) / n)
sigma_p
## [1] 0.04472136
Aplicamos la transformación:
\[ Z = \frac{\hat{p} - p}{\sigma_{\hat{p}}} \]
# Cálculo de los valores Z
z1 <- (p1 - p) / sigma_p
z2 <- (p2 - p) / sigma_p
# Mostrar valores Z
z1; z2
## [1] -1.118034
## [1] 1.118034
Usamos la función pnorm
para encontrar la probabilidad
acumulada:
# Probabilidad acumulada
p_z1 <- pnorm(z1) # P(Z <= z1)
p_z2 <- pnorm(z2) # P(Z <= z2)
# Probabilidad del intervalo
probabilidad <- p_z2 - p_z1
# Mostrar resultado
probabilidad
## [1] 0.7364475
La probabilidad de que entre el 35% y el 45% de las solicitudes correspondan a mujeres es 73.64%, lo que equivale aproximadamente a 73.64%.
Esto indica que en el 73.64% de las muestras de 120 solicitudes, el porcentaje de mujeres estará dentro del rango especificado. Como esta probabilidad es alta, podemos concluir que es bastante probable que en una selección aleatoria de solicitudes se cumpla esta condición.
# Instalar y cargar los paquetes necesarios
if(!require(rriskDistributions)) install.packages("rriskDistributions")
## Cargando paquete requerido: rriskDistributions
## Warning: package 'rriskDistributions' was built under R version 4.4.3
if(!require(fitdistrplus)) install.packages("fitdistrplus")
## Cargando paquete requerido: fitdistrplus
## Cargando paquete requerido: MASS
## Cargando paquete requerido: survival
library(rriskDistributions)
library(fitdistrplus)
# Establecer semilla para reproducibilidad
set.seed(42)
# Tamaño de muestra pequeño (n < 30)
n <- 20
nsim <- 1000 # Número de simulaciones Monte Carlo
# Parámetros poblacionales conocidos
sigma_pop <- 1
mu_pop <- 0.5
# Simulación de datos log-normales
lognorm_data <- rlnorm(n, meanlog = mu_pop, sdlog = sigma_pop)
# Ajuste de la distribución log-normal
fit_ln <- fitdist(lognorm_data, "lnorm")
summary(fit_ln)
## Fitting of the distribution ' lnorm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## meanlog 0.691920 0.2860807
## sdlog 1.279392 0.2022891
## Loglikelihood: -47.14487 AIC: 98.28974 BIC: 100.2812
## Correlation matrix:
## meanlog sdlog
## meanlog 1 0
## sdlog 0 1
plot(fit_ln)
# Ajuste con diferentes distribuciones para comparar
ln_dist_names <- c("lnorm", "gamma", "weibull", "exp")
ln_results <- list()
for (dist in ln_dist_names) {
ln_results[[dist]] <- try(fitdist(lognorm_data, dist), silent = TRUE)
}
# Comparar ajustes usando AIC (menor es mejor)
ln_aic <- sapply(ln_results, function(x) {
if (inherits(x, "try-error")) return(NA)
return(x$aic)
})
print("Comparación de AIC para datos log-normales:")
## [1] "Comparación de AIC para datos log-normales:"
print(ln_aic)
## lnorm gamma weibull exp
## 98.28974 97.26801 97.17829 95.39267
cat("Mejor distribución:", names(ln_aic)[which.min(ln_aic)], "\n")
## Mejor distribución: exp
# Simulación Monte Carlo para la distribución de medias
medias_lognorm <- replicate(nsim, mean(sample(lognorm_data, n, replace = TRUE)))
###GUI rriskDistributions. fit.cont
fit.cont(lognorm_data)
##
## Begin fitting distributions ---------------------------------------
## * fitting normal distribution ... OK
## * fitting Cauchy distribution ... OK
## * fitting logistic distribution ... OK
## * fitting beta distribution ... failed
## * fitting exponential distribution ... OK
## * fitting chi-square distribution ... OK
## * fitting uniform distribution ... OK
## * fitting gamma distribution ... OK
## * fitting lognormal distribution ... OK
## * fitting Weibull distribution ... OK
## * fitting F-distribution ... OK
## * fitting Student's t-distribution ... OK
## * fitting Gompertz distribution ... failed
## * fitting triangular distribution ... failed
## End fitting distributions -----------------------------------------
## logL AIC BIC Chisq(value) Chisq(p) AD(value) H(AD)
## Normal -56.88 117.76 119.75 22.42 0.00 1.69 rejected
## Cauchy -53.12 110.25 112.24 11.06 0.01 2.12 not rejected
## Logistic -55.55 115.1 117.09 17.24 0.00 1.33 rejected
## Exponential -46.7 95.39 96.39 7.21 0.12 0.45 not rejected
## Chi-square -49.44 100.88 101.87 12.00 0.02 1.26 NULL
## Uniform NULL NULL NULL 27.28 0.00 Inf NULL
## Gamma -46.63 97.27 99.26 6.85 0.08 0.37 NA
## Lognormal -47.14 98.29 100.28 5.44 0.14 0.44 not rejected
## Weibull -46.59 97.18 99.17 6.66 0.08 0.35 not rejected
## F -49.41 102.81 104.81 6.48 0.09 1.63 NULL
## Student -64.21 130.42 131.41 31.23 0.00 12.30 NULL
## KS(value) H(KS)
## Normal 0.27 NULL
## Cauchy 0.22 NULL
## Logistic 0.20 NULL
## Exponential 0.17 NULL
## Chi-square 0.20 NULL
## Uniform 0.25 NULL
## Gamma 0.15 NULL
## Lognormal 0.15 NULL
## Weibull 0.14 NULL
## F 0.29 NULL
## Student 0.57 NULL
##
## Chosen continuous distribution is: NA
## Fitted parameters are:
## [1] NA
# Parámetros de la distribución beta
shape1 <- 2
shape2 <- 5
beta_data <- rbeta(n, shape1, shape2)
# Ajuste de la distribución beta
fit_b <- fitdist(beta_data, "beta")
summary(fit_b)
## Fitting of the distribution ' beta ' by maximum likelihood
## Parameters :
## estimate Std. Error
## shape1 2.018552 0.5985295
## shape2 3.552095 1.1098809
## Loglikelihood: 6.086597 AIC: -8.173194 BIC: -6.181729
## Correlation matrix:
## shape1 shape2
## shape1 1.0000000 0.8267198
## shape2 0.8267198 1.0000000
plot(fit_b)
# Ajuste con otras distribuciones
beta_dist_names <- c("beta", "norm", "gamma", "unif")
beta_results <- list()
for (dist in beta_dist_names) {
beta_results[[dist]] <- try(fitdist(beta_data, dist), silent = TRUE)
}
# Comparar ajustes usando AIC y BIC
beta_aic <- sapply(beta_results, function(x) {
if (inherits(x, "try-error")) return(NA)
return(x$aic)
})
beta_bic <- sapply(beta_results, function(x) {
if (inherits(x, "try-error")) return(NA)
return(x$bic)
})
print("Comparación de AIC para datos beta:")
## [1] "Comparación de AIC para datos beta:"
print(beta_aic)
## beta norm gamma unif
## -8.173194 -4.125469 -7.416622 -14.870173
cat("Mejor distribución (AIC):", names(beta_aic)[which.min(beta_aic)], "\n")
## Mejor distribución (AIC): unif
print("Comparación de BIC para datos beta:")
## [1] "Comparación de BIC para datos beta:"
print(beta_bic)
## beta norm gamma unif
## -6.181729 -2.134005 -5.425157 -12.878708
cat("Mejor distribución (BIC):", names(beta_bic)[which.min(beta_bic)], "\n")
## Mejor distribución (BIC): unif
# Simulación Monte Carlo para la distribución de medias
medias_beta <- replicate(nsim, mean(sample(beta_data, n, replace = TRUE)))
###GUI rriskDistributions. fit.cont
fit.cont(beta_data)
##
## Begin fitting distributions ---------------------------------------
## * fitting normal distribution ... OK
## * fitting Cauchy distribution ... OK
## * fitting logistic distribution ... OK
## * fitting beta distribution ... OK
## * fitting exponential distribution ... OK
## * fitting chi-square distribution ... OK
## * fitting uniform distribution ... OK
## * fitting gamma distribution ... OK
## * fitting lognormal distribution ... OK
## * fitting Weibull distribution ... OK
## * fitting F-distribution ... failed
## * fitting Student's t-distribution ... OK
## * fitting Gompertz distribution ... OK
## * fitting triangular distribution ... failed
## End fitting distributions -----------------------------------------
## logL AIC BIC Chisq(value) Chisq(p) AD(value) H(AD)
## Normal 4.06 -4.13 -2.13 3.38 0.34 0.62 not rejected
## Cauchy -0.57 5.14 7.13 2.68 0.44 0.93 not rejected
## Logistic 3.26 -2.52 -0.53 3.41 0.33 0.58 not rejected
## Beta 6.09 -8.17 -6.18 2.59 0.46 0.38 NULL
## Exponential 0.42 1.16 2.15 2.54 0.64 1.76 rejected
## Chi-square -10.06 22.13 23.12 13.30 0.01 3.96 NULL
## Uniform NULL NULL NULL 1.81 0.61 Inf NULL
## Gamma 5.71 -7.42 -5.43 3.17 0.37 0.31 not rejected
## Lognormal 5.41 -6.83 -4.84 3.79 0.28 0.31 not rejected
## Weibull 5.66 -7.31 -5.32 2.74 0.43 0.36 not rejected
## Student -20.11 42.21 43.21 36.67 0.00 6.22 NULL
## Gompertz 5.04 -6.07 -4.08 2.16 0.54 0.45 NULL
## KS(value) H(KS)
## Normal 0.16 NULL
## Cauchy 0.18 NULL
## Logistic 0.14 NULL
## Beta 0.12 NULL
## Exponential 0.23 NULL
## Chi-square 0.41 NULL
## Uniform 0.13 NULL
## Gamma 0.13 NULL
## Lognormal 0.13 NULL
## Weibull 0.12 NULL
## Student 0.54 NULL
## Gompertz 0.13 NULL
##
## Chosen continuous distribution is: NA
## Fitted parameters are:
## [1] NA
# Configuración de ventana gráfica
par(mfrow = c(3, 2))
# Histograma de datos log-normales con curva ajustada
hist(lognorm_data, main = "Caso 5: Log-normal (datos)", col = "lightblue", prob = TRUE,
xlab = "Valor", ylab = "Densidad")
lines(density(lognorm_data), col = "red", lwd = 2)
x_ln <- seq(min(lognorm_data), max(lognorm_data), length.out = 100)
lines(x_ln, dlnorm(x_ln, meanlog = fit_ln$estimate[1], sdlog = fit_ln$estimate[2]),
col = "blue", lwd = 2, lty = 2)
legend("topright", legend = c("Densidad empírica", "Log-normal ajustada"),
col = c("red", "blue"), lwd = 2, lty = c(1, 2), cex = 0.7)
# Distribución de medias (Monte Carlo) para log-normal
hist(medias_lognorm, main = "Caso 5: Monte Carlo medias", col = "skyblue", prob = TRUE,
xlab = "Media", ylab = "Densidad")
lines(density(medias_lognorm), col = "darkblue", lwd = 2)
media_tcl <- mean(medias_lognorm)
sd_tcl <- sd(medias_lognorm)
x_tcl <- seq(min(medias_lognorm), max(medias_lognorm), length.out = 100)
lines(x_tcl, dnorm(x_tcl, mean = media_tcl, sd = sd_tcl),
col = "purple", lwd = 2, lty = 2)
legend("topright", legend = c("Densidad empírica", "Normal teórica (TCL)"),
col = c("darkblue", "purple"), lwd = 2, lty = c(1, 2), cex = 0.7)
# Histograma de datos beta con curva ajustada
hist(beta_data, main = "Caso 7: Beta (datos)", col = "lightgreen", prob = TRUE,
xlab = "Valor", ylab = "Densidad")
lines(density(beta_data), col = "darkgreen", lwd = 2)
x_beta <- seq(0, 1, length.out = 100)
lines(x_beta, dbeta(x_beta, shape1 = fit_b$estimate[1], shape2 = fit_b$estimate[2]),
col = "green4", lwd = 2, lty = 2)
legend("topright", legend = c("Densidad empírica", "Beta ajustada"),
col = c("darkgreen", "green4"), lwd = 2, lty = c(1, 2), cex = 0.7)
# Distribución de medias (Monte Carlo) para beta
hist(medias_beta, main = "Caso 7: Monte Carlo medias", col = "palegreen", prob = TRUE,
xlab = "Media", ylab = "Densidad")
lines(density(medias_beta), col = "darkgreen", lwd = 2)
media_tcl_beta <- mean(medias_beta)
sd_tcl_beta <- sd(medias_beta)
x_tcl_beta <- seq(min(medias_beta), max(medias_beta), length.out = 100)
lines(x_tcl_beta, dnorm(x_tcl_beta, mean = media_tcl_beta, sd = sd_tcl_beta),
col = "purple", lwd = 2, lty = 2)
legend("topright", legend = c("Densidad empírica", "Normal teórica (TCL)"),
col = c("darkgreen", "purple"), lwd = 2, lty = c(1, 2), cex = 0.7)
# QQ-plot para verificar normalidad en las medias log-normal
qqnorm(medias_lognorm, main = "QQ-plot medias log-normal")
qqline(medias_lognorm, col = "red", lwd = 2)
# QQ-plot para verificar normalidad en las medias beta
qqnorm(medias_beta, main = "QQ-plot medias beta")
qqline(medias_beta, col = "red", lwd = 2)