Un ejercicio sobre extracción de bolas con reemplazamiento

Se calcula una probabilidad mediante simulación y se compara con el valor exacto. Para conocer cómo puede obtenerse este último valor puede verse el ejemplo 4.20 del interesante texto Dados, Monedas y Urnas (pág. 134) de Víctor Hernández y Ricardo Vélez (UNED, 1ª ed.).

#INICIO -------------------------
rm(list=ls(all=TRUE))

#Se extraen n bolas con reemplazamiento de una urna que contiene 3 bolas marcadas con los números 1, 2 y 3.
#Calcular la probabilidad de que en k ocasiones, 0<=k<=n-1, el número de una bola sea mayor que el de la primera.

bola<-1:3

#Se eligen los valores de n y k
n<-8;k<-3

#Definimos el número de veces N que simulamos la experiencia
N<-1000000

#Realizamos las extracciones
x<-replicate(N,
{
a<-sample(bola,n,replace=T);
b<-a-a[1];
length(b[b>0])
}
)

#Calculamos el valor aproximado VA obtenido por simulación
VA<-length(x[x==k])/length(x)

#Cálculo del valor exacto VE (expresiones diferentes para k=0 y para k>0)
if(k==0) {
VE<-(1+2^(n-1)+3^(n-1))/3^n
} else {
VE<-choose(n-1,k)*(2^k+2^(n-1-k))/3^n
}

paste("El","valor","aproximado","para k=",k,"es",round(VA,5))
paste("El","valor","exacto","para k=",k,"es",round(VE,5))
#FIN -------------------------