DISTRIBUCIONES DE PROBABILIDAD EN R

El paquete stats de R (que se instala por defecto al instalar R, y se carga en memoria siempre que iniciamos sesión) implementa numerosas funciones para la realización de cálculos asociados a distintas distribuciones de probabilidad. Entre las utilizadas más comunmente podemos citar:

Algunas Distribuciones de probabilidad

Ejecutando:

help(Distributions)

se nos muestra el listado de distribuciones de probabilidad disponibles en el paquete stats. Otras distribuciones están disponibles en otros paquetes.

Para cada distribución, R dispone de cuatro funciones. Se puede acceder a cada una de ellas simplemente precediendo el nombre de la distribución que figura en la tabla anterior por la letra que se indica a continuación:

d: función de densidad o de probabilidad.

p: función de distribución.

q: función para el cálculo de cuantiles.

r: función para simular datos con dicha distribución.

Así, por ejemplo, para la distribución normal, la función de densidad se obtiene como dnorm(), la función de distribución como pnorm(), los cuantiles se calculan mediante qnorm() y se pueden generar valores aleatorios con distribución normal mediante rnorm(). Puede consultarse la ayuda, help(dnorm) para conocer la sintaxis específica de estas funciones.

Distribuciones discretas de probabilidad

Se denomina distribución de variable discreta a aquella cuya función de probabilidad solo toma valores positivos en un conjunto de valores de \(X\) finito o infinito numerable. A dicha función se le llama función de masa de probabilidad.

Distribución Binomial

La distribución de Bernoulli es una distribución de probabilidad que consiste en la observación de un experimento aleatorio con dos posibles resultados. A uno de los resultados del experimento se le denomina “éxito”, con probabilidad de ocurrencia \(\pi\), y al otro se le llama “fracaso”, con probabilidad de ocurrencia \(1-\pi\), donde \(\pi\) es un número real tal que \(0<\pi<1\).

La distribución binomial está conformada por \(n\) repeticiones independientes de un experimento de Bernoulli.

Se dice que una v.a.d. \(X\) que asume valores \(0,1,\ldots,n\) tiene una distribución binomial con parámetros \(n\) y \(\pi\), lo que se escribe:

\(X\sim Bin(n,\pi)\), si la f.m.p. de \(X\) está dada por:

