Distribución binomial en R: dbinom(), pbinom() y qbinom()

Normalmente se empieza a estudiar las distribuciones binomiales con los procesos de Bernoulli. Tales procesos pueden basarse, simple y llanamente, en tirar una moneda al aire y ver si sale cara o cruz. Sus características principales son las siguientes:

1. Constan de ensayos repetidos. Por ejemplo, tirar una moneda al aire 100 veces.
2. Cada ensayo produce un resultado que se puede clasificar dicotómicamente como éxito o fracaso. En el caso de la moneda, que sea cara o cruz.
3. La probabilidad de un éxito, que se denota como p, permanece constante de un ensayo a otro. Si hay un 0.5 de probabilidades de obtener cara en una moneda no trucada, y queremos saber cuantas probabilidades hay de que eso ocurra entonces p = 0.5
4. Los ensayos repetidos son independientes. No cambiará la probabilidad del primer ensayo al segundo ensayo por culpa de lo ocurrido en el primer ensayo.

Podemos lanzar 3 monedas al aire y calcular la probabilidad de obtener 2 caras y elaborar la probabilidad de todos los resultados que puedan darse. Este proceso constaría de 3 ensayos repetidos (tres lanzamientos de moneda) cada ensayo puede clasificarse dicotómicamente como éxito o fracaso (cara sería éxito y cruz fracaso), la probabilidad de éxito y fracaso en este caso sería la misma (p=0.5) — a no ser que la moneda esté trucada — y los ensayos son independientes (la probabilidad de que salga cara en el segundo lanzamiento no dependerá de lo que haya salido en el primero).

Aquí va una tabla en R donde se ilustran los posibles resultados y sus respectivas probabilidades, donde X será el símbolo que hace referencia a cruz y C a cara:

Ncaras <- c(0,1,1,1,2,2,2,3)
Combinaciones <- c('XXX', 'CXX', 'XCX', 'XXC', 'CCX', 'XCC', 'CXC', 'CCC')
Probabilidad <- c(0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125)
Bernoullimonedas <- data.frame(Combinaciones, Ncaras, Probabilidad)

Básicamente se han creado dos vectores y un data.frame. El resultado es el que vemos. Las probabilidades resultan de multiplicar la probabilidad de cada éxito o fracaso en una combinación concreta. Por ejemplo, la probabilidad de que no salga ninguna cara es XXX = 0.5*0.5*0.5=0.125. Y así son todas, ya que éxito y fracaso tienen la misma probabilidad (p=q).

La ecuación para el cálculo de probabilidades de una variable binomial es la siguiente:

p: Probabilidad de éxito

q: Probabilidad de fracaso

n: número de ensayos

x: Valor de la variable aleatoria (en el ejemplo que hemos puesto antes, 0, 1, 2 ó 3 caras)

El elemento

es una partición que vendría a representar el número de combinaciones para un determinado valor de x. Si x es 2, es decir, si queremos saber cuántas combinaciones hay en las que aparezcan 2 caras, esto se puede calcular operando con factoriales:

Pero ese cálculo de probabilidad se puede hacer mucho más rápido en R, lo cual nos permite construir mucho más rápidamente una distribución de probabilidad, es decir, conocer las probabilidades de todos los valores posibles de una variable aleatoria X. Aquí es cuando debemos usar la función dbinom(). Esta función de R nos permite conocer la probabilidad de que nuestra variable aleatoria tome un valor concreto dado un número de ensayos y la probabilidad de éxito. Lo que nos devuelve es un vector con los resultados:

> dbinom(1,3,0.5)
[1] 0.375

Si queremos conocer la probabilidad de más valores, basta con añadir en el parámetro de los valores de X una secuencia con la función seq(), un nuevo vector c(), o simplemente el número mínimo de una secuencia separado por dos puntos del número máximo:

> dbinom(seq(0,3,1),3,0.5)
[1] 0.125 0.375 0.375 0.125
> dbinom(c(0,1,2,3),3,0.5)
[1] 0.125 0.375 0.375 0.125
> dbinom(0:3,3,0.5)
[1] 0.125 0.375 0.375 0.125

Y con la creación de estos vectores podemos graficar la distribución de probabilidad. Para graficarlo con ggplot2, por ejemplo, podemos crear el data.frame “monedas” conformado por dos vectores, uno con una secuencia de 0 a 3, y otro con los valores de la probabilidad de que la variable aleatoria tome valores de 0, 1, 2 y 3 por separado:

vectorbinom <-dbinom(0:100,100,0.5)
posiblesvalores <-c(seq(0,100,1))
monedas<-data.frame(vectorbinom,posiblesvalores)
library(ggplot2)
ggplot(data=monedas, aes(x=posiblesvalores, y=vectorbinom))+
  geom_bar(stat="identity", position="stack")+
  labs(x="Valores de la variable aleatoria", y="Probabilidad")

Con esto podemos visualizar también la aproximación de la distribución de la binomial a una distribución normal conforme el número de ensayos aumenta. Este es el gráfico para n = 100 con una probabilidad de éxito de 0.5:

Siguiendo con el cálculo de probabilidades, R nos permite también calcular la probabilidad de que una variable binomial tome un valor menor o igual a un x determinado. Esto se hace con la función pbinom(). Los argumentos son muy parecidos a la función dbinom(), lo que cambia es el primer argumento, donde ahora debemos introducir el valor que R captará para calcular la probabilidad de que la variable aleatoria tome un valor menor o igual a dicho valor. Por ejemplo, podemos calcular la probabilidad de obtener 2 caras o menos:

> pbinom(2,3,0.5)
[1] 0.875

Cualquiera que sepa un poco de probabilidad, le sonará esto a la distribución de probabilidad acumulada. Y sí, en efecto, el cálculo que acabamos hacer se obtendría de sumar las probabilidades desde x=0 hasta x=2. Es decir, sería el tramo que va desde x=0 a x=2 en una distribución de probabilidad acumulada. De hecho, el mismo resultado se obtendría utilizando la función sum():

> sum(dbinom(0:2,3,0.5))
[1] 0.875

Para graficar una distribución acumulada de probabilidad de un variable binomial, podemos utilizar la función pbinom() poniendo el primer parámetro como variable. Veamos el ejemplo con n=80 y p=0.2:

plot(pbinom(x,50,0.5), type='s', xlab='Probabilidades', ylab='valores de x')

Y, por supuesto, se pueden calcular tramos. Si queremos calcular P(30≤X≤40) podemos hacerlo restando P(X≤40)-P(X≥30):

> pbinom(40,50,0.5) - pbinom(30,50,0.5)
[1] 0.059

Y si queremos calcular P(X≥x) simplemente tendremos que restar el resultado de P(≤X) a 1. Por ejemplo, podemos calcular la probabilidad de P(X≥20) con n=50 y p=0.5:

> 1 - pbinom(20,50,0.5)
[1] 0.898

R también nos permite hacer cálculos a la inversa, es decir, en vez de calcular probabilidades, que nos calcule hasta qué punto concreto de los valores que podría tomar X se acumula una probabilidad concreta. Esto se hace con la función qbinom(). Imaginemos que tenemos un proceso binomial que consta de 50 ensayos y una probabilidad de éxito de 0.5. Podemos pedirle a R que nos diga hasta qué punto concreto de los valores que podría tomar X se acumula el 80% de la probabilidad:

> qbinom(0.8, 50, 0.5)
[1] 28

Esto significa que desde el menor valor que pueda tomar la variable X, hasta el valor 28 tenemos el 0.8 de probabilidad acumulada. Es decir, que hay 0.8 de probabilidades de que la variable obtenga un valor menor o igual a 28.