Simulación convergencia en Probabilidad

Consigna
¿Qué funciones emplearemos y cómo lo haremos?
Convergencia de \(X_n\)
Convergencia de \(Y_n\)
Convergencia de \(Y_n / X_n\)

Consigna

Sea \(X_1, X_2, \dots\) una sucesión de variables aleatorias independientes distribuidas \(U(0 ,1)\). Pruebe que:

\[ \dfrac{\sum_{i=1}^n X_i^2}{\sum_{i=1}^n X_i} \xrightarrow{p} \dfrac{2}{3} \]

Definiendo las sucesiones \(X_n\) y \(Y_n\) como: \[ X_n = \dfrac{\sum_{i=1}^n X_i}{n} \qquad Y_n = \dfrac{\sum_{i=1}^n X_i^2}{n} \]

Se puede demostrar (en base a la LDGN) que: \[ X_n \xrightarrow{p} E(X_1) = \dfrac{1}{2} \qquad Y_n \xrightarrow{p} E(Y_1) = E(X_1^2) = \dfrac{1}{3}\]

Por propiedades, \[ \dfrac{Y_n }{X_n} \xrightarrow{p} \dfrac{1/3}{1/2} = \dfrac{2}{3}\]

A continuación generaremos muestras aleatorias de distintos tamaños provenientes de una distribución uniforme \((0, 1)\), para ver si efectivamente se presentan las convergencias en probabilidad.

¿Qué funciones emplearemos y cómo lo haremos?

Funciones que permiten tomar muestras aleatorias de distribuciones de probabilidad conocidas (uniforme, normal, poisson)

runif(10)  # Muestra de tamaño n=10, de una distribución uniforme (0, 1). Estos parámetros pueden cambiarse.

##  [1] 0.52898787 0.74404879 0.08569978 0.26545157 0.56062061 0.78822620
##  [7] 0.02930056 0.15437680 0.46361223 0.07162891

rnorm(10)  # Muestra de tamaño n=10, de una distribución Normal (0, 1). Estos parámetros pueden cambiarse.

##  [1] -0.9544586 -0.3019802  0.4265645 -0.5220278 -0.9462120 -1.2528616
##  [7]  1.1065187  0.8706713  0.2595833 -0.3557741

rpois(10, lambda = 5)  # Muestra de tamaño n=10, de una distribución poisson con lambda = 5 (en este caso, la función no asume una valor del parámetro or defecto, por lo que hay que incorporarlo).

##  [1] 3 5 5 5 4 5 2 6 4 3

Simulamos \(m\) muestras de tamaño \(n\) provenientes de una distribución uniforme. Dispondremos las realizaciones en una matriz \(X\), donde cada fila representa una muestra, y en la cantidad de columnas está dada por el tamaño de las muestras (\(n\)).

m <- 100 # cantidad de muestras a tomar
n <- 10  # cantidad de muestras a tomar

# Matriz donde 'guardaremos' las muestras
X <- matrix(runif(n*m), nrow = m, ncol = n)

Podemos calcular la media de una muestra (fila) en particular, o de todas las muestras (todas las filas) simultáneamente con la función apply

# media de la muestra '1'
mean(X[1, ])

## [1] 0.542163

# Calculamos la media de cada una de las cada muestras de tamanio n
apply(X, 1, mean)

##   [1] 0.5421630 0.4658526 0.4744550 0.4143791 0.5488333 0.4673041 0.4979160
##   [8] 0.4442514 0.5043181 0.3641588 0.6312693 0.6270850 0.3974034 0.4133549
##  [15] 0.5471144 0.4063190 0.3952005 0.5217043 0.4063593 0.6145535 0.4292541
##  [22] 0.4827278 0.6521018 0.5525213 0.5534718 0.4994029 0.4381122 0.6534748
##  [29] 0.5617021 0.4770255 0.6347221 0.5108623 0.5364520 0.5196190 0.5745253
##  [36] 0.4807536 0.4055315 0.3809835 0.5190311 0.5856824 0.4314840 0.4541317
##  [43] 0.5127003 0.7378451 0.4158828 0.6154631 0.5948265 0.5395949 0.5067913
##  [50] 0.4949260 0.4857456 0.3675386 0.3850265 0.4241930 0.4572151 0.4728636
##  [57] 0.3929693 0.5644671 0.5251211 0.5126370 0.3655075 0.4059024 0.4822333
##  [64] 0.5392568 0.4132684 0.5038262 0.5306348 0.4119788 0.6280295 0.5746457
##  [71] 0.4720785 0.4536338 0.4138830 0.4630058 0.4603245 0.5422760 0.5477354
##  [78] 0.4988872 0.4272931 0.4195251 0.5756989 0.6291569 0.5256334 0.3628874
##  [85] 0.4585950 0.3930238 0.5989263 0.3305258 0.5010114 0.3945380 0.5103372
##  [92] 0.5122826 0.4738913 0.3482415 0.4713247 0.6090851 0.4952459 0.5636527
##  [99] 0.3942912 0.6497366

Como nuestro objetivo calcular la media de muestras de distintos tamaños, generaremos la siguiente función para poder hacerlo de manera más eficiente.

## Creamos una función que nos permitirá tomar muestras de determinado tamaño e inmediatamente calcular las medias
mean_samples <- function(n = 10){
     X = matrix(runif(n*m), nrow = m, ncol = n)
     return(apply(X, 1, mean))
     }

Convergencia de \(X_n\)

Tener en cuenta que cuando aplicamos la función apply, estamos obteniendo la siguiente expresión \[ X_n = \dfrac{\sum_{i=1}^n X_i}{n} \]

# generamos una matriz (por ahora vacía) donde 'guardaremos' las medias de las distintas muestras de cada uno de los tamaño definidos (en este caso, las 20 columnas implica que trabajaremos con 20 tamaños de muestras distintos) 
B <- matrix(NA, m, 20)  # (m filas por 20 columnas)

# aplicamos la funcion tomando muestra de tamanio multiplos de 10 (de 10 a 200), incorporándolas en cada columna de la matriz B
for(i in 1:20){
   B[, i] = mean_samples(i*10)
}

# cambiamos nombres de columnas asignando el tamaño de muestra considerado en cada caso
colnames(B) <- as.character(seq(10, 200, by = 10))
B[1:5, ]  # pedimos que "imprima" las primeros 5 filas de la matriz¨

##             10        20        30        40        50        60        70
## [1,] 0.5408227 0.4744287 0.5357460 0.4739029 0.4856007 0.5077449 0.5296567
## [2,] 0.4016694 0.4561577 0.4621632 0.5678537 0.4624273 0.5144947 0.4634063
## [3,] 0.4788536 0.4915172 0.5346564 0.4348796 0.4616510 0.4936077 0.4426194
## [4,] 0.4999822 0.4228648 0.5371678 0.5152422 0.5558272 0.4626386 0.5363853
## [5,] 0.6440038 0.4252950 0.5311732 0.4890920 0.5497460 0.5394863 0.4566588
##             80        90       100       110       120       130       140
## [1,] 0.4937269 0.4561497 0.4740208 0.4649205 0.5017377 0.4793454 0.4750540
## [2,] 0.4818711 0.5194483 0.4625957 0.4670672 0.5053425 0.5435783 0.4631577
## [3,] 0.4975883 0.5564486 0.5279661 0.4908409 0.5424980 0.4841970 0.5187819
## [4,] 0.5187226 0.4871835 0.5410851 0.4797304 0.4752718 0.5501181 0.4938288
## [5,] 0.4308495 0.4153100 0.4853791 0.5044069 0.5158308 0.4993881 0.4796981
##            150       160       170       180       190       200
## [1,] 0.5174869 0.4886395 0.4834059 0.5013229 0.5004346 0.4835970
## [2,] 0.4873156 0.5262874 0.4684504 0.4936174 0.5037208 0.5086636
## [3,] 0.5054434 0.4950989 0.5150393 0.5171002 0.5205278 0.4843264
## [4,] 0.5214524 0.4957667 0.5118778 0.4807865 0.5290852 0.4839215
## [5,] 0.4857198 0.4954199 0.5192821 0.4626065 0.4877617 0.4741605

Realizando box plots para cada tamaño de muestra, podemos ver que la varianza se reduce a medida que \(n\) crece, y por lo tanto \(X_n\) se encuentra “cerca” de \(1/2\) con probabilidad muy alta.

# generamos boxplot
boxplot(B, xlab = 'n', ylab = 'mean Xn')
abline(h = 1/2, col = "red", lwd = 2)  # incorporamos una línea que representa al valor a donde converge Xn en probabilidad

A continuación pocederemos de forma análoga para analizar la convergencia de \(Y_n\)

Convergencia de \(Y_n\)

Tener en cuenta que cuando aplicamos la función apply, estamos obteniendo la siguiente expresión \[ Y_n = \dfrac{\sum_{i=1}^n X_i^2}{n} \]

Para proceder debemos modificar la fucnión que calcule las medias sobre los cuadrados de las realizaciones de las variables aleatorias.

## Creamos una función que nos permitirá tomar muestras de determinado tamaño e inmediatamente calcular las medias
mean_samples <- function(n = 10){
     X = matrix(runif(n*m), nrow = m, ncol = n)
     return(apply(X^2, 1, mean))  # X^2 eleva cada valor de las muestras al cuadrado
     }

El resto del código es análogo

# generamos una matriz (por ahora vacía) donde 'guardaremos' las medias de las distintas muestras de cada uno de los tamaño definidos (en este caso, las 20 columnas implica que trabajaremos con 20 tamaños de muestras distintos) 
B <- matrix(NA, m, 20)  # (m filas por 20 columnas)

# aplicamos la funcion tomando muestra de tamanio multiplos de 10 (de 10 a 200), incorporándolas en cada columna de la matriz B
for(i in 1:20){
   B[, i] = mean_samples(i*10)
}

# cambiamos nombres de columnas asignando el tamaño de muestra considerado en cada caso
colnames(B) <- as.character(seq(10, 200, by = 10))
B[1:5, ]  # pedimos que "imprima" las primeros 5 filas de la matriz¨

##             10        20        30        40        50        60        70
## [1,] 0.4663185 0.3512836 0.3265745 0.3154977 0.3005180 0.2781387 0.3208049
## [2,] 0.1948251 0.3463984 0.3073680 0.3188864 0.3401818 0.4426400 0.3877116
## [3,] 0.3254424 0.3804446 0.3448004 0.3575052 0.3598972 0.3401637 0.3025192
## [4,] 0.3197110 0.3319935 0.3237937 0.3479445 0.3634923 0.3471807 0.3363946
## [5,] 0.2931878 0.3127539 0.3223880 0.3009040 0.3380005 0.4015190 0.3354067
##             80        90       100       110       120       130       140
## [1,] 0.2929952 0.3146176 0.3177320 0.3689642 0.3592090 0.3274235 0.3364902
## [2,] 0.3548219 0.2731926 0.3540486 0.3105725 0.2773328 0.2983826 0.3420231
## [3,] 0.2980701 0.3531487 0.4015017 0.3184492 0.3025360 0.3067532 0.3333816
## [4,] 0.2616330 0.3377635 0.3269771 0.3697253 0.3373079 0.3642094 0.2860134
## [5,] 0.3319824 0.3136484 0.3517903 0.3112413 0.3393840 0.3227350 0.3218281
##            150       160       170       180       190       200
## [1,] 0.3339262 0.3253702 0.3172511 0.3535975 0.3317104 0.3382012
## [2,] 0.3406054 0.2937628 0.3298657 0.3567397 0.3649143 0.3415607
## [3,] 0.3543865 0.3438667 0.3586830 0.3509904 0.3470578 0.3508280
## [4,] 0.3674478 0.3386235 0.2921754 0.3455451 0.3446752 0.3705752
## [5,] 0.3357418 0.2629140 0.3036886 0.3265668 0.3162326 0.3473090

Realizando box plots para cada tamaño de muestra, podemos ver que la varianza se reduce a medida que \(n\) crece, y por lo tanto \(Y_n\) se encuentra “cerca” de \(1/3\) con probabilidad muy alta.

# generamos boxplot
boxplot(B, xlab = 'n', ylab = 'mean Yn')
abline(h = 1/3, col = "red", lwd = 2)  # incorporamos una línea que representa al valor a donde converge Xn en probabilidad

Convergencia de \(Y_n / X_n\)

Finalmente, corroboraremos mediante las simulaciones que: \[ \dfrac{Y_n }{X_n} \xrightarrow{p} \dfrac{1/3}{1/2} = \dfrac{2}{3}\]

Nuevamente, lo único que adaptaremos es la función a efectos de que para cada muestra obtengamos: \[ \dfrac{\sum_{i=1}^n X_i^2}{\sum_{i=1}^n X_i} \]

mean_samples <- function(n = n){
  X = matrix(runif(n*m), nrow = m, ncol = n)
  X2 = X^2  # elevamos al cuadrado las realizaciones de las distintas muestras
  return(apply(X^2, 1, mean) / apply(X, 1, mean))  # aquí se realiza la modificación mencionada para analizar la convergencia buscada
}