Esperanza de la suma de dados de probabilidades variables

Se lanza un dado n veces, luego dos dados n veces y así sucesivamente hasta d dados n veces. En cada caso se calcula la media de la suma de los valores obtenidos (siendo la probabilidad de aparición de cada una de las caras proporcional a un valor aleatorio de una distribución de Poisson de parámetro variable entre 1 y s). El objetivo es ver cómo aumenta la esperanza de la suma a medida que se van lanzando más dados, para los casos 1, 2, ... s. Parece que, a medida que aumenta el parámetro de Poisson, el incremento (variación de la esperanza matemática cada vez que se incrementa en una unidad el número de dados lanzados) se estabiliza alrededor de 3,5.

#INICIO -------------------------

rm(list=ls(all=TRUE))

dado<-1:6

n<-1000 #nº de veces que se lanza cada dado

d<-100 #d es el número máximo de dados lanzados

#Establecemos una semilla para que los resultados de las simulaciones sean reproducibles

set.seed(42)

s<-1000 #s es el nº de casos que estudiamos con dados de diferentes probabilidades

#Inicializamos los vectores A y B
A<-rep(0,s);B<-rep(0,s)

for (j in 1:s)

{

#Las probabilidades de las 6 caras del dado son proporcionales a los elementos (pesos) del vector prdado

prdado<-rpois(6,j)

M<-matrix(rep(0,n*d),nrow=n)

for (i in 1:n) 

{

M[i,]<-sample(dado,d,replace=T,prob=prdado)

M[i,]<-cumsum(M[i,])

}

x<-apply(M,2,mean);y<-diff(x)

A[j]<-mean(y);B[j]<-var(y)

}

#Los valores de la matriz A son las medias de los incrementos de las sumas cada vez que se lanza un dado más
#En la matriz B aparecen las correspondientes varianzas

A;B

plot(1:s,A,type='l',lwd=2)

abline(h=3.5,lwd=3,col='red')

#FIN -------------------------