InicioEstad铆stica con RDistribuci贸n normal en R

Distribuci贸n normal en R

Crea y dibuja distribuciones Normales en R

La distribuci贸n normal o gaussiana es la distribuci贸n m谩s conocida e importante en Estad铆stica. En este tutorial aprender谩s qu茅 son y c贸mo utilizar las funciones dnorm, pnorm, qnorm y rnorm en R y las diferencias entre ellas. En consecuencia, aprender谩s c贸mo crear y dibujar la distribuci贸n normal en R o RStudio, calcular las probabilidades bajo las curvas, los cuantiles, generaci贸n de n煤meros pseudoaleatorios normales e incluso c贸mo sombrear un 谩rea espec铆fica bajo una curva normal.

La distribuci贸n gaussiana

Entre las variables aleatorias continuas, la m谩s importante es la distribuci贸n normal o gaussiana. Esta variable fue introducida por Carl Friedrich en el siglo XIX para estudiar las medidas de error.

Sea X \sim N(\mu, \sigma), es decir, una variable aleatoria con distribuci贸n normal de media \mu y desviaci贸n t铆pica \sigma:

  • La funci贸n de densidad (o PDF, por sus siglas en ingl茅s), tambi茅n conocida como campana de Gauss, de x es f(x) = \frac{1}{\sqrt{2\pi \sigma^{2}}} e^{\frac{1}{2} (\frac{x - \mu}{\sigma})^2}.
  • La funci贸n de distribuci贸n acumulada (o CDF, por sus siglas en ingl茅s) es F(x) = P(X \leq x).
  • La funci贸n cuantil es Q(p) = F^{-1}(p).
  • La esperanza y la varianza es E(X) = \mu y Var(X) = \sigma^2, respectivamente.

En R existen las funciones dnorm, pnorm y qnorm, que permiten calcular la funci贸n de densidad, la de distribuci贸n y la cuantil normal para una serie de valores. Adem谩s, la funci贸n rnorm permite obtener observaciones aleatorias siguiendo una distribuci贸n normal. Las funciones relacionadas con la distribuci贸n normal se describen en la siguiente tabla:

Funci贸nDescripci贸n
dnormDensidad normal
(Funci贸n de densidad de probabilidad)
pnormDistribuci贸n normal
(Funci贸n de distribuci贸n acumulada)
qnormFunci贸n cuantil de la distribuci贸n normal
rnormGeneraci贸n de n煤meros
pseudoaleatorios normales
Por defecto, estas funciones consideran la distribuci贸n normal est谩ndar, que tiene media 0 y desviaci贸n t铆pica 1.

Aunque trataremos en detalle cada funci贸n en su correspondiente secci贸n, en la siguiente ilustraci贸n puedes ver la relaci贸n entre las funciones dnorm, pnorm y qnorm:

Funciones dnorm pnorm y qnorm en R

La funci贸n dnorm

En R, puedes hacer uso de la funci贸n dnorm para calcular la funci贸n de densidad normal de par谩metros \mu y \sigma para cualquier valor de x.

dnorm(x,           # Valores del eje X (rejilla)
      mean = 0,    # N煤mero o vector representando la/s media/s
      sd = 1,      # N煤mero o vector representando la/s desviaci贸n/es t铆pica/s
      log = FALSE) # Si TRUE, las probabilidades se devuelven como logaritmos

Considera, por ejemplo, que quieres crear la funci贸n de densidad de probabilidad normal para x \in (-4, 4), con media 1 y desviaci贸n t铆pica de 3. Para calcularla, puedes escribir:

x <- -4:4
# x <- seq(-4, 4, length = 100) # M谩s observaciones
dnorm(x, mean = 1, sd = 3) 
0.03315905 0.05467002 0.08065691 0.10648267 0.12579441
0.13298076 0.12579441 0.10648267 0.08065691

Tambi茅n puedes especificar vectores en los argumentos mean y sd de la funci贸n. En el siguiente ejemplo, el valor de los elementos impares corresponde a \mu = 1, \sigma = 3 y los pares a \mu = 2, \sigma = 4.

dnorm(x, mean = c(1, 2), sd = c(3, 4))
0.03315905 0.04566227 0.08065691 0.07528436 0.12579441
0.09666703 0.12579441 0.09666703 0.08065691

Dibujar la distribuci贸n normal en R

Crear un gr谩fico de distribuci贸n normal en R es f谩cil. Solo necesitas crear una rejilla (grid) para el eje X para el primer argumento de la funci贸n plot y pasar como entrada del segundo la funci贸n dnorm para la rejilla correspondiente. En el siguiente ejemplo mostramos c贸mo dibujar distribuciones normales para diferentes medias y varianzas.

par(mfrow = c(1, 2))

# Rejilla de valores para el eje X
x <- seq(-4, 8, 0.1)

#-----------------------------------------
# Misma desviaci贸n t铆pica, distinta media
#-----------------------------------------
# Media 0, desviaci贸n t铆pica 1
plot(x, dnorm(x, mean = 0, sd = 1), type = "l",
     ylim = c(0, 0.6), ylab = "", lwd = 2, col = "red")
# Media 3, desviaci贸n t铆pica 1
lines(x, dnorm(x, mean = 3, sd = 1), col = "blue", lty = 1, lwd = 2)

# A帽adimos una leyenda
legend("topright", c(expression(paste(, mu, " ", sigma)), "0 1", "3 1"),
       lty = c(0, 1, 1), col = c("blue", "red"), box.lty = 0, lwd = 2)

#-----------------------------------------
# Misma media, distinta desviaci贸n t铆pica
#-----------------------------------------
# Media 1, desviaci贸n t铆pica 1
plot(x, dnorm(x, mean = 1, sd = 1), type = "l",
     ylim = c(0, 1), ylab = "", lwd = 2, col = "red")
# Media 1, desviaci贸n t铆pica 0.5
lines(x, dnorm(x, mean = 1, sd = 0.5), col = "blue", lty = 1, lwd = 2)

# A帽adimos una leyenda
legend("topright", c(expression(paste(, mu, " ", sigma)), "1 1", "1 0.5"),
       lty = c(0, 1, 1), col = c("blue", "red"), box.lty = 0, lwd = 2)

par(mfrow = c(1, 1))
Dibujar la funci贸n de densidad normal en R

La funci贸n pnorm

La funci贸n pnorm permite calcular la funci贸n de distribuci贸n acumulada (CDF, por sus siglas en ingl茅s) de la distribuci贸n normal en R, que es la probabilidad de que la variable X tome valores menores o iguales que x. La sintaxis de la funci贸n es la siguiente:

pnorm(q,
      mean = 0,
      sd = 1,
      lower.tail = TRUE, # Si TRUE, las probabilidades son P(X <= x), o P(X > x) si FALSE
      log.p = FALSE)     # Si TRUE, las probabilidades vienen dadas como logaritmos

Como ejemplo, teniendo en cuenta que la distribuci贸n normal es sim茅trica, la probabilidad de que la variable tome un valor inferior a la media es 0.5:

pnorm(0, mean = 0, sd = 1) # 0.5

Ejemplo de la funci贸n pnorm

Ahora, sup贸n que tienes una m谩quina que empaqueta arroz dentro de cajas. El proceso sigue una distribuci贸n normal y se sabe que la media del peso de cada caja es de 1000 gramos y la desviaci贸n t铆pica es 10 gramos. Puedes dibujar la funci贸n de densidad normal en R escribiendo:

# Media y desviaci贸n t铆pica
mu <- 1000
sigma <- 10

# Grid para una distribuci贸n normal no est谩ndar
x <- seq(-3, 3, length = 100) * sigma + mu

# Funci贸n de densidad
f <- dnorm(x, mu, sigma)

plot(x, f, type = "l", lwd = 2, col = "blue", ylab = "", xlab = "Weight")
abline(v = mu) # L铆nea vertical en la media
Densidad normal y su media

Primero, si quieres calcular la probabilidad de que una caja pese menos de 1010 gramos (P(X < 1010) = P(X \leq 1010)), puedes escribir lo que sigue:

pnorm(1010, mu, sigma) # 0.8413447 o 84.13%
1 - pnorm(1010, mu, sigma, lower.tail = FALSE) # Equivalente

Por lo que la probabilidad de que la caja pese menos de 1010 gramos es 0.8413 o 84.13%, que corresponde a el 谩rea de la siguiente ilustraci贸n:

lb <- min(x) # L铆mite inferior
ub <- 1010   # L铆mite superior

x2 <- seq(min(x), ub, length = 100) # Nueva rejilla
y <- dnorm(x2, mu, sigma) # Densidad

plot(x, f, type = "l", lwd = 2, col = "blue", ylab = "", xlab = "Peso")
abline(v = ub) 

polygon(c(lb, x2, ub), c(0, y, 0), col = rgb(0, 0, 1, alpha = 0.5))
text(995, 0.01, "84.13%")
脕rea bajo la curva normal en R

Como sombrear el 谩rea debajo de la curva normal puede ser complicado y requiere varias l铆neas de c贸digo, hemos creado una funci贸n muy sencilla para lograrlo en una sola l铆nea:

# mean: media de la variable normal
# sd: desviaci贸n t铆pica de la variable normal
# lb: l铆mite inferior del 谩rea
# ub: l铆mite superior del 谩rea
# acolor: color del 谩rea
# ...: argumentos adicionales para ser pasados a la funci贸n lines

normal_area <- function(mean = 0, sd = 1, lb, ub, acolor = "lightgray", ...) {
    x <- seq(mean - 3 * sd, mean + 3 * sd, length = 100) 
    
    if (missing(lb)) {
       lb <- min(x)
    }
    if (missing(ub)) {
        ub <- max(x)
    }

    x2 <- seq(lb, ub, length = 100)    
    plot(x, dnorm(x, mean, sd), type = "n", ylab = "")
   
    y <- dnorm(x2, mean, sd)
    polygon(c(lb, x2, ub), c(0, y, 0), col = acolor)
    lines(x, dnorm(x, mean, sd), type = "l", ...)
}

Como ejemplo, si quieres sombrear el 谩rea entre -1 y 2 de una distribuci贸n normal est谩ndar puedes escribir:

normal_area(mean = 0, sd = 1, lb = -1, ub = 2, lwd = 2, acolor = "white")
Funci贸n para sombrear el 谩rea baja la curva normal en R

Segundo, en caso de que quieras calcular la probabilidad de que una caja pese m谩s de 980 gramos (P(X > 980) = P(X \geq 980)) puedes usar el argumento lower.tail.

pnorm(980, mu, sigma, lower.tail = FALSE) # 0.9772499 o 97.72%
1 - pnorm(980, mu, sigma) # Equivalente
pnorm(1020, mu, sigma)    # Equivalente por simetr铆a

La probabilidad calculada corresponde al 谩rea:

normal_area(mean = mu, sd = sigma, lb = 980, acolor = rgb(0, 0, 1, alpha = 0.5))
text(1000, 0.01, "97.72%")
Ejemplo con la funci贸n pnorm en R

Por 煤ltimo, si quieres calcular la probabilidad de que una caja pese m谩s de 990 gramos y menos de 1000 tienes que calcular P(X \leq 1000) - P(x \leq 990) = P(X < 1000) - P(x <990) y por tanto puedes ejecutar:

pnorm(1000, mu, sigma) - pnorm(990, mu, sigma) # 0.3413447 o 34.13%

Puedes dibujar el 谩rea con el siguiente c贸digo:

normal_area(mean = mu, sd = sigma, lb = 990, ub = 1000,
            acolor = rgb(0, 0, 1, alpha = 0.5))
text(995, 0.01, "34.13%", srt = 90)
Calculando la probabilidad bajo una curva normal
Para variables continuas, como P(X = x) = 0, P(X \geq x) = P(X > x) y P(X \leq x) = P(X < x).

Graficar la funci贸n de distribuci贸n normal acumulada en R

Con la funci贸n pnorm tambi茅n puedes dibujar la funci贸n de densidad acumulada de la distribuci贸n normal en R:

par(mfrow = c(1, 2))

# Rejilla de valores para el eje X
x <- seq(-4, 8, 0.1)

#-----------------------------------------
# Misma desviaci贸n t铆pica, distinta media
#-----------------------------------------
# Media 0, desviaci贸n t铆pica 1
plot(x, pnorm(x, mean = 0, sd = 1), type = "l",
     ylim = c(0, 1), ylab = "", lwd = 2, col = "red")
# Media 3, desviaci贸n t铆pica 1
lines(x, pnorm(x, mean = 3, sd = 1), col = "blue", lty = 1, lwd = 2)

# Leyenda
legend("topleft", c(expression(paste(, mu, " ", sigma)), "0 1", "3 1"),
       lty = c(0, 1, 1), col = c("blue", "red"), box.lty = 0, lwd = 2)

#-----------------------------------------
# Misma media, distinta desviaci贸n t铆pica
#-----------------------------------------
# Media 1, desviaci贸n t铆pica 1
plot(x, pnorm(x, mean = 1, sd = 1), type = "l",
     ylim = c(0, 1), ylab = "", lwd = 2, col = "red")
# Media 0, desviaci贸n t铆pica 0.5
lines(x, pnorm(x, mean = 1, sd = 0.5), col = "blue", lty = 1, lwd = 2)

# Leyenda
legend("topleft", c(expression(paste(, mu, " ", sigma)), "1 1", "1 0.5"),
       lty = c(0, 1, 1), col = c("blue", "red"), box.lty = 0, lwd = 2)

par(mfrow = c(1, 1))
Gr谩fico de la funci贸n de distribuci贸n normal acumulada en R

Recuerda que P(X < 0) = 0.5 para una distribuci贸n normal est谩ndar:

x <- seq(-4, 4, 0.1)
plot(x, pnorm(x, mean = 0, sd = 1), type = "l",
     ylim = c(0, 1), ylab = "P(X < x)", lwd = 2, col = "red")
segments(0, 0, 0, 0.5, lwd = 2, lty = 2)
segments(-4, 0.5, 0, 0.5, lwd = 2, lty = 2)
Distribuci贸n normal acumulada o CDF en R

La funci贸n qnorm

La funci贸n qnorm permite encontrar el cuantil (percentil) Q para cualquier probabilidad p.

En consecuencia, la funci贸n qnorm es la inversa de la funci贸n pnorm. La sintaxis de la funci贸n qnorm es como sigue:

qnorm(p,                 # N煤mero o vector de probabilidades
      mean = 0,          # N煤mero o vector representando la/s media/s
      sd = 1,            # N煤mero o vector representando la/s desviaci贸n/es t铆pica/s
      lower.tail = TRUE, # Si TRUE, las probabilidades son P(X <= x), o P(X > x) si FALSE
      log.p = FALSE)     # Si TRUE, las probabilidades vienen dadas como logaritmos

Como primer ejemplo, el cuantil para la probabilidad 0.5 (Q(0.5)) en una distribuci贸n sim茅trica es igual a la media:

qnorm(0.5, mean = 0, sd = 1) # 0

Adem谩s, puedes obtener el cuantil para cualquier probabilidad dada. N贸tese la relaci贸n entre las funciones pnorm y qnorm:

x <- pnorm(-1.5, mean = 0, sd = 1) # 0.0668072
qnorm(x, mean = 0, sd = 1) # -1.5

normal_area(mean = 0, sd = 1, ub = -1.5,
            lwd = 2, acolor = rgb(0, 0, 1, alpha = 0.5))
arrows(-0.5, 0.1, -1.45, 0, lwd = 2, length = 0.2)
text(-0.25, 0.13, "-1.5", cex = 1.5)
Calculando los cuantiles de una distribuci贸n normal en R con la funci贸n qnorm

Si quieres calcular, por ejemplo, el cuantil Q(P(X > 1.5)) = Q(1 - P(X \leq 1.5)) = Q(0.067), puedes establecer el argumento lower.tail de la funci贸n como TRUE.

x <- pnorm(1.5, mean = 0, sd = 1, lower.tail = FALSE) # 0.0668072
qnorm(x, mean = 0, sd = 1, lower.tail = FALSE) # 1.5

# Equivalente a:
# x <- 1- pnorm(1.5, mean = 0, sd = 1) 
# qnorm(1 - x, mean = 0, sd = 1)

normal_area(mean = 0, sd = 1, lb = 1.5, lwd = 2,
            acolor = rgb(0, 0, 1, alpha = 0.5))
arrows(0, 0.1, 1.45, 0, lwd = 2, length = 0.2)
text(0, 0.13, "1.5", cex = 1.5)
qnorm cola derecha de la distribuci贸n normal

Dibujar la funci贸n cuantil normal

Puedes dibujar la funci贸n cuantil de una distribuci贸n normal est谩ndar escribiendo lo siguiente:

plot(qnorm, pnorm(-4), pnorm(4), lwd = 2, xlab = "p", ylab = "Q(p)")

# Equivalente a:
# x <- seq(pnorm(-4), pnorm(4), length = 100)
# plot(x, qnorm(x, mean = 0, sd = 1), type = "l",
#      lwd = 2, xlab = "p", ylab = "Q(p)")
Gr谩fic贸 de la funci贸n cuantil normal en R

La gr谩fica anterior representa los posibles resultados de la funci贸n qnorm para la distribuci贸n normal est谩ndar. Recuerda que Q(0.5) = 0:

plot(qnorm, pnorm(-4), pnorm(4), lwd = 2,
     xlab = "p", ylab = "Q(p)")
segments(0.5, -4, 0.5, 0, lty = 2, lwd = 2)
segments(0, 0, 0.5, 0, lty = 2, lwd = 2)
Gr谩fico de la inversa de la funci贸n de distribuci贸n en R

La funci贸n rnorm

La funci贸n rnorm genera n observaciones de la distribuci贸n normal con media \mu y desviaci贸n t铆pica \sigma . La sintaxis de la funci贸n rnorm en R es la siguiente:

rnorm(n,        # N煤mero de observaciones a ser generadas
      mean = 0, # N煤mero o vector representando la/s media/s
      sd = 1)   # N煤mero o vector representando la/s desviaci贸n/es t铆pica/s

Por lo tanto, puedes generar 10 observaciones de una distribuci贸n normal est谩ndar en R con el siguiente c贸digo:

rnorm(10)
0.3320065 -0.6454991 -0.3692699 -0.4160260  0.5756340
-0.8793241 -0.8095590  0.4926698 -0.8038873 -0.3446890

Sin embargo, debe tenerse en cuenta que si no se especifica una semilla, la salida no ser谩 reproducible. Puedes utilizar la funci贸n set.seed para hacer que la salida sea reproducible:

# Semilla para reproducibilidad
set.seed(1)

rnorm(1) # -0.6264538

Adem谩s, en la siguiente gr谩fica puedes observar c贸mo al aumentar el n煤mero de observaciones, los histogramas de los datos se acercan a la funci贸n de densidad normal real:

# Dividimos la ventana gr谩fica en una fila y tres columnas
par(mfrow = c(1, 3))

x <- seq(-10, 10, length = 200)

# Semilla
set.seed(3)

# n = 10
hist(rnorm(10, mean = 0, sd = 1), main = "n = 10",
     xlab = "", prob = TRUE)
lines(x, dnorm(x), col = "red", lwd = 2)

# n = 100
hist(rnorm(100, mean = 0, sd = 1), main = "n = 100",
     xlab = "", prob = TRUE)
lines(x, dnorm(x), col = "red", lwd = 2)

# n = 1000
hist(rnorm(1000, mean = 0, sd = 1), main = "n = 1000",
     xlab = "", prob = TRUE)
lines(x, dnorm(x), col = "red", lwd = 2)

# Volvemos a la ventana original
par(mfrow = c(1, 1))
Aproximaci贸n de la funci贸n de densidad con observaciones aleatorias de la normal con rnorm