• En este tema estudiaremos diversos tipos de experimentos que son muy frecuentes y algunas de las variables aleatorias asociadas a ellos.

• Estas variables reciben distintos nombres que aplicaremos sin distinción al tipo de población del experimento a la variable o a su función de probabilidad, densidad o distribución.

• Empezaremos con las variables aleatorias discretas que se presentan con frecuencia ya que están relacionadas con situaciones muy comunes como el número de caras en varios lanzamiento de una moneda, el número de veces que una maquina funciona hasta que se estropea, el numero de clientes en una cola,…

Distribución Bernoulli

• Consideremos un experimento con dos resultados posibles éxito (E) y fracaso (F). El espacio de sucesos será \(\Omega=\{E,F\}\).
• Supongamos que la probabilidad de éxito es \(P(E)=p\), y naturalmente \(P(F)=1-p=q\) con \(0<p<1\).
• Consideremos la aplicación

\[ X:\Omega=\{E,F\}\to \mathbb{R} \]

definida por

\[ X(E)=1\mbox{, }X(F)=0. \]

Su función de probabilidad es

\[ P_{X}(x)= \left\{ \begin{array}{ll} 1-p=q & \mbox{si } x=0\\ p & \mbox{si } x=1\\ 0 & \mbox{en cualquier otro caso} \end{array} \right.. \]

Su función de distribución es

\[ F_{X}(x)=P(X\leq x)= \left\{ \begin{array}{ll} 0 & \mbox{si } x<0\\ 1-p=q & \mbox{si } 0\leq x <1\\ 1 & \mbox{si } 1\leq x \\ \end{array} \right.. \]

• Bajo estas condiciones diremos que \(X\) es una v.a. Bernoulli o que sigue una ley de distribución de probabilidad Bernoulli de parámetro \(p\).
• Lo denotaremos por \[X\equiv Ber(p)\mbox{ o también } X\equiv B(1,p).\]
• A este tipo de experimentos (éxito/fracaso)se les denomina experimentos Bernoulli.
• Fue su descubridor un científico suizo Jacob Bernoulli, uno más de la de la conocida familia de científicos suizos Bernoulli

Su valor esperado es

\[E(X)=\displaystyle\sum_{x=0}^1 x\cdot P(X=x)= 0\cdot(1-p)+1\cdot p=p.\]

Calculemos también \(E(X^2)\)

\[E(X^2)=\displaystyle\sum_{x=0}^1 x^2\cdot P(X=x)= 0^2\cdot(1-p)+1^2\cdot p=p.\]

Su varianza es

\[Var(X)=E(X^2)-\left(E(X)\right)^2=p-p^2=p\cdot (1-p)=p\cdot q.\]

Su desviación típica es

\[ \sqrt{Var(X)}=\sqrt{p \cdot (1-p)}. \]

Veamos los cálculos básicos \(Ber(p=0.25)\) en R.

dbinom(0,size=1,prob=0.25)

## [1] 0.75

dbinom(1,size=1,prob=0.25)

## [1] 0.25

rbinom(n=20,size = 1,prob=0.25)

##  [1] 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0

El siguiente código dibuja las función de probabilidad y la de distribución de una \(Ber(p=0.25)\)

par(mfrow=c(1,2))
plot(x=c(0,1),y=dbinom(c(0,1),size=1,prob=0.25),
     ylim=c(0,1),xlim=c(-1,2),xlab="x",
     main="Función de probabilidad\n Ber(p=0.25)")
lines(x=c(0,0,1,1),y=c(0,0.75,0,0.25), type = "h", lty = 2,col="blue")
curve(pbinom(x,size=1,prob=0.25),
      xlim=c(-1,2),col="blue",
      main="Función de distribución\n Ber(p=0.25)")
par(mfrow=c(1,1))

Distribución binomial

Si repetimos \(n\) veces de forma independiente un experimento Bernoulli de parámetro \(p\).

El espacio muestral \(\Omega\) estará formado por cadenas de \(E\)’s y \(F\)’s de longitud \(n\) Consideremos la v.a.

\[X(\overbrace{EFFF\ldots EEF}^{n})=\mbox{número de éxitos en la cadena}.\] A la variable aleatoria anterior se le conoce como distribución binomial de parámetros \(n\) y \(p\), y lo denotaremos por \(X\equiv B(n,p).\)

Entonces su función de probabilidad es

\[ P_{X}(x)=\left\{ \begin{array}{ll} {n\choose x}\cdot p^x \cdot(1-p)^{n-x} &\mbox{ si } x=0,1,\ldots,n\\ 0 & \mbox{ en otro caso} \end{array}\right.. \]

Su función de distribución no tiene una fórmula cerrada. Hay que acumular la función de probabilidad:

\[ \begin{eqnarray*} F_{X}(x)=P(X\leq x) & = & \sum_{i=0}^x P_X(i)\\ & = & \left\{ \begin{array}{ll} 0 & \mbox{ si } x<0\\\displaystyle \sum_{i=0}^k {n\choose i}\cdot p^i \cdot (1-p)^{n-i} & \mbox{ si } \left\{ \begin{array}{l} k\leq x< k+1\\ k=0,1,\ldots,n. \end{array} \right.\\ 1 & \mbox{ si } n\leq x \end{array} \right. \end{eqnarray*} \]

Los números binomiales calculan el número de equipos de baloncesto distintos que (\(k=5\) jugadores) se pueden hacer con 6 jugadores (\(n=6\)).

Es decir cuántas maneras distintas hay para elegir (choose) 5 jugadores en un conjunto de 6 jugadores. Todo el mundo diría ¡¡¡6!!!. Efectivamente con R es

choose(6,5)

## [1] 6

Con 10 jugadores el número de equipos de 5 distintos es bastante más grande

choose(10,5)

## [1] 252

Y, por ejemplo, con un equipo de fútbol profesional que tiene en plantilla 22 jugadores (quitando los guardametas) se pueden formar ¡¡nada menos que!!

choose(22,10)

## [1] 646646

un bonito número capicúa que nos da el número de equipos distintos que se pueden formar.

Obviamente se tiene que una v.a. Bernoulli es una binomial con \(n=1\)

\[B(1,p)=Ber(p).\]

• La probabilidad de fracaso se suele denotar con \(q=1-p\), sin ningún aviso adicional, con el fin de acortar y agilizar la escritura de las fórmulas.
• Su función de distribución no tienen una formula general, hay que calcularla con una función de R o python… En el siglo pasado se tabulaban en los libros de papel :-).
• En el material adicional os pondremos unas tablas de esta distribución para distintos valores de \(n\) y \(p\) para que disfrutéis de tan ancestral método de cálculo.
• Cualquier paquete estadístico, hoja de cálculo dispone de funciones para el cálculo de estas probabilidades, así que el uso de las tablas queda totalmente anticuado.

Su esperanza es \[E(X)=\displaystyle\sum_{k=0}^n k \cdot {n \choose k }\cdot p^k\cdot q^{n-k} = n\cdot p.\]

La esperanza de \(X^2\) es

\[ \begin{eqnarray*} E(X^2)&=& \displaystyle\sum_{k=0}^n k^2 \cdot {n \choose k }\cdot p^k\cdot q^{n-k}\\ &=& n\cdot p\cdot q+(n\cdot p)^2. \end{eqnarray*} \]

Su varianza es

\[Var(X)=E(X^2)-\left(E(X)\right)^2=n\cdot p \cdot q=n\cdot p\cdot (1-p).\]

Su desviación típica es

\[\sqrt{n\cdot p\cdot q}=\sqrt{n\cdot p\cdot (1-p)}.\]

En temas posteriores veremos una forma sencilla del cálculo de la esperanza y varianza de una \(B(n,p)\) como las suma de \(n\) v.a. \(Ber(p)\) independientes.

Veamos los cálculos básicos con funciones de R para una v.a \(X\) con distribución binomial \(B(n=10,p=0.25)\).

Si queremos calcular con R algún valor de la función de distribución como por ejemplo \(F_X(0)=P(X\leq 0)\), tenemos que hacer:

pbinom(0,size=10,prob=0.25)

## [1] 0.05631351

y si queremos por ejemplo \(F_X(4)=P(X\leq 4)\), tenemos que hacer:

pbinom(4,size=10,prob=0.25)

## [1] 0.9218731

Sin embargo, si queremos calcular algún valor de la función de probabilidad como por ejemplo \(P(X=0)\), tenemos que hacer:

dbinom(0,size=10,prob=0.25)

## [1] 0.05631351

o por ejemplo para \(P(X=4)\):

dbinom(4,size=10,prob=0.25)

## [1] 0.145998

Generaremos una muestra aleatoria de 100 valores de una población con distribución \(B(20,0.5)\)

set.seed(2019)
rbinom(100,size = 20,prob=0.5)

##   [1] 12 11  9 11  6  6 12  5  7 11 12 11  8  8 11 11  7 11  9 10  9 10 14  8  8
##  [26]  5 11 14 11 10 11  5 12  8  6  7  9 10  5 12 11  9 12 11 12 10 13 13  8  8
##  [51]  9  7  6  9 10  9 16 13  6  6  8  8 11  9 12 15  9  7 12 11  9  8  9  8 11
##  [76] 15  7 10  9 12  6 13 14  8 10  8 10 11 11  9 10 11 12  8 10 12  9 13  9 13

Veamos los cálculos básicos con funciones de python para una v.a \(X\) con distribución binomial \(B(n=10,p=0.25)\).

Primero importamos la función binom de la librería scipy.stat

from scipy.stats import binom

En general en el paquete scipy, la función de probabilidad se invocará con el método pmf, la de distribución con el método cdf mientras que una muestra aleatoria que siga esta distribución con el método rvs. En todos ellos aparecerá siempre el parámetro loc que se utiliza para desplazar el dominio de la variable aleatoria. Por ejemplo, en este caso

binom.pmf(k, n, p, loc) =  binom.pmf(k - loc, n, p)

Para calcular los valores de la función de distribución como por ejemplo \(F_X(0)=P(X\leq 0)\) y \(F_X(4)=P(X\leq 4)\) utilizamos la función cdf

binom.cdf(0,n=10,p=0.25)

## 0.056313514709472656

binom.cdf(4,n=10,p=0.25)

## 0.9218730926513672

Notemos que al no indicar el valor de loc, se le asume que toma el valor 0.

Para calcular los valores de la función de probabilidad \(P(X=0)\) y \(P(X=4)\) utilizamos la función pmf:

binom.pmf(0,n=10,p=0.25)

## 0.056313514709472684

binom.pmf(4,n=10,p=0.25)

## 0.14599800109863295

Notemos que al no indicar el valor de loc, se le asume que toma el valor 0.

Si queremos generar una muestras aleatorias que siga una distribución binomial, podemos usar la función rvs. En este caso, generaremos una muestra aleatoria de 100 valores de una población \(B(20,0.5)\)

binom.rvs(n=20,p=0.25,size = 100)

## array([ 2,  3,  7,  2,  5,  2,  7,  7,  5,  4,  6,  7,  4,  3,  4,  4,  6,
##         6,  5,  1,  3,  3,  3,  9,  6,  3,  8,  3,  7,  5,  6,  3,  6,  7,
##         5,  1,  4,  8,  6,  7,  3,  3,  7,  4,  4,  3,  8,  4,  6,  2,  4,
##         6,  4,  5,  4,  4,  4,  4,  6,  6,  4,  5,  6,  7,  3,  3,  4,  3,
##         6,  6,  9,  2,  4,  3,  7,  5,  5,  7,  6,  6,  6,  3,  5,  5,  5,
##         5,  5,  4,  2,  4,  5,  5,  5,  6,  6,  7,  5,  6, 10,  3])

Observación Notemos que la secuencia aleatoria generada no es la misma que con R. De hecho, si volvemos a ejecutar esta función obtendremos una muestra aleatoria distinta.

binom.rvs(n=20,p=0.25,size = 100)

## array([ 3,  4,  4,  5,  6,  4,  6,  2,  6, 10,  6,  3,  3,  2,  6,  6,  4,
##         5,  6,  5,  7,  5,  8,  8,  6,  8,  3,  4,  8,  8,  7,  4,  4,  4,
##         3, 10,  5,  3,  2,  2,  3,  2,  8,  3,  2,  5,  7,  2,  7,  2,  3,
##         2,  7,  5,  4,  7, 10,  2,  6,  4,  4,  5, 11,  8,  2,  7,  4,  4,
##         5,  7,  6,  5,  6,  5,  5,  8,  3,  8,  3,  2,  7,  8,  6,  6,  3,
##         6,  4,  5,  4,  3,  5,  5,  6,  6,  7,  8,  6,  4,  5,  3])

Veamos algunos cálculos básicos con funciones de python para la binomial \(B(n=10,p=0.25)\).

binom.cdf(5,n=10,p=0.25)

## 0.9802722930908203

binom.pmf(1,n=10,p=0.25)

## 0.18771171569824247

binom.rvs(n=20,p=0.25,size=10)

## array([10,  7,  5,  5,  5,  6,  3,  6,  3,  5])

El siguiente código de R dibuja las función de probabilidad y la de distribución de una \(B(n=10,p=0.25)\)

par(mfrow=c(1,2))
aux=rep(0,22)
aux[seq(2,22,2)]=dbinom(c(0:10),size=10,prob=0.25)
plot(x=c(0:10),y=dbinom(c(0:10),size=10,prob=0.25),
  ylim=c(0,1),xlim=c(-1,11),xlab="x",
  main="Función de probabilidad\n B(n=10,p=0.25)")
lines(x=rep(0:10,each=2),y=aux, type = "h", lty = 2,col="blue")
curve(pbinom(x,size=10,prob=0.25),
  xlim=c(-1,11),col="blue",
  main="Función de distribución\n B(n=10,p=0.25)")
par(mfrow=c(1,1))

El siguiente código de R dibuja las función de probabilidad y la de distribución de una \(B(n=10,p=0.25)\)

n, p = 10, 0.25
x = np.arange(binom.ppf(0.01, n, p),binom.ppf(0.99, n, p))
fig =plt.figure(figsize=(5, 2.7))
ax = fig.add_subplot(1,2,1)
ax.plot(x, binom.pmf(x, n, p), 'bo', ms=8, label='binom pmf')
ax.vlines(x, 0, binom.pmf(x, n, p), colors='b', lw=5, alpha=0.5)
for tick in ax.xaxis.get_major_ticks():
  tick.label.set_fontsize(5)
for tick in ax.yaxis.get_major_ticks():
  tick.label.set_fontsize(5) 
ax = fig.add_subplot(1,2,2)
ax.plot(x, binom.cdf(x, n, p), 'bo', ms=8, label='binom pmf')
ax.vlines(x, 0, binom.cdf(x, n, p), colors='b', lw=5, alpha=0.5)
for tick in ax.xaxis.get_major_ticks():
  tick.label.set_fontsize(5)
for tick in ax.yaxis.get_major_ticks():
  tick.label.set_fontsize(5)
fig.suptitle('Distribucion Binomial')
plt.show()

• Todos hemos jugado a, por ejemplo, tirar una moneda hasta que obtengamos la primera cara.

• O también tirar una pelota a una canasta de baloncesto hasta obtener la primera canasta.

• Desde otro punto de vista también podemos intentar modelar el número de veces que accionamos una interruptor y la bombilla se ilumina hasta que falla.

• O también el número de veces que un cajero automático nos da dinero hasta que falla.

La modelización de este tipo de problemas se consigue con la llamada distribución geométrica.

Distribución geométrica

• Repitamos un experimento Bernoulli, de parámetro \(p\), de forma independiente hasta obtener el primer éxito.
• Sea \(X\) la v.a. que cuenta el número de fracasos antes del primer éxito. Por ejemplo que hayamos tenido \(x\) fracasos será una cadena de \(x\) fracasos culminada con un éxito. Más concretamente

\[P(\overbrace{FFF\ldots F}^{x}E)=P(F)^{x}\cdot P(E)=(1-p)^{x}\cdot p=q^{x}\cdot p.\]

Su función de probabilidad es

\[ P_X(x)=P(X=x)=\left\{\begin{array}{ll} (1-p)^{x}\cdot p & \mbox{ si } x=0,1,2,\ldots\\ 0 &\mbox{ en otro caso} \end{array}\right.. \]

• La v.a. definida anteriormente diremos que sigue una distribución geométrica de parámetro \(p\).
• La denotaremos por \(Ge(p)\).
• Su dominio es \(D_X=\{0,1,2,\ldots\}\).

Calculemos P(\(X\leq 3\)).

Por la propiedad de la probabilidad del suceso complementario tenemos que

\[ P(X\leq 3 )=1-P(X> 3)=1-P(X\geq 4) \]

Efectivamente, el complementario del evento \(X\leq 3\) nos dice que hemos fracasado más de tres veces hasta conseguir el primer éxito, es decir, hemos fracasado 4 o más veces. Podemos simbolizar dicho evento de la forma siguiente: \[ \{X>3\}=\{X\geq 4\}= \{FFFF\} \]

Ahora, al ser los intentos independientes, tenemos que:

\[ \begin{eqnarray*} P(X>3) & = & P(\{FFFF\})= P(F)\cdot P(F)\cdot P(F)\cdot P(F)\\ &=& (1-p)\cdot (1-p)\cdot (1-p)\cdot (1-p)= (1-p)^{3+1}\\\ &=&(1-p)^{4}. \end{eqnarray*} \]

El valor de la función de distribución de \(X\) en \(x=3\) será, pues: \[F_X(3)=P(X\leq 3)=1-P(X>3)=1-(1-p)^{3+1}.\] Generalizando el resultado anterior a cualquier entero positivo \(k=0,1,2,\ldots\), tenemos: \[F_X(k)=P(X\leq k)=1-(1-p)^{k+1},\mbox{ si } k=0,1,2,\ldots\]

En general, tendremos que: \[ F_X(x)=P(X\leq x)= \left\{\begin{array}{ll} 0, & \mbox{ si } x<0,\\ 1- (1-p), & \mbox{ si } k=0\leq x <1,\\ 1- (1-p)^2, & \mbox{ si } k=1\leq x <2,\\ 1- (1-p)^3, & \mbox{ si } k=2\leq x <3,\\ 1- (1-p)^{k+1}, & \mbox{ si } \left\{ \begin{array}{l}k\leq x< k+1,\\\mbox{para } k=0,1,2,\ldots\end{array} \right.\end{array}\right. \]

De forma más compacta, tendremos que \[ F_X(x)=P(X\leq x)= \left\{\begin{array}{ll} 0, & \mbox{ si } x<0,\\ 1- (1-p)^{k+1}, & \mbox{ si } \left\{ \begin{array}{l}k\leq x< k+1,\\\mbox{para } k=0,1,2,\ldots\end{array} \right.\end{array} \right. \]

Notemos que el límite de la función de distribución es: \[ \displaystyle\lim_{k\to +\infty } F_X(k)=\lim_{k\to +\infty } 1-(1-p)^{k+1}= 1, \] ya que \(0<1-p<1\).

Recordemos del tema de variables aleatorias que

Propiedades

• Si \(|r|<1\) también son convergentes las derivadas, respecto de \(r\), de la serie geométrica y convergen a la derivada correspondiente. Así tenemos que \[ \begin{eqnarray*} \left(\sum_{k=0}^{+\infty} r^k\right)'&= & \sum_{k=1}^{+\infty}k\cdot r^{k-1} &=& \left(\frac1{1-r}\right)'=\frac1{(1-r)^2}\\ \left(\sum_{k=0}^{+\infty} r^k\right)^{''}&=& \sum_{k=2}^{+\infty}k \cdot(k-1)\cdot r^{k-2}&=&\left(\frac1{1-r}\right)^{''}=\frac2{(1-r)^3} \end{eqnarray*} \]

Recordemos que \(P(X=x)=(1-p)^x\cdot p\) si \(x=0,1,2,\ldots\) y aplicado la fórmula anterior con \(r=1-p\)

\[ \begin{eqnarray*} E(X)&=&\sum_{x=0}^{+\infty} x\cdot P_x(x)=\sum_{x=0}^{+\infty} x\cdot (1-p)^x\cdot p\\ &=& p\cdot (1-p) \cdot \sum_{x=1}^{+\infty} x\cdot (1-p)^{x-1}\\ &=& p\cdot (1-p)\cdot \frac{1}{(1-(1-p))^2}=p\cdot (1-p)\cdot \frac{1}{p^2}=\frac{1-p}{p} \end{eqnarray*} \]

\[ \begin{eqnarray*} E(X^2)&=&\sum_{x=0}^{+\infty} x^2\cdot P_X(x)=\sum_{x=1}^{+\infty} x^2\cdot (1-p)^x\cdot p\\ &=& \sum_{x=1}^{+\infty} (x\cdot (x-1)+x)\cdot (1-p)^{x}\cdot p\\ &=& \sum_{x=1}^{+\infty} x\cdot (x-1)\cdot (1-p)^{x}\cdot p+\sum_{x=1}^{+\infty} x \cdot (1-p)^{x}\cdot p\\ &=& (1-p)^{2}\cdot p\cdot \sum_{x=2}^{+\infty} x\cdot (x-1)\cdot (1-p)^{x-2}\\ & +& (1-p)\cdot p\sum_{x=1}^{+\infty} x \cdot (1-p)^{x-1} = \ldots \end{eqnarray*}. \]

\[ \begin{eqnarray*} E(X^2)&=&\ldots\\ &=& (1-p)^{2}\cdot p\cdot \sum_{x=2}^{+\infty} x\cdot (x-1)\cdot (1-p)^{x-2}\\ & +& (1-p)\cdot p\sum_{x=1}^{+\infty} x \cdot (1-p)^{x-1}\\ &=& p\cdot (1-p)^2 \frac{2}{(1-(1-p))^3}+ (1-p)\cdot p \frac{1}{(1-(1-p))^2}\\ &=& p\cdot (1-p)^2 \frac{2}{p^3}+ (1-p)\cdot p \frac{1}{p^2}\\ &=&\frac{2\cdot (1-p)^2}{p^2}+\frac{1-p}{p}. \end{eqnarray*} \]

\[ \begin{eqnarray*} Var(X)&=&E(X^2)-E(X)^2=\frac{2\cdot (1-p)^2}{p^2}+\frac{1-p}{p}-\left(\frac{1-p}{p}\right)^2\\ &=& \frac{2\cdot (1-p)^2+p\cdot(1-p)-(1-p)^2}{p^2}=\frac{(1-p)^2+p\cdot(1-p)}{p^2}\\ &=& \frac{1-2\cdot p + p^2+p-p^2}{p^2}\\ &=& \frac{1-p}{p^2}. \end{eqnarray*} \] Y su desviación típica será

\[\sqrt{Var(X)}=\sqrt{\frac{1-p}{p^2}}.\]

• Supongamos que sólo estamos interesados en el número de intentos para obtener el primer éxito.
• Si definimos \(Y\)= número de intentos para obtener el primer éxito. Entonces \(Y=X+1\) donde \(X\equiv Ge(p)\).
• Su dominio es \(D_Y=\{1,2,\ldots\}\)
• La media se incrementa en un intento debido al éxito \(E(Y)=E(X+1)=E(X)+1=\frac{1-p}{p}+1=\frac1{p}\).
• La varianza es la misma \(Var(Y)=Var(X+1)=Var(X)=\frac{1-p}{p^2}\).

Propiedad de la falta de memoria

Sea \(X\) una v.a. discreta con dominio \(D_X=\{0,1,2,\ldots\}\), con \(P(X=0)=p\).

Entonces \(X\) sigue una ley \(Ge(p)\) si, y sólo si, \[ P\left(X> k+j\big| X\geq j\right)=P(X> k) \] para todo \(k,j=0,1,2,3\ldots\).

Observación: Interpretación de la propiedad

La propiedad de la falta de memoria \[ P(X> k+j\big|X \geq j)=P(X > k), \]
significa que, aunque ya llevemos al menos \(j\) fracasos, la probabilidad de que fracasemos \(k\) veces más no disminuye, es la misma que era cuando empezamos el experimento.

A este efecto se le suele etiquetar con la frase el experimento carece de memoria o es un experimento sin memoria (Memoryless Property).

Un ejemplo muy sencillo nos aclarará el alcance de esta propiedad

Veamos los cálculos básicos con R para la distribución geométrica \(Ge(p=0.25)\). R implementa la geométrica que cuenta el número de fracasos.

\(P(X=0)=(1-0.25)^0\cdot 0.25^1=0.25\)

dgeom(0,prob=0.25)

## [1] 0.25

\(P(X\leq 0)=1- (1-0.25)^{0+1}=1-0.75=0.25\)

pgeom(0,prob=0.25)

## [1] 0.25

\(P(X\leq 4)=1-(1-0.25)^{4+1}=1-0.75=1-0.75^5=0.7626953.\)

pgeom(4,prob=0.25)

## [1] 0.7626953

Una muestra aleatoria de tamaño 25 de una \(Ge(0.25)\)

rgeom(n=25,prob=0.25)

##  [1]  5  4  1  6 10  0  0 10  7  0  6  2  1  3  0  2  5  0  0  5  5  3  3  2  2

par(mfrow=c(1,2))
x=c(0:10)
plot(x=x,y=dgeom(x,prob=0.25),
  ylim=c(0,1),xlim=c(-1,11),xlab="x",
  main="Función de probabilidad\n Ge(p=0.25)")
lines(x=rep(0:10,each=2),y=aux, type = "h", lty = 2,col="blue")
aux0=dgeom(c(0:10),prob=0.25)
ceros=rep(0,21)
ceros
aux=ceros
aux[2*(c(1:11))]<-aux0
curve(pgeom(x,prob=0.25),
  xlim=c(-1,10),col="blue",
  main="Función de distribución\n Ge(p=0.25)")
par(mfrow=c(1,1))

Veamos los cálculos básicos con python para la distribución geométrica \(Ge(p=0.25)\). scipy.stats implementa la distribución geométrica que cuenta el número intentos así que empieza en 1

Cargamos la función de la librería

from scipy.stats import geom

La función de probabilidad es geom.pmf(x,p,loc=0)=geom.pmf(x,p) es un geométrica que cuenta el número de intentos para obtener el primer éxito el valor por defecto del último parámetro es loc=0.

Si queremos la que cuenta el número de fracasos para obtener el primer éxito (la geométrica que empieza en 0) tenemos que usar geom.pmf(x,p,loc=-1).

Es decir geom.pmf(x,p,loc=-1)=geom.pmf(x-1,p,loc=0)

Veamos pues los cálculos para la \(Ge(p)\) que empieza en \(0\).

\(P(X=0)=(1-0.25)^0\cdot 0.25^1=0.25\)

geom.pmf(0,p=0.25,loc=-1)

## 0.25

\(P(X\leq 0)=1- (1-0.25)^{0+1}=1-0.75=0.25\)

geom.cdf(0,p=0.25,loc=-1)

## 0.24999999999999997

\(P(X\leq 4)=1-(1-0.25)^{4+1}=1-0.75=1-0.75^5=0.7626953.\)

geom.cdf(4,p=0.25,loc=-1)

## 0.7626953125

Una muestra aleatoria de tamaño 25 de una \(Ge(0.25)\)

geom.rvs(p=0.25, size=20, loc=-1)

## array([ 0,  3,  0,  5,  2,  3,  0,  1, 10,  0,  8,  2,  4,  4,  0,  4,  5,
##         6,  4,  1])

geom.cdf(range(5),p=0.3,loc=0)

## array([0.    , 0.3   , 0.51  , 0.657 , 0.7599])

geom.cdf(range(5),p=0.3,loc=-1)

## array([0.3    , 0.51   , 0.657  , 0.7599 , 0.83193])

Con python también podemos calcular directamente algunos parámetros asociado a una función de distribución predefinida

geom.stats(p=0.25, loc=0, moments='mv')

## (array(4.), array(12.))

geom.stats(p=0.25, loc=-1, moments='mv')

## (array(3.), array(12.))

p = 0.25
x = np.arange(geom.ppf(0.01, p),geom.ppf(0.99, p))
fig =plt.figure(figsize=(5, 2.7))
ax = fig.add_subplot(1,2,1)
ax.plot(x, geom.pmf(x, p), 'bo', ms=5, label='geom pmf')
ax.vlines(x, 0, geom.pmf(x, p), colors='b', lw=2, alpha=0.5)
for tick in ax.xaxis.get_major_ticks():
  tick.label.set_fontsize(5)
for tick in ax.yaxis.get_major_ticks():
  tick.label.set_fontsize(5) 
ax = fig.add_subplot(1,2,2)
ax.plot(x, geom.cdf(x, p), 'bo', ms=5, label='geom pmf')
ax.vlines(x, 0, geom.cdf(x, p), colors='b', lw=2, alpha=0.5)
for tick in ax.xaxis.get_major_ticks():
  tick.label.set_fontsize(5)
for tick in ax.yaxis.get_major_ticks():
  tick.label.set_fontsize(5)
fig.suptitle('Distribucion Geometrica')
plt.show()

Supongamos que disponemos de 10 llaves distintas y tenemos que abrir una puerta con dos cerraduras.

Comenzamos por la primera cerradura, de tal forma que cada vez olvidamos qué llave hemos probado.

Una vez abierta la primera cerradura probamos de igual forma con la segunda hasta que también la abrimos.

Sea \(X=\) la v.a. que cuenta el número de fracasos hasta abrir la puerta.

Acertar una llave de la puerta es un experimento Bernoulli con probabilidad de éxito \(p=0.1\). Lo repetiremos hasta obtener 2 éxitos.

En general tendremos un experimento de Bernoulli con probabilidad de éxito \(0<p<1\) tal que:

• Repetimos el experimento hasta obtener el \(n\)-ésimo éxito ¡¡abrir la maldita puerta!!.
• Sea \(X\) la v.a. que cuenta el número fallos hasta abrir la puerta, es decir, hasta conseguir el n-ésimo éxito. Notemos que no contamos los éxitos, solo contamos los fracasos

Si representamos como es habitual un suceso como una cadena de F’s y E’s, para \(n=2\), algunos sucesos elementales serán: \[\small{\{EE,FEE,EFE, FFEE,FEFE,EFFE,FFFEE,FFEFE,FEFFE,EFFFE\}.}\]

Calculemos algunas probabilidades para \(n=2\): \[ \small{ \begin{array}{rl} P(X=0) & =P(\{EE\})=p^2, \\ P(X=1) & =P(\{FEE,EFE\})=2\cdot (1-p)\cdot p^2, \\ P(X=2) & =P(\{FFEE,FEFE,EFFE\})=3\cdot (1-p)^2\cdot p^2, \\ P(X=3) & =P(\{FFFEE,FFEFE,FEFFE,EFFFE\})=4\cdot (1-p)^3\cdot p^2. \end{array} } \]

En general su función de probabilidad es

\[ P_{X}(k)=P(X=k)=\left\{\begin{array}{ll} {k+n-1\choose n-1} \cdot (1-p)^{k}\cdot p^n & \mbox{si } k=0,1,\ldots\\ 0 & \mbox{en otro caso}\end{array}\right. \]

Una v.a. con este tipo de distribución recibe el nombre de binomial negativa y la denotaremos por \(BN(n,p)\).

Notemos que \(BN(1,p)=Ge(p)\).

Números binomiales negativos

Dados dos enteros positivos \(n\) y \(k\) se define el número binomial negativo como

\[\binom{-n}{k}=\frac{(-n)(-n-1)\cdots (-n-k+1)}{k!}.\]

Los números binomiales negativos generalizan la fórmula de Newton para exponentes negativos: \[ (t+1)^{-n}=\sum_{k=0}^{+\infty}\left(\begin{array}{c} -n \\ k\end{array}\right) t^{k} \]

R usa la función choose para calcular números binomiales, sean negativos o no. Veámoslo con un ejemplo:

\[ \begin{array}{rl} {-6\choose 4}&=\frac{-6\cdot (-6-1)\cdot \cdot (-6-2)\cdot (-6-3) }{4!}\\ &= \frac{-6\cdot(-7)\cdot (-8)\cdot (-9)}{24}\\ &= \frac{3024}{24}=126. \end{array} \]

Si realizamos el cálculo con R obtenemos el mismo resultado:

choose(-6,4)

## [1] 126

Su esperanza es

\[E(X)=\sum_{k=0}^{+\infty} k\cdot {k+n-1\choose n-1} \cdot (1-p)^{k}\cdot p^n=n\cdot\frac{1-p}{p}.\]

La esperanza de \(X^2\) es

\[E(X^2)=\sum_{k=0}^{+\infty} k^2\cdot {k+n-1\choose n-1} \cdot (1-p)^{k}\cdot p^n=n\cdot\frac{1-p}{p^2}+\left(n\cdot \frac{1-p}{p}\right)^2.\]

Por último la varianza es

\[ Var(X)=E(X^2)-E(X)^2= \]

\[=n\cdot \frac{1-p}{p^2}+\left(n\cdot \frac{1-p}{p}\right)^2-\left(n\cdot \frac{1-p}{p}\right)^2= n\cdot \frac{1-p}{p^2}.\]

y por tanto la desviación típica es

\[\sqrt{Var(X)} = \frac{\sqrt{n(1-p)}}{p}\]

La función de R que calcula la función de probabilidad de la binomial negativa con sus parámetros básicos es:

dnbinom(x, size, prob,...)`

donde size (\(n\)) es el número de éxitos y prob (\(p\)), la probabilidad de éxito.

Así en el ejemplo de la puerta con dos cerraduras, \(X\) es una \(BN(n=size=2,p=prob=0.1)\). Por ejemplo, \(P(X=5)\) que hemos calculado en el ejemplo anterior, vale:

dnbinom(5,size=2,p=0.1)

## [1] 0.0354294

De forma similar calculamos calculamos \(P(X\leq 4)\), \(P(X>4)=1-P(X\leq 4)\) y \(P(X>4)\).

pnbinom(4,size=2,p=0.1)

## [1] 0.114265

1-pnbinom(4,size=2,p=0.1)

## [1] 0.885735

pnbinom(4,size=2,p=0.1,lower.tail=FALSE)

## [1] 0.885735

La función con python es nbinom.pmf(k, n, p, loc). Hay que cargarla desde scpi.stats

from scipy.stats import nbinom

Recordemos que de nuevo se cumple que

nbinom.pmf(k, n, p, loc) = nbinom.pmf(k-loc, n, p)`

nbinom.pmf(k=5,n=2,p=0.1)

## 0.03542940000000002

nbinom.pmf(k=5,n=2,p=0.1,loc=0)

## 0.03542940000000002

nbinom.cdf(k=4,n=2,p=0.1)

## 0.11426500000000003

1-nbinom.cdf(k=4,n=2,p=0.1)

## 0.8857349999999999

Generemos 100 observaciones aleatorias de una \(BN(n=2,0.1)\). Es decir serán las veces que hemos fallado hasta abrir la puerta 100 veces.

nbinom.rvs(n=2, p=0.1, size=100)

## array([35, 14, 19,  1,  5,  1, 41, 26, 36,  2, 12, 15, 16, 10,  3,  2, 36,
##        44, 15,  7, 21,  9, 18, 18, 12,  9, 50, 30,  3,  8,  9, 15, 27, 35,
##         2,  7, 11,  6, 14, 45,  8,  7,  5, 29, 70, 28, 34, 11, 43,  1,  6,
##        15,  1, 16,  9, 23, 21, 13,  3, 18,  5,  3,  4, 34,  5,  2, 14, 12,
##         3, 15, 21, 31, 12, 21,  6,  3, 38, 13, 14,  8, 14, 12, 15,  7, 49,
##        14,  8, 28,  4,  5, 26, 19, 12, 17,  4, 15, 70, 13, 26, 20])

La esperanza y la varianzade una \(BN(n=2,0.1)\) valen:

n, p=2,0.1
params = nbinom.stats(n,p,moments='mv')
print("E(X)={m}".format(m=params[0]))

## E(X)=18.0

print("Var(X)={v}".format(v=params[1]))

## Var(X)=180.0

El siguiente código de R dibuja las función de probabilidad y la de distribución de una \(BN(n=2,p=0.1)\)

par(mfrow=c(1,2))
aux=rep(0,22)
aux[seq(2,22,2)]=dnbinom(c(0:10),size=2,prob=0.1)
plot(x=c(0:10),y=dnbinom(c(0:10),size=2,prob=0.1),
  ylim=c(0,1),xlim=c(-1,11),xlab="x",
  main="Función de probabilidad\n BN(n=2,p=0.1)")
lines(x=rep(0:10,each=2),y=aux, type = "h", lty = 2,col="blue")
curve(pnbinom(x,size=2,prob=0,1),
  xlim=c(-1,11),col="blue",
  main="Función de distribución\n BN(n=2,p=0.1)")
par(mfrow=c(1,1))

El siguiente código de R dibuja las función de probabilidad y la de distribución de una \(BN(n=2,p=0.1)\)

n, p = 10, 0.25
x = np.arange(0,nbinom.ppf(0.99, n, p))
fig =plt.figure(figsize=(5, 2.7))
ax = fig.add_subplot(1,2,1)
ax.plot(x, nbinom.pmf(x, n, p), 'bo', ms=5, label='nbinom pmf')
ax.vlines(x, 0, nbinom.pmf(x, n, p), colors='b', lw=2, alpha=0.5)
for tick in ax.xaxis.get_major_ticks():
  tick.label.set_fontsize(5)
for tick in ax.yaxis.get_major_ticks():
  tick.label.set_fontsize(5) 
ax = fig.add_subplot(1,2,2)
ax.plot(x, nbinom.cdf(x, n, p), 'bo', ms=5, label='nbinom pmf')
ax.vlines(x, 0, nbinom.cdf(x, n, p), colors='b', lw=2, alpha=0.5)
for tick in ax.xaxis.get_major_ticks():
  tick.label.set_fontsize(5)
for tick in ax.yaxis.get_major_ticks():
  tick.label.set_fontsize(5)
fig.suptitle('Distribucion Binomial Negativa')
plt.show()

Diremos que una v.a. discreta \(X\) con \(X(\Omega)=\mathbf{N}\) tiene distribución de Poisson con parámetro \(\lambda>0\), y lo denotaremos por \(Po(\lambda)\) si su función de probabilidad es:

\[ P_{X}(x)=P(X=x)= \left\{\begin{array}{ll} \frac{\lambda^x}{x!} e^{-\lambda}& \mbox{ si } x=0,1,\ldots\\ 0 & \mbox{en otro caso}\end{array}\right.. \]

Usando que el desarrollo en serie de Taylor de la función exponencial es \[ e^{\lambda}=\sum_{x=0}^{+\infty} \frac{\lambda^x}{x!}, \] es fácil comprobar que la suma de la función de probabilidad en todos los valores del dominio de \(X\), o sea, los enteros positivos, vale 1.

Además recordemos que dado \(x\in\mathbb{R}-\{0\}\) se tiene que

\[ \lim_{n\to\infty} \left(1+\frac{x}{n}\right)^n=e^x. \]

Usando la expresión anterior para \(x=-\lambda\), tenemos:

\[ \lim_{n\to\infty} \left(1-\frac{\lambda}{n}\right)^n=\lim_{n\to\infty} \left(1+\frac{-\lambda}{n}\right)^n=e^{-\lambda}. \]

La distribución de Poisson (Siméon Denis Poisson) aparece en el conteo de determinados eventos que se producen en un intervalo de tiempo o en el espacio.

Supongamos que nuestra variable de interés es \(X\), el número de eventos en el intervalo de tiempo \((0,t]\), como por ejemplo el número de llamadas a un call center en una hora donde suponemos que se cumplen las siguientes condiciones:

1. El número promedio de eventos en el intervalo \((0,t]\) es \(\lambda>0\).
2. Es posible dividir el intervalo de tiempo en un gran número de subintervalos (denotemos por \(n\) al número de intervalos) de forma que:
   • La probabilidad de que se produzcan dos o más eventos en un subintervalo es despreciable.
   • El número de ocurrencias de eventos en un intervalo es independiente del número de ocurrencias en otro intervalo.
   • La probabilidad de que un evento ocurra en un subintervalo es \(p_n=\frac{\lambda}{n}\)·

Bajo estas condiciones, podemos considerar que el número de eventos en el intervalo \((0,t]\) será el número de “éxitos” en \(n\) repeticiones independientes de un proceso Bernoulli de parámetro \(p_n\)

Entonces si \(n\to\infty\) y \(p_n\cdot n\) se mantiene igual a \(\lambda\) resulta que la función de probabilidad de \(X\) se puede escribir como

\[ \begin{array}{rl} P(X_n=k)&=\left(\begin{array}{c} n\\ k\end{array}\right) \cdot p_n^k\cdot (1-p_n)^{n-k} \\ &= {n\choose k}\cdot \left(\frac{\lambda}{n}\right)^{k}\cdot \left(1-\frac{\lambda}{n}\right)^{n-k}\\ &= \frac{\lambda^k}{k!}\cdot\frac{n!}{(n-k)!\cdot n^k}\cdot \left(1-\frac{\lambda}{n}\right)^{n}\cdot \left(1-\frac{\lambda}{n}\right)^{-k}. \end{array} \]

Si hacemos tender \(n\) hacia \(\infty\), obtenemos: \[ \lim_{n\to \infty} P(X_n=k) = \lim_{n\to \infty} \frac{\lambda^k}{k!}\cdot\frac{n!}{(n-k)!\cdot n^k} \cdot \left(1-\frac{\lambda}{n}\right)^{n}\cdot \left(1-\frac{\lambda}{n}\right)^{-k}. \]

Calculemos el límite de algunos de los factores de la expresión

\[ \displaystyle\lim_{n\to \infty}\frac{n!}{(n-k)!\cdot n^k}= \lim_{n\to \infty}\frac{n\cdot (n-1)\cdots (n-k-1)}{n^k} =\lim_{n\to \infty}\frac{n^{k}+\cdots}{n^k}=1. \]

\[ \lim_{n\to \infty} \left(1-\frac{\lambda}{n}\right)^{n}=e^{-\lambda} \]

Y también teniendo en cuanta que \(k\) es constante.

\[ \lim_{n\to \infty} \left(1-\frac{\lambda}{n}\right)^{-k}=\lim_{n\to \infty} 1^{-k}=\lim_{n\to \infty} 1=1. \]

Para acabar

\[ \displaystyle\lim_{n\to\infty} P(X_n=k)= \lim_{n\to\infty} \left(\begin{array}{c} n\\ k\end{array}\right) \cdot p_n^k \cdot (1-p_n)^{n-k}= \frac{\lambda^k}{k!}\cdot 1 \cdot e^{-\lambda}\cdot 1=\frac{\lambda^k}{k!}\cdot e^{-\lambda}. \]

Lo que confirma que límite de una serie de variables \(B(n,p_n=\frac{\lambda}{n})\) sigue una ley \(Po(\lambda)\).

Lo interesante de las variables Poisson es que podemos modificar (si el modelo lo permite) el intervalo de tiempo \((0,t]\) en el que contamos los eventos.

Claro que esto no tiene que poder ser así.

Pero en general si la variable es poisson en \((0,t]\) también lo será en cualquier subintervalo \((0,t']\) para todo \(t'\) tal que \(0<t'<t\).

Así que podremos definir una serie de variables \(X_t\) de distribución \(Po(\lambda\cdot t)\).

Definición procesos de Poisson

Consideremos un experimento Poisson con \(\lambda\) igual al promedio de eventos en una unidad de tiempo (u.t.).

Si \(t\) es una cantidad de tiempo en u.t., la v.a. \(X_{t}\)=numero de eventos en el intervalo \((0,t]\) es una \(Po(\lambda\cdot t)\).

El conjunto de variables \(\{X_t\}_{t>0}\) recibe el nombre de proceso de Poisson.

Bajo el punto de vista anterior y si \(p\) es pequeño y \(n\) suficientemente grande la distribución \(B(n,p)\) se aproxima a una \(Po(\lambda=n\cdot p)\).

Existen distintos criterios (ninguno perfecto) de cuando la aproximación es buena.

Por ejemplo si

\[n\geq 20\mbox{ o mejor }n\geq 30, n\cdot p < 10 \mbox{ y } p\leq 0.05,\]

la aproximación de una \(B(n,p)\) por una \(Po(n\cdot p)\) es buena. Sobre todo para los valores cercanos a \(E(X)=\lambda\).

Condición deseable \(n\geq 20\), \(n\cdot p < 10\), \(p\leq 0.05\)

Consideremos por ejemplo una v.a. \(X\) con distribución \(Po(\lambda=3)\). Calculemos \(P_X(0)=P(X=0), P_X(1)=P(X=1)\) con R:

dpois(0,lambda = 3)

## [1] 0.04978707

dpois(1,lambda = 3)

## [1] 0.1493612

Si quisiéramos hallar la función de distribución en los mismos valores anteriores, \(F_X(0)=P(X\leq 0), F_X(1)=P(X\leq 1)\), haríamos lo siguiente:

ppois(0,lambda = 3)

## [1] 0.04978707

ppois(1,lambda = 3)

## [1] 0.1991483

dpois(0,lambda = 3)+dpois(1,lambda = 3) ## es igual a ppois(1,lambda=3)

## [1] 0.1991483

A continuación, comprobemos que \(F_X(10)=\sum\limits_{x=0}^{10} P_X(x)\):

dpois(0:10,3)

##  [1] 0.0497870684 0.1493612051 0.2240418077 0.2240418077 0.1680313557
##  [6] 0.1008188134 0.0504094067 0.0216040315 0.0081015118 0.0027005039
## [11] 0.0008101512

sum(dpois(0:10,3))

## [1] 0.9997077

ppois(10,3)

## [1] 0.9997077

Si quisiéramos generar una secuencia de \(100\) observaciones para una distribución de Poisson de parámetro \(\lambda=3\), \(Po(3)\), tendríamos que hacer:

rpois(n=100,lambda = 3)

##   [1] 2 5 3 3 2 2 5 2 4 4 2 3 2 2 2 2 2 3 3 5 3 3 2 4 2 3 2 1 1 3 4 6 2 5 3 4 1
##  [38] 1 6 3 4 1 4 3 4 3 0 2 1 4 3 0 2 4 2 3 5 2 1 3 3 4 2 5 0 3 1 1 4 6 4 5 0 4
##  [75] 0 3 3 3 4 1 2 6 2 2 2 2 1 2 5 2 5 3 7 3 5 2 3 2 1 3

lambda=20
par(mfrow=c(1,2))
n=qpois(0.99,lambda=lambda)
aux=rep(0,(n+1)*2)
aux[seq(2,(n+1)*2,2)]=dpois(c(0:n),lambda=lambda)
ymax=max(ppois(0:n,lambda=lambda))
plot(x=c(0:n),y=dpois(c(0:n),lambda=lambda),
     ylim=c(0,ymax),xlim=c(-1,n+1),xlab="x",ylab="Función de probabilidad",
     main=paste0(c("Función de probabilidad\n  Po(lambda=",lambda,")"),collapse = ""))
lines(x=rep(0:n,each=2),y=aux,pch=21, type = "h", lty = 2,col="blue")
curve(ppois(x,lambda=lambda),
      xlim=c(-1,n+1),col="blue",ylab="Función de Distribución",
      main=paste0(c("Función de distribución \n Po(lambda=",lambda,")"),collapse = ""))
par(mfrow=c(1,1))

Sea \(X\) un una v.a. \(Po(\lambda=3)\). Entonces

\(P_X(0)=P(X=0), P_X(1)=P(X=1)\) en este orden son

from scipy.stats import poisson
poisson.pmf(0,mu = 3)

## 0.049787068367863944

poisson.pmf(1,mu = 3)

## 0.14936120510359185

Sea \(X\) un una v.a. \(Po(\lambda=3)\). Entonces

\(F_X(0)=P(X\leq 0), F_X(1)=P(X\leq 1)\) en este orden son

poisson.cdf(0,mu = 3)

## 0.04978706836786395

poisson.cdf(1,mu = 3)

## 0.1991482734714558

poisson.pmf(0,mu = 3)+poisson.pmf(1,mu= 3) ## es igual a poisson.cdf(1,lambda=3)

## 0.1991482734714558

Por ejemplo podemos comprobar que \(F_X(10)=\displaystyle\sum_{0}^{10} P_X(x)\)

range(0,10)

## range(0, 10)

poisson.pmf(range(0,10),mu=3)

## array([0.04978707, 0.14936121, 0.22404181, 0.22404181, 0.16803136,
##        0.10081881, 0.05040941, 0.02160403, 0.00810151, 0.0027005 ])

sum(poisson.pmf(range(0,10),mu=3))

## 0.9988975118698846

poisson.cdf(10,mu=3)

## 0.9997076630493527

Como ya hemos visto con scipy.stats podemos pedir los momentos de una variable aleatoria \(Po(3)\)

poisson.stats(mu=3, moments='mv')

## (array(3.), array(3.))

Y también generar secuencias de observaciones aleatorias de una población \(Po(3)\)

poisson.rvs(mu=3,size=40)

## array([1, 3, 2, 5, 3, 2, 0, 3, 5, 1, 2, 2, 2, 3, 6, 4, 0, 2, 6, 3, 4, 0,
##        1, 1, 5, 3, 5, 2, 2, 3, 6, 2, 5, 4, 3, 3, 3, 1, 3, 5])

Supongamos que disponemos de una urna de de sorteos que contiene \(m\) bolas blancas y \(n\) bolas rojas.

En total en esta urna hay \(m+n\) bolas, \(m\) blancas y \(n\) rojas. Si extraemos dos bolas de la urna lo podemos hacer de dos formas:

• Extraer una anotar su color y reponerla. Sacar otra y anotar su color. Hemos extraído la bola con reposición.
• Extraer simultáneamente dos bolas (sin reposición) y contar el número de bolas blancas.

Sea \(X\) es la v.a. que cuenta el número de bolas blancas extraídas.

• En el primer caso, \(X\) es una \(B(n=2,p=\frac{m}{m+n})\) ya que consiste en repetir dos veces el mismo experimento de Bernoulli.
• En el segundo caso, \(X\) sigue una distribución hipergeométrica que estudiaremos en esta sección.

Distribución hipergeométrica

Sean \(n\), \(m\) y \(k\) tres número enteros positivos y tales que \(k<m+n\).

Consideremos una urna que contiene \(m+n\) bolas de las que \(m\) son blancas y las restantes \(n\) no (son no blancas).

El número total de bolas es \(m+n\). Extraemos de forma aleatoria \(k\) bolas de la urna sin reemplazarlas.

Sea \(X\) la v.a. que cuenta el número de bolas blancas extraídas. Diremos que la distribución de \(X\) es hipergeométrica de parámetros \(m\), \(n\) y \(k\) y la denotaremos por \(H(m,n,k)\).

Su dominio es

\[D_X=\left\{x\in\mathbf{N}\mid \max\{0,k-n\}\leq x \leq \min\{m,k\}\right\}\]

Para explicarlo, veamos varios ejemplos:

• \(H(m=5,n=2,k=3)\). Tenemos \(m=5\) bolas blancas, \(n=2\) no blancas y sacamos \(k=3\) bolas sin reposición.
  • En este caso el mínimo de bolas blancas extraídas es \(1=k-n=3-2\), ya que sólo hay dos no blancas.
  • En cambio, el máximo si es \(k=3\), ya que tenemos bolas blancas de “sobra”.

\[D_X=\left\{x\in\mathbf{N}\mid \max\{0,k-n\}\leq x \leq \min\{m,k\}\right\}\]

• \(H(m=2,n=5,k=3)\). Tenemos \(m=2\) bolas blancas, \(n=5\) no blancas y sacamos \(k=3\) bolas sin reposición.
  • En este caso el mínimo de bolas blancas es \(0\) ya que puedo sacar 3 no blancas.
  • En cambio, el máximo si es \(m=2\), ya que aunque saquemos \(k=3\) bolas, al llegar a 2 ya hemos extraído todas las bolas blancas de la urna.
• \(H(m=10,n=10,k=3)\). Tenemos \(m=10\) bolas blancas, \(n=10\) no blancas y sacamos \(k=3\) bolas sin reposición.
  • En este caso podemos obtener desde \(0\) blancas hasta \(k=3\) blancas.

Su función de probabilidad es:

Su función de probabilidad es: \[ P_{X}(x)=\left\{ \begin{array}{ll} \frac{\binom{m}{x}\cdot \binom{n}{k-x}}{\binom{m+n}{k}}, & \mbox{ si } \max\{0,k-n\}\leq x \leq \min\{m,k\}, \mbox { para } x\in \mathbf{N},\\ 0, & \mbox{en otro caso.}\end{array}\right. \]

Observación: otras parametrizaciones

En ocasiones se parametriza una v.a. hipergeométrica mediante \(N=m+n\), número total de bolas, \(k\), número de extracciones y \(p\), probabilidad de extraer una bola blanca.

Así podemos parametrizar alternativamente la distribución hipergeométrica así

\[H(N,k,p)\mbox{ donde } p=\frac{m}{N}.\]

Sea \(X\) una v.a. \(H(m,n,k)\). La función de R para calcular la función de probabilidad en un valor \(x\), \(P(X=x)\), es dhyper(x,m,n,k) y para calcular la función de distribución en un valor \(q\), \(P(X\leq q)\), es phyper(q,m,n,k). Para generar una muestra de valores que siga la distribución \(H(m,n,k)\), hay que usar la función rhyper(nn,m,n,k) donde nn es el número de observaciones aleatorias deseado de la muestra.

Por ejemplo, si \(X\) es una \(H(m=15,n=10,k=3)\), los valores de \(P(X=2)\) y que \(P(X>1)=1-P(X\leq 1)\) son:

dhyper(x=2,m=15,10,k=3)

## [1] 0.4565217

phyper(q=1,m=15,n=10,k=3)# sí, le han puesto q ya veremos el porqué

## [1] 0.3456522

1-phyper(q=1,m=15,n=10,k=3)

## [1] 0.6543478

Una muestra aleatoria de este experimento de tamaño 200 sería:

rhyper(nn=200,m=15,n=10,k=3)

##   [1] 2 3 1 3 1 2 2 3 2 2 1 2 1 2 2 3 3 1 1 1 1 0 2 3 2 1 3 2 2 2 2 3 2 3 3 2 0
##  [38] 1 2 1 3 2 2 3 2 3 2 2 3 2 3 1 2 2 2 2 3 2 2 1 3 2 2 3 1 2 2 2 2 2 3 0 2 0
##  [75] 3 2 2 2 1 2 2 3 1 1 1 2 2 2 2 1 1 3 2 2 3 2 2 1 1 1 3 3 2 2 2 1 3 2 2 2 1
## [112] 1 2 3 2 2 1 2 2 2 2 2 2 3 1 2 3 3 1 1 2 2 1 1 3 2 1 1 2 2 3 1 1 1 2 1 1 3
## [149] 1 2 2 3 3 2 3 1 2 1 2 2 2 1 2 3 1 3 3 3 2 2 1 3 3 1 1 2 2 2 2 2 3 2 1 2 1
## [186] 1 1 1 2 1 1 2 2 2 2 3 3 1 0 2

Sea \(X\) una \(H(m,n,k)\), las funciones de scipy.stats cambian los parámetros

• \(M\) es el número total de bolas. Con nuestra parametrización \(M=m+n\).
• \(n\) es el número de bolas blancas. Con nuestra parametrización \(n=m\).
• \(N\) es el número de extracciones. Con nuestra parametrización \(N=k\).

from scipy.stats import hypergeom

hypergeom.pmf(1,M=15+10,n=15,N=3)

## 0.29347826086956585

hypergeom.cdf(1,M=15+10,n=15,N=3)

## 0.3456521739130442

1-hypergeom.cdf(1,M=15+10,n=15,N=3)

## 0.6543478260869557

Una muestra aleatoria de este experimento sería…

hypergeom.rvs(M=15+10,n=15,N=3,size=100)

## array([1, 2, 2, 2, 3, 1, 3, 1, 2, 3, 3, 2, 2, 2, 3, 2, 2, 2, 2, 3, 2, 3,
##        2, 1, 3, 3, 2, 2, 3, 2, 2, 3, 3, 3, 2, 2, 3, 0, 2, 3, 1, 3, 3, 0,
##        2, 3, 2, 2, 1, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 1, 1, 3, 2, 3, 0,
##        1, 0, 1, 2, 3, 3, 2, 3, 2, 1, 2, 2, 3, 3, 1, 3, 0, 1, 3, 3, 2, 1,
##        2, 3, 2, 1, 2, 2, 2, 2, 2, 3, 1, 2])

from scipy.stats import hypergeom
[M, n, N] = [20, 7, 12] ##20 elementos, 7 del tipo, extraemos 12
x = np.arange(max(0, N-M+n),min(n, N))
fig =plt.figure(figsize=(5, 2.7))
 =ax = fig.add_subplot(1,2,1)
 =ax.plot(x, hypergeom.pmf(x, M, n, N), 'bo', ms=5, label='hypergeom pmf')
 =ax.vlines(x, 0, hypergeom.pmf(x, M, n, N), colors='b', lw=2, alpha=0.5)
 =ax.set_ylim([0, max(hypergeom.pmf(x, M, n, N))*1.1])
for tick in ax.xaxis.get_major_ticks():
   =tick.label.set_fontsize(5)
for tick in ax.yaxis.get_major_ticks():
  =tick.label.set_fontsize(5) 
ax = fig.add_subplot(1,2,2)
 =ax.plot(x, hypergeom.cdf(x, M, n, N), 'bo', ms=5, label='hypergeom cdf')
 =ax.vlines(x, 0, hypergeom.cdf(x, M, n, N), colors='b', lw=2, alpha=0.5)
for tick in ax.xaxis.get_major_ticks():
   =tick.label.set_fontsize(5)
for tick in ax.yaxis.get_major_ticks():
   =tick.label.set_fontsize(5)
 =fig.suptitle('Distribucion Hipergeometrica')
 =plt.show()