\[f_X(x;n,\pi) = \left\{ \begin{array}{ll} \binom{n}{x} \pi^x (1-\pi)^{n-x}, & \hbox{si $x=0,1,\ldots,n$;} \\ 0, & \hbox{en otro caso.} \end{array} \right.\]

donde \(n\) es un en tero positivo y \(\pi\) es un número real tal que \(0<\pi<1\).

Cuando \(n = 1\) la distribución Binomial coincide con la distribución de Bernoulli de parámetro \(\pi\), lo que se escribe: \(X \sim Ber(\pi)\).

Si \(X\) sigue una distribución binomial \(B(n,p)\), entonces:

  • \(P\left(X=k\right)\) = dbinom(k,n,p)

  • \(P\left(X\le k\right)\) = pbinom(k,n,p)

  • \(q_{a}=\min\left\{ x:\, P\left(X\le x\right)\ge a\right\}\) = qbinom(a,n,p)

  • rbinom(m,n,p) genera \(m\) valores aleatorios con esta distribución

Parámetros

Si \(x\) es una v.a. tal que \(X \sim Bin(n,\pi)\), entonces:

  • \(\textsf{E}[X] = n\pi\).

  • \(\textsf{V}[X] = n\pi(1-\pi)\).

Si \(X\approx B(10,0.6)\), tenemos entonces:

  • \(P\left(X=8\right)\)
dbinom(8,10,0.6)
## [1] 0.1209324
  • \(P\left(X\le 8\right)\)
pbinom(8,10,0.6)
## [1] 0.9536426

o también:

sum(dbinom(0:8,10,0.6))
## [1] 0.9536426
  • \(q_{0.95}=\min\left\{ x:\, P\left(X\le x\right)\ge 0.95\right\}\)
qbinom(0.95,10,0.6)
## [1] 8

Podemos obtener simultáneamente varios cuantiles:

qbinom(c(0.05,0.95),10,0.6)
## [1] 3 8
  • Simulamos 15 valores de esta distribución:
rbinom(15,10,0.6)
##  [1] 6 5 7 4 5 6 5 6 7 6 6 7 6 4 7

Podemos representar fácilmente la función de probabilidad de la distribución binomial:

plot(dbinom(0:10,10,0.6),type="h",xlab="k",ylab="P(X=k)",
     main="Función de Probabilidad B(10,0.6)")

También podemos representar su función de distribución:

plot(stepfun(0:10,pbinom(0:11,10,0.6)),xlab="k",ylab="F(k)",
     main="Función de distribución B(10,0.6)")

Si simulamos una muestra muy grande de esta distribución (por ejemplo, 10.000 valores), podemos comprobar como las frecuencias relativas son muy similares a las probabilidades teóricas:

X=rbinom(10000, 10, 0.6)
freqAbs=table(X)
freqAbs
## X
##    0    1    2    3    4    5    6    7    8    9   10 
##    2   12  115  424 1160 2024 2471 2134 1222  376   60
freqRel=prop.table(freqAbs)
freqRel
## X
##      0      1      2      3      4      5      6      7      8      9     10 
## 0.0002 0.0012 0.0115 0.0424 0.1160 0.2024 0.2471 0.2134 0.1222 0.0376 0.0060

Vamos a mostrar lado a lado las frecuencias relativas y las probabilidades teóricas. Para ello, generamos primero un dataframe con la tabla de probabilidades teóricas:

probsTeo=data.frame(X=0:10,Prob=dbinom(0:10,10,0.6))
probsTeo
##     X         Prob
## 1   0 0.0001048576
## 2   1 0.0015728640
## 3   2 0.0106168320
## 4   3 0.0424673280
## 5   4 0.1114767360
## 6   5 0.2006581248
## 7   6 0.2508226560
## 8   7 0.2149908480
## 9   8 0.1209323520
## 10  9 0.0403107840
## 11 10 0.0060466176

Para presentar en una única tabla las probabilidades teóricas y las frecuencias relativas de nuesta simulación, podemos utilizar la función merge. Esta función combina dataframes que tengan un campo en común (en este caso el valor de la variable \(X\)). El objeto probsTeo es obviamente un dataframe. Ahora bien, mediante class(freqRel) podemos comprobar que freqRel es un objeto de clase table. Para poder combinarlo con probsTeo hemos de convertirlo primero en dataframe:

freqRel=as.data.frame(freqRel)
str(freqRel)
## 'data.frame':    11 obs. of  2 variables:
##  $ X   : Factor w/ 11 levels "0","1","2","3",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Freq: num  0.0002 0.0012 0.0115 0.0424 0.116 ...

Vemos aquí que la variable \(X\) es de tipo factor; para poder combinar bien este dataframe con probsTeo hemos de convertir \(X\) a la clase integer (ya que en ese dataframe, \(X\) es entera, como puede comprobarse mediante str(probsTeo)):

freqRel$X=as.integer(as.character(freqRel$X))

Por último procedemos a combinar ambos objetos mediante la función merge y mostramos el resultado:

compara=merge(freqRel,probsTeo,all=TRUE)
compara
##     X   Freq         Prob
## 1   0 0.0002 0.0001048576
## 2   1 0.0012 0.0015728640
## 3   2 0.0115 0.0106168320
## 4   3 0.0424 0.0424673280
## 5   4 0.1160 0.1114767360
## 6   5 0.2024 0.2006581248
## 7   6 0.2471 0.2508226560
## 8   7 0.2134 0.2149908480
## 9   8 0.1222 0.1209323520
## 10  9 0.0376 0.0403107840
## 11 10 0.0060 0.0060466176

Podemos hacer un gráfico que nos muestre la similitud entre probabilidad y frecuencia relativa:

with(compara,{
  plot(X,Freq, type="b")
  points(X,Prob,col="red",pch=4)
  lines(X,Prob,col="red",lty=2,lwd=2)
  legend("topleft",c("frec. relativa","probabilidad"),col=c("black","red"),lty=1:2,pch=c(1,4))
})

Ejemplo

Una empresa farmacéutica desarrolló un nuevo medicamento y lo suministró a 10 enfermos elegidos aleatoriamente. La experiencia ha demostrado que 30% de las personas que padecen la enfermedad se recupera al tomar dicho medicamento. ¿Cuál es la probabilidad de que por lo menos nueve de las 10 personas que toman el medicamento se recuperen?

En este caso se tiene que la v.a. de estudio es el “número de personas en la muestra de 10 pacientes que se recupera de la enfermedad” y el éxito consiste en recuperarse de la enfermedad y esto ocurre con una probabilidad de \(0.3\). En consecuencia, la f.m.p. de \(X\) es:

\[f_X(x;10,0.3) = \left\{ \begin{array}{ll} \binom{10}{x} (0.3)^x (0.7)^{10-x}, & \hbox{si $x=0,1,\ldots,10$;} \\ 0, & \hbox{en otro caso.} \end{array} \right.\]

Además, se pide calcular la probabilidad de que por lo menos nueve de las 10 personas que toman el medicamento se recuperen, esto es, \(\textsf{Pr}(X \geq 9)\):

\[\begin{align*} \textsf{Pr}(X \geq 9) &= \textsf{Pr}(X=9;X=10)\\ &= \textsf{Pr}(X=9) + \textsf{Pr}(X=10) \\ &= \binom{10}{9} (0.3)^9 (0.7)^{10-9} + \binom{10}{10} (0.3)^{10} (0.7)^{10-10} \\ %&= 0.000138 + 0.000006 \\ &= 0.000144. \end{align*}\]

De otra parte, se observa que:

  • \(\textsf{E}[X] = 10(0.3)=3\)

Este valor indica que se espera la recuperación de 3 enfermos de una muestra aleatoria de 10 pacientes.

  • \(\textsf{V}[X] = 10(0.3)(0.7) = 2.1\)
# parametros
p <- 0.3
n <- 10
x <- 0:n
# P(X >= 9)
sum(dbinom(x = c(9, 10), size = n, prob = p))
## [1] 0.0001436859
# valor esperado
n*p
## [1] 3
# varianza 
n*p*(1-p)
## [1] 2.1

En la siguiente figura se presenta el gráfico de la f.m.p. y de la f.d.a. de una variable con distribución binomial con parámetros \(n = 10\) y \(\pi = 0.3\).

# f.m.p.
fx <- dbinom(x = x, size = n, prob = p)
# f.d.a.
Fx <- pbinom(q = x, size = n, prob = p)
# graficos
par(mfrow = c(1,2))
# f.m.p
plot(x = x, y = fx, xlab = "x", ylab = "f(x)", pch = 15, col = "blue")
segments(x0 = x, y0 = 0, x1 = x, y1 = fx, lwd = 2, col = "blue")
# f.d.a.
plot(x = c(0, x), y = c(0, Fx), type = "s", xlab = "x", ylab = "F(x)", col = "blue", lwd = 2)
points(x, Fx, col = "blue", pch = 15)

Distribución Hipergeométrica

La distribución hipergeométrica surge a partir del “número de éxitos en \(n\) ensayos dependientes de un experimento de Bernoulli”.

Un experimento hipergeométrico con parámetros \(n\), \(M\), y \(N\) está basado en las siguientes condiciones:

  1. Se elige una muestra sin reemplazo de \(n\) elementos de un conjunto compuesto por \(N\) elementos, de los cuales \(M\) tienen una característica de interés.

  2. Cada elemento se puede caracterizar como un “éxito (el elemento tiene la característica de interés) o como un”fracaso” (el elemento no tiene la característica de interés).

Se dice que una v.a. \(X\) tiene una distribución hipergeométrica con parámetros \(n\), \(M\), y \(N\), lo que se escribe: \(X\sim Hg(n,M,N)\), si la f.m.p. de \(X\) está dada por:

\[f_X(x;n,M,N) = \left\{ \begin{array}{ll} \frac{\binom{M}{x}\binom{N-M}{n-x}}{\binom{N}{n}}, & \hbox{si $x=\max\{0, n+M-N\},\ldots,\min\{n,M\}$;} \\ 0, & \hbox{en otro caso.} \end{array} \right.\]

donde \(n\), \(M\) y \(N\) son números enteros no negativos tales que \(n\leq N\) y \(M\leq N\).

Parámetros

Si \(X\) es una v.a. tal que: \(X \sim H(n,M,N)\), entonces:

  • \(\textsf{E}[X] = n \frac{M}{N}\).

  • \(\textsf{V}[X] = n \frac{M}{N} \left(1-\frac{M}{N}\right) \left(\frac{N-n}{N-1}\right)\).

La expresión anterior evidencia que el valor esperado de la distribución binomial y de la distribución hipergeométrica coinciden, mientras que la varianza de las dos distribuciones difieren por el factor \((N-n)/(N-1)\), denominado factor de corrección por población finita.

Ejemplo

Un equipo de trabajo establecido por el Ministerio de Medio Ambiente, programó visitas a dos fábricas para investigar posibles violaciones a los reglamentos para el control de contaminación ambiental. Sin embargo, los recortes presupuestales han reducido drásticamente el tamaño del equipo de trabajo por lo que solamente se podrán investigar cinco de las 25 fábricas. Si se sabe que 10 de las fábricas están operando sin cumplir los reglamentos, calcular la probabilidad de que al menos una de las fábricas muestreadas esté operando en contra del reglamento.

Se define la v.a. \(X\) como el “número de fábricas en la muestra seleccionada que operan sin cumplir los reglamentos”; de acuerdo con las características del problema se supone que el muestreo se hace sin reemplazo y por lo tanto se sigue que: \(X \sim H(5,10,25)\). La probabilidad pedidad es:

\[\begin{align*} \textsf{Pr}(X \geq 1) &= \sum_{i=1}^5 \textsf{Pr}(X=i) \\ &= 1 - \textsf{Pr}(X=0) \\ &= 1 - \frac{\binom{10}{0}\binom{15}{5}}{\binom{25}{5}} \\ &= 0.9434. \end{align*}\]

En consecuencia, la probabilidad de que al menos una de las fábricas muestreadas esté operando en contra al reglamento es 0.9434.

# parametros
n <- 5
M <- 10
N <- 25
# P(X >= 1)
# la parametrizacion de esta rutina es diferente a la presentada en la formula
sum(dhyper(x = 1:5, m = M, n = N-M, k = n))
## [1] 0.9434783

De otra forma:

1 - dhyper(x = 0, m = M, n = N-M, k = n)
## [1] 0.9434783

La siguiente figura presenta el gráfico de la f.m.p. y de la f.d.a. de una variable con distribución hipergeométrica con parámetros \(n = 5\), \(M = 10\) y \(N = 25\).

# parametros
n <- 5
M <- 10
N <- 25
x <- 0:5
# f.m.p.
fx <- dhyper(x = x, m = M, n = N-M, k = n)
# f.d.a.
Fx <- phyper(q = x, m = M, n = N-M, k = n)
# graficos
par(mfrow = c(1,2))
# f.m.p
plot(x = x, y = fx, xlab = "x", ylab = "f(x)", pch = 15, col = "blue")
segments(x0 = x, y0 = 0, x1 = x, y1 = fx, lwd = 2, col = "blue")
# f.d.a.
plot(x = c(0, x), y = c(0, Fx), type = "s", xlab = "x", ylab = "F(x)", col = "blue", lwd = 2)
points(x, Fx, col = "blue", pch = 15)

Distribución Poisson

Diseñada por el matemático francés Possion en 1837. Para denotar que una variable aleatoria \(𝑋\) sigue una distribución de Poisson, se usa 𝑋 ~ 𝑃𝑜𝑖(𝜆𝑖 ), 𝑖 = 1, ⋯ , 𝑁.La distribución de Poisson se usa para describir el número de ocurrencias de eventos en un espacio de observación limitado. Por ejemplo, la distribución de Poisson puede describir la cantidad de defectos en el sistema mecánico de una aeronave o la cantidad de llamadas a un centro de llamadas en una hora. La distribución de Poisson se usa comúnmente para control de calidad, investigación de confiabilidad / vida útil y seguros. La ocurrencia de eventos de un cierto tipo en el tiempo o un espacio sigue un Proceso de Poisson con tasa de ocurrencia λ>0 por unidad de tiempo o espacio,

La distribución Poisson se utiliza para caracterizar probabilísticamente el número de veces que ocurre un evento en relación con una unidad de medida bien definida (como una unidad de tiempo o espacio, por ejemplo), de forma que:

  1. La probabilidad de que el evento ocurra en una unidad de medida dada es igual para todas las unidades.

  2. El número de eventos que ocurren en una unidad de medida es independiente del número de eventos que ocurren en otras unidades.

Se dice que una v.a. \(X\) tiene distribución de Poisson de parámetro \(\lambda\), se escribe \(X\sim Pois(\lambda).\), si la f.m.p. de \(X\) está dada por:

\[f_X(x;\lambda) = \left\{ \begin{array}{ll} \frac{e^{-\lambda} \lambda^x}{x!}, & \hbox{si $x=0,1,2,\ldots$;} \\ 0, & \hbox{en otro caso.} \end{array} \right.\]

donde \(\lambda\) es un número real positivo.

Parámetros

Si \(X\) es una v.a. tal que \(X \sim P(\lambda)\), entonces:

  • \(\textsf{E}[X] = \lambda\)

  • \(\textsf{V}[X] = \lambda\)

Ejemplo

Los pacientes que entran a un centro de salud lo hacen a una tasa esperada de 0.50 clientes por minuto. Hallar la probabilidad de que el número de clientes que entran en un intervalo específico de 10 minutos sea a lo más 3.

Las hipótesis del proceso de Poisson parecen ser razonables en este contexto. Se da por sentado que los pacientes no llegan en grupos (o es posible contar al grupo entero como un solo paciente) y que la entrada de un paciente no aumenta ni disminuye la probabilidad de que llegue otro.

Para obtener \(\lambda\), se observa que a una tasa media de 0.50 por minuto durante un periodo de 10 minutos, se sigue que \(\lambda = (0.50)(10) = 5\) entradas. Sea \(X\) la v.a. dada por el “número de pacientes que entran en un intervalo de 10 minutos”. Así, se tiene que \(X \sim Pois(5)\), por lo que la f.m.p de \(X\) es:

\[f_X(x;5) = \left\{ \begin{array}{ll} \frac{e^{-5} 5^x}{x!}, & \hbox{si $x=0,1,2,\ldots$;} \\ 0, & \hbox{en otro caso.} \end{array} \right.\]

Se pide calcular:

\[\begin{align*} \textsf{Pr}(X \leq 3) &= \textsf{Pr}(X=0;X=1;X=2;X=3) \\ &= \textsf{Pr}(X=0) + \textsf{Pr}(X=1) + \textsf{Pr}(X=2) + \textsf{Pr}(X=3) \\ &= \frac{e^{-5} 5^0}{0!} + \frac{e^{-5} 5^1}{1!} + \frac{e^{-5} 5^2}{2!} + \frac{e^{-5} 5^3}{3!} \\ &= 0.2650. \end{align*}\]

Se observa además que:

\[\textsf{Pr}(X > 3) = 1 - \textsf{Pr}(X \leq 3) = 1 - 0.2650 = 0.7350.\]

# parametros
lambda <- 5
x <- 0:20
# P(X <= 3)
sum(dpois(x = 0:3, lambda = lambda))
## [1] 0.2650259
# otra manera
ppois(q = 3, lambda = lambda)
## [1] 0.2650259

En siguiente figura presenta el gráfico de la f.m.p. y la f.d.a. de una variable con distribución de poisson de parámetro \(\lambda=5\).

# f.m.p.
fx <- dpois(x = x, lambda = lambda)
# f.d.a.
Fx <- ppois(q = x, lambda = lambda)
# graficos
par(mfrow = c(1,2))
# f.m.p
plot(x = x, y = fx, xlab = "x", ylab = "f(x)", pch = 15, col = "red")
segments(x0 = x, y0 = 0, x1 = x, y1 = fx, lwd = 2, col = "red")
# f.d.a.
plot(x = c(0, x), y = c(0, Fx), type = "s", xlab = "x", ylab = "F(x)", col = "red", lwd = 2)
points(x, Fx, col = "red", pch = 15)

Distribución Geométrica

Para denotar que una variable aleatoria 𝑋 sigue una distribución geométrica, se usa 𝑋 ~ 𝐺𝑒𝑜𝑚(𝑝), 𝑖 = 1, ⋯ , 𝑁. La fórmula para hallar la distribución geométrica es:

\[P(X = x)=p(1-p)^{x-1}\]

Donde: 𝑥 = 0,1,2, …, es el número de fallas en una secuencia antes de que ocurra el primer éxito. \(𝑝\) es la probabilidad de éxito en cada prueba.

\(X\) es el número de ensayos hasta obtener el primer éxito. Pruebas independientes.

En R tenemos las siguientes funciones para la distribución Geométricas:

Funciones distribución geométrica

Los argumentos de estas funciones son:

  • x: vector de cuantiles que representa el número de fallos antes del primer éxito. prob: Probabiidad de éxito en cada ensayo.

  • p: vector de probabilidades.

  • n: número de valores aleatorios por devolver. log,log.p: parámetro booleano. Si es TRUE, las probabilidades p son devueltas como log(p). lower.tail: parámetro booleano, por defecto es TRUE. Las probabilidades son P(X < x), de lo contrario P(X > x).

Parámetros

  • \(\mu=\frac{1}{p}\)

  • \(V(X)=\frac{(1-p)}{p^2}\)

Ejemplo

Un experto tirador acierta en el blanco el 95 % de veces. ¿Cuál es la probabilidad de que falle por primera vez en su decimoquinto disparo?

\(P(X=15) = 0.05(1-0.05)^{14} = 0.05(.95)^{14} = 0.0244 = 2.44%\)

En R

x<-14
p<-0.05
pgeomt<-dgeom(x,p)
paste0("La probabilidad de que falle en la ",x+1,"° vez es ", 
       format(100*pgeomt, digits = 3), "%.")
## [1] "La probabilidad de que falle en la 15° vez es 2.44%."

Distribución Binomial negativa

La distribución binomial negativa es una distribución discreta, que simula el número de ensayos necesarios para producir un número específico de eventos. Cada prueba tiene dos posibles resultados. La distribución binomial negativa también puede modelar el número de no eventos que deben ocurrir para observar un número específico de resultados. Para representar que una variable aleatoria \(𝑋\) sigue una distribución binomial negativa, se usa: \(X \thicksim BN(k,p)\)

Donde \(𝑋\) es el número de ensayos independientes hasta obtener el 𝒌-ésimo éxito y \(𝑝\), la probabilidad de éxito en los ensayos.

\[P(X = x) = {x-1 \choose k-1}p^k(1-p)^{x-k}\]

\[0 < p <1\]

Parámetros

  • \(\mu=\frac{k(1-p)}{p}\) si se piensa en el número de fracasos.

  • \(\mu=\frac{k}{p}\) si se cuenta también los \(k-1\) éxitos.

  • \(V(X)=\frac{k(1-p)}{p^2}\)

En R tenemos las siguientes funciones para la distribución binomiales negativas:

Funciones distribución binomial negativa

Donde: - X = es el número de pruebas falladas.

  • size = es el número de ensayos.

  • prob = es la probabilidad, debe estar entre 0 y 1.

  • mu = es la parametrización alternativa por la media.

  • q = vector de cuantiles.

  • p = vector de probabilidades.

  • log, log.p = parámetro booleano. Si es TRUE, las probabilidades p son devueltas como log(p).

  • lower.tail = parámetro booleano. Por defecto es TRUE. Las probabilidades son \(P(X<x)\) de lo contrario \(P(X>x)\).

Ejemplo

Una refrigeradora tiene una probabilidad de 0.95 % de pasar un control de calidad. Se asume que los resultados del control de calidad de diferentes refrigeradoras son independientes. Calcular la probabilidad de que sea necesario revisar 20 refrigeradoras para que 15 pasen el control de calidad.

\(X\): número de refrigeradoras hasta que 15 pasen el control de calidad. Los datos que se tienen son \(x=20\), \(k=15\), \(p=0.95\).

\[P(X=20)={20-1 \choose 15-1}(0.95)^{15}(1-0.95)^{20-15} \\ P(X=20)={19 \choose 14}(0.95)^{15}(1-0.95)^{5} \\ P(X=20)=\frac{19!}{14!5!}(0.95)^{15}(1-0.95)^5 \\ P(X=20)=0.001683=0.1683%\]

En R

x<-20
k<-15
p<-0.95
pbn<-dnbinom(x-k,size = k, prob = p)
paste0("La probabilidad de que será necesaria es ",format(100*pbn,digits = 3),"%.")
## [1] "La probabilidad de que será necesaria es 0.168%."