Unidad II Números Pseudo-Aleatorios. Introducción Los números aleatorios son un ingrediente básico en la simulación de la mayoría de los sistemas discretos. La mayoría m ayoría de los lenguajes de computadora tienen una subrutina o función que permite generar un número aleatorio. De manera análoga, los lenguajes de simulación generan números aleatorios que son usados para generar tiempos de evento y otr as variables aleatorias. En este capitulo se utiliza el término pseudo para implicar que el acto mismo de la generación de números aleatorios por un método dado, remueve el potencial para la aleatoriedad verdadera. Si el método es conocido, el conjunto de números aleatorios puede ser repetido. Las secuencias de números pseudoaleatorios no muestran ningún patrón o regularidad aparente desde un punto de vista estadístico, aunque estos hayan sido generados por un algoritmo completamente determinista, en el que las mismas condiciones iniciales producen siempre el mismo resultado. Los mecanismos de generación de números aleatorios que se utilizan en la l a mayoría de los sistemas informáticos son en realidad procesos pseudo-aleatorios. Una de las utilidades principales de los números pseudo-aleatorios se lleva a cabo en el llamado método de Montecarlo, con múltiples utilidades, por ejemplo para hallar áreas / volúmenes encerradas en una gráfica y cuyas integrales son muy difíciles de hallar o irresolubles; mediante la generación de puntos basados en estos números, podemos hacer una buena aproximación de la superficie /volumen total, encerrándolo en un cuadrado / cubo, aunque no lo suficientemente buena. Asimismo, también destacan en el campo de la criptografía. Por ello se sigue investigando en la generación de dichos números, empleando por ejemplo medidores de ruido blanco o analizadores atmosféricos, ya que experimentalmente se ha comprobado que tienen una aleatoriedad bastante alta.
2.1 Métodos de generación de números Pseudoaleatorios Pse udoaleatorios Existen multitud de métodos para la generación de números pseudo-aleatorios que tengan valor entre 0 y 1, cada uno de estos métodos puede considerarse adecuado siempre y cuando el conjunto de números generados cumpla con las pruebas de aleatoriedad, que se verán más adelante. Entre los métodos mas utilizados podemos mencionar: 2.1.1 Métodos No Congruenciales Entre los métodos no congruenciales, se encuentran el de cuadrado de medios, el de producto de medios y el de multiplicador constante.
Algoritmo de cuadrado de medios Este método fue propuesto en la década de los 40’s por John von Neumann y Metropolis [1]. El procedimiento consiste en lo siguiente: 1. Seleccionar un numero numero de inicio (X0), llamado semilla con d dígitos (d > 3) 2. Elevar X0 al cuadrado, completar el número de dígitos con ceros a la izquierda, de manera que sea el doble de d ; sea X1 los d dígitos del centro y sea r i 0.d . 3. Sea Yi el resultado de elevar al cuadrado a Xi (los dígitos del centro) y r i 1 0.d de los dígitos obtenidos. 4. Repetir el paso paso 3 hasta hasta obtener los n números números
r i deseados.
Ejemplo: Generar los primeros 5 números
r i a
partir de una semilla X0 =5735,
observando que se utilizan 4 dígitos. Y 0 (5735) 2 32890225
X 1 8902
r 1 0.8902
Y 1 (8902) 2 79245604 X 2 2456 Y 2 (2456) 2 06031936 X 3 0319
r 2 0.2456 r 3 0.0319
Y 3 (0319)2 00101761 X 4 1017
r 4 0.1017
Y 4 (1017) 2 01034289 X 5 0342
r 5 0.0342
Este algoritmo generalmente es incapaz de generar una secuencia de periodo de vida de n grande.
r i con
Algoritmo de producto de medios Este método es semejante al anterior, la diferencia entre ambos radica en que este último requiere de dos semillas para iniciar. Ambas semillas deben ser del mismo número de dígitos. El procedimiento consiste en lo siguiente: 1. Seleccionar un número número de inicio (X0), llamado semilla con d dígitos (d > 3). 2. Seleccionar otra semilla (X1), con d dígitos (d > 3). 3). 3. Sea Y0 = X0*X1, sea X2 = los dígitos d del centro, completando el numero de dígitos con ceros a la izquierda, de manera que sea el doble de d ; y sea r i 0.d . 4. Sea Yi = Xi*Xi+1; sea Xi+2 (los dígitos del centro) y r i 1 0.d de los dígitos obtenidos. 5. Repetir el paso 4 hasta obtener obtener los n números números r i deseados. Ejemplo:
Generar los primeros 5 números
r i a
partir de las semillas X0 =5015 y X1
=5734, observando que se utilizan 4 dígitos. Y 0 (5015)(5734) 28756010
X 1 7560
r 1 0.7560
Y 1 (5734)(7560) 43349040 X 2 3490 Y 2 (7560)(3490) 26384400 X 3 3844
r 2 0.3490 r 3 0.3844
Y 3 (3490)(3844) 13415560 X 4 4155
r 4 0.4155
Y 4 (3844)(4155) 15971820
r 5 0.9718
X 5 9718
Algoritmo de multiplicador constante Este algoritmo es semejante al anterior, la diferencia ahora radica en que este último considera que una de las semillas permanezca constante. Ambas semillas deben ser del mismo número de dígitos. El procedimiento consiste en lo siguiente: 1. Seleccionar un número de inicio (X0), llamado semilla con d dígitos (d > 3). 2. Seleccionar una constante (K), con d dígitos (d > 3). 3. Sea Y0 = K*X0, sea X1 = los dígitos d del centro, completando el numero de dígitos con ceros a la izquierda, de manera que sea el doble de d ; y sea r i 0.d . 4. Sea Yi = K*Xi; sea Xi+1 (los dígitos d del centro) y r i 1 0.d de los dígitos obtenidos. 5. Repetir el paso 4 hasta obtener los n números r i deseados. Ejemplo: Generar los primeros 5 números r i a partir de las semillas X 0 =9803 y la constante K =6965, observando que se utilizan 4 dígitos para cada una. Y 0 (6965)(9803) 68277895 X 1 2778
r 1 0.2778
Y 1 (6965)(2778) 19348770 X 2 3487 Y 2 (6965)(3487) 24286955 X 3 2869
r 2 0.3487 r 3 0.2869
Y 3 (6965)(2869) 19982585
X 4 9825
r 4 0.9825
Y 4 (6965)(9825) 68431125 X 5 4311
r 5 0.4311
Los mecanismos de generación de números aleatorios que se utilizan en la mayoría de los sistemas informáticos son en realidad procesos pseudoaleatorios.
2.1.2 Métodos Congruenciales
Entre los métodos congruenciales, se encuentran el lineal, el multiplicativo, el aditivo y el congruencial no lineal cuadrático. Algoritmo lineal Este algoritmo fue propuesto en 1951 por D.H. Lehmer [2]. Se considera que es el algoritmo mas ampliamente utilizado para generar números aleatorios. Muchos otros métodos han sido propuestos y son revisados por Bratley, Fox and Schrage [3], Law and Kelton [4] y por Ripley [5]. Este método produce una secuencia de enteros X 1, X2, … entre cero y m-1, de acuerdo a la siguiente relación recursiva: X i 1 (aX i c) mod m, i 0,1,2,...
El valor inicial X 0 es llamado la semilla, a es la constante multiplicadora, c es el incremento y m es el módulo. Para lograr que el algoritmo tenga un largo periodo de vida n, es necesario que estos parámetros cumplan con lo siguiente: m 2 g ,
a 1 4k ,
k y g son enteros y c primo relativo a m
Es importante señalar que este algoritmo recursivo genera una secuencia de números enteros {0, 1, 2, 3, … , m-1}, y para obtener los números pseudoaleatorios en el intervalo (0, 1) se requiere hacer: r i
X i m 1
, i 0,1,2,..., n
Para comprender mejor la mecánica de este algoritmo, el procedimiento consiste en lo siguiente: 1. Seleccionar un número de inicio (X0), llamado semilla con d dígitos (d 1). 2. Sea Y0 = (a*X0+c) mod m . 3. Sea Yi = (a*Xi+c) mod m; sea Xi+1 = Yi y r i 1 X i /(m 1) . 4. Repetir el paso 3 hasta obtener los n números
>
r i deseados.
Ejemplo: Generar los primeros 5 números
r i a
partir de las semillas X 0 =37, a=17, c =33,
y m=64. X 1 (17 * 37 33) mod 64 22
r 1 22 / 63 0.3492
X 2 (17 * 22 33) mod 64 23 X 3 (17 * 23 33) mod 64 40
r 2 23 / 63 0.3651 r 3 40 / 63 0.6349
X 4 (17 * 40 33) mod 64 9
r 4 9 / 63 0.1429
X 5 (17 * 9 33) mod 64 58
r 5 58 / 63 0.9206
Algoritmo Multiplicativo Es el mismo que el algoritmo congruencial lineal, solo que con c =0, de manera que la ecuación recursiva se reduce a: X i 1 (aX i ) mod m, i 0,1,2,..., n
A diferencia del algoritmo congruencial lineal, este tiene la ventaja que se puede reducir una operación en el procedimiento, por lo que la secuencia es la misma. Los parámetros de arranque son X 0, a y m y deben cumplir con las condiciones siguientes: m 2 g ,
a 3 8k ,
k y g son enteros y X 0 debe ser impar.
Ejemplo: Generar los primeros 5 números
r i a
partir de las semillas X 0 =10, a=7, y m=13.
X 1 (7 *10) mod13 5
r 1 5 / 13 0.3846
X 2 (7 * 5) mod13 9 X 3 (7 * 9) mod13 11
r 2 9 / 13 0.6923 r 3 11 / 13 0.8462
X 4 (7 * 11) mod13 12
r 4 12 / 13 0.9231
X 5 (7 *12) mod13 6
r 5 6 / 13 0.4615
Algoritmo Congruencial Aditivo Es semejante al algoritmo congruencial lineal, solo que utiliza una secuencia previa de n números enteros x 1, x 2, x 3, …, x n, como entrada, de manera que la nueva secuencia de números generada x n+1, x n+2 , x n+3, …, y la ecuación recursiva utilizada es: X i ( X i 1 X i n ) mod m, i n 1, n 2, n 3,...
Los números r i generados se obtienen: r i
X i m 1
Ejemplo: Generar los 10 números pseudo-aleatorios r i a partir de las semillas siguientes X 1 =65, X 2 =89, X Para generarlos se sigue el 3 =98, X 4 =03, X 5 =69. procedimiento:
X 6 ( X 5 X 1 ) mod100 (69 65) mod100 34
r 1 34 / 99 0.3434
X 7 ( X 6 X 2 ) mod100 (34 89) mod100 23 X 8 ( X 7 X 3 ) mod100 (23 98) mod100 21
r 2 23 / 99 0.2323 r 3 21 / 99 0.2121
X 9 ( X 8 X 4 ) mod100 (21 3) mod100 24
r 4 24 / 99 0.2424
X 10 ( X 9 X 5 ) mod100 (24 69) mod100 93
r 5 93 / 99 0.9393
X 11 ( X 10 X 6 ) mod100 (93 34) mod100 27
r 1 27 / 99 0.2727
X 12 ( X 11 X 7 ) mod100 (27 23) mod100 50 X 13 ( X 12 X 8 ) mod100 (50 21) mod100 71
r 2 50 / 99 0.5050 r 3 71 / 99 0.7171
X 14 ( X 13 X 9 ) mod100 (71 24) mod100 95
r 4 95 / 99 0.9595
X 15 ( X 14 X 10 ) mod100 (95 93) mod100 88
r 5 93 / 99 0.8888
Algoritmo No lineal Cuadrático Es semejante al algoritmo congruencial lineal, solo que se tiene una constante más por el término cuadrático, la del término lineal y el t érmino independiente. La ecuación recursiva se define como: X i 1 (aX i2 bX i c) mod m, i 0,1,2,..., n
A diferencia del algoritmo congruencial lineal, este tiene la necesidad de incrementar una operación en el procedimiento, pero la secuencia es semejante. Los parámetros de arranque son: X 0, a, b, c y m y deben cumplir con las condiciones siguientes: a numero _ par , m 2 g , entero y X 0 debe ser impar.
(b a) mod 4 1,
c numero _ impar ,
g es
Ejemplo: Generar los primeros 5 números r i a partir de la semilla X 0 =27, y las constantes a=12, b=13, c =43, y m=64. X 1 (12 * 27 2 13 * 27 43) mod 64 54
r 1 54 / 63 0.8571
X 2 (12 * 27 2 13 * 27 43) mod 64 25 X 3 (12 * 272 13 * 27 43) mod 64 60
r 2 25 / 63 0.3968 r 3 60 / 63 0.9524
X 4 (12 * 272 13 * 27 43) mod 64 55
r 4 55 / 63 0.8730
X 5 (12 * 27 2 13 * 27 43) mod 64 2
r 5 2 / 63 0.0317
2.2 Pruebas estadísticas de aleatoriedad. Una secuencia de números aleatorios, R 1, R2,…, debe tener dos propiedades estadísticas importantes de manera de validar si los números obtenidos son aptos para utilizarse en un estudio de simulación. Las pruebas que deben realizarse para esta validación, se agrupan en las siguientes:
1. Prueba de media y varianza. 2. Pruebas de uniformidad. 3. Pruebas de independencia
2.2.1 Prueba de media y varianza Cada numero aleatorio R i es una muestra independiente extraída de una distribución continua uniforme entre cero y uno. Esto es, la función de densidad de probabilidad (pdf) esta dada por:
1, 0 x 1 f ( x) 0, 0 x 1 La pdf de la función es mostrada en la siguiente figura:
El valor esperado de Ri esta dado por: 1 1 x2 1 E ( R) xdx 0 2 2 0
Y para la varianza está dado por: 1
1
2
x3
V ( R) E ( R 2 ) E ( R) x 2 dx 0 3 2 2
1
0
1 4
1 3
1 4
1 12
La prueba de media consiste en determinar el promedio de los n números que contiene el conjunto R , mediante la definición, es decir: R
1
n
R n
i
i 1
Posteriormente se calculan los límites de aceptación inferior y superior con las ecuaciones siguientes: 1 1 1 1 LI R z / 2 LS R z / 2 2 2 12n 12n Si el valor promedio de R se encuentra entre los limites de aceptación, concluimos que no se puede rechazar que el conjunto R , tiene un valor esperado de 0.5 con un nivel de aceptación de 1- α. En caso contrario se rechaza el conjunto R . Para el cálculo de los limites de aceptación se utiliza el estadístico Z α/2 , el cual se determina por medio de la tabla de distribución normal estándar (también se puede usar las funciones DISTR.NORM.INV de Excel y norminv de Matlab, como se verá en la práctica)
Ejemplo: Considere los 30 números del conjunto R que se presentan a continuación, y determine si tienen un valor esperado de 0.5 con un nivel de aceptación de 95%. 0.587301587 0.031746032 0.746031746 0.698412698 0.904761905 0.349206349
0.365079365 0.317460317 0.523809524 0.984126984 0.682539683 0.634920635
0.142857143 0.603174603 0.301587302 0.253968254 0.46031746 0.920634921
0.349206349 0.047619048 0.396825397 0.206349206 0.619047619 0.571428571
0.841269841 0.285714286 0.555555556 0.952380952 0.777777778 0.222222222
El conjunto R contiene 30 números, por lo tanto n=30. Un nivel de aceptación
del 95% implica que α=5%. Enseguida procedemos a calcular el promedio de los números y los limites de aceptación. R
n
1
1
1
R 30 R 30 (0.587301587 0.365079365 ... 0.571428571 0.222222222) n i
i
i 1
R
30
1 30
i 1
(15.3333333) 0.51111111
Ahora se calculan los límites de aceptación inferior y superior con las ecuaciones siguientes:
1 1 z / 2 z 2 12 n 2
LI R
1
LI R
1
LS R
1
/
2
1 12(30)
1 0.5 0.10330107 0.39669893 (1.96) 2 12 ( 30 ) 1 1 z / 2 z 2 12n 2
1 / 2 12 ( 30 )
1 0.5 0.10330107 0.60330107 (1.96) 2 12(30) Como el valor promedio: R 0.51111111 se encuentra entre los límites de LS R
1
aceptación, se concluye que no se puede rechazar que el conjunto de 30 números R tiene un valor esperado de 0.5 con una aceptación de 95%. La otra propiedad, que se menciono debe satisfacer el conjunto R, es que sus números tengan una varianza de 1/12. La prueba consiste en pr imero determinar la varianza de los n números que contiene el conjunto el conjunto R mediante la siguiente ecuación: 1 n V ( R) ( Ri R)2 n 1 i 1
Después se calculan los límites de aceptación inferior y superior con las ecuaciones siguientes:
2
LI V ( R )
2
/ 2,n1
LS V ( R )
12(n 1)
1 / 2 , n 1
12(n 1)
Si el valor de la varianza del conjunto de números se encuentra entre los limites de aceptación, decimos que no se puede rechazar que el conjunto R tiene una varianza de 1/12 con un nivel de aceptación de 1- α ; de lo contrario, se rechaza que el conjunto r tiene una varianza de 1/12. Los estadísticos correspondientes a la función chi-cuadrada se obtienen de las tablas correspondientes al nivel de significancia y los grados de libertad correspondientes.(En Excel se tiene la función PRUEBA.CHI.INV y en Matlab el valor se obtiene con la función chiinv) Ejemplo: Considere los 30 números del conjunto R y el nivel de significancia que se presentaron en el ejemplo anterior. El conjunto R contiene 30 números, por lo tanto n=30. Un nivel de aceptación
del 95% implica que α=5%. Enseguida procedemos a calcular la varianza de los números y los limites de aceptación correspondientes. 1 n 1 30 2 V ( R) ( Ri R) ( Ri 0.5111111) 2 n 1 i 1 30 1 i 1
V ( R)
30
1
(0.587301587 0.5111111) 29
2
... (0.2222222 - 0.5111111)2
i 1
V ( R) 0.07352499 2
LI V ( R)
1 / 2,n1
LS V ( R )
/ 2, n 1
12(n 1)
2
10.05 / 2, 29
0.05 / 2, 29
2
12(n 1)
12(29)
2
2
12(29)
0.975, 29
12(29)
2
0.025, 29
12(29)
16.0471 348
45.7223 348
0.04733045
0.13138591
Dado que el valor de la varianza esta entre los limites de aceptación, podemos decir que no se puede rechazar que el conjunto de 30 números R tiene una varianza de 1/12=0.08333.
2.2.3 Pruebas de uniformidad Una de las propiedades más importantes que deben satisfacer los números pseudo-aleatorios de un conjunto R es la uniformidad. Algunas de las pruebas que tratan de corroborar si los números del conjunto R son uniformes, es decir, si se pueden considerar aleatorios, son varias, aquí se van a considerar las siguientes: 1. Chi-cuadrada 2. Kolmogorov-Smirnov Prueba Chi-cuadrada
Con esta prueba se busca determinar si los números del conjunto R se encuentran uniformemente distribuidos en el intervalo (0, 1). Para ello es necesario dividir el intervalo en un número s de subintervalos, que se recomienda que sea s N . Una vez hecha la división, se clasifican cada uno de los números del conjunto R en los s subintervalos. La cantidad de números del conjunto R que se encuentran en cada subintervalo se le conoce como frecuencia observada (Oi ), y a la cantidad de números que se espera encontrar en cada subintervalo se le llama frecuencia esperada (E i ); en teoría esta última seria N/s. Con los valores de obtenidos de Oi y E i se determina el estadístico de chi-cuadrada de la muestra con la expresión: 2 s ( E i Oi ) 2 0 E i i 1
Si el valor estadístico de
2
0
es menor al valor en la tabla para
2
, s 1 ,
entonces
no se puede rechazar que el conjunto de números R sigue una distribución uniforme. En caso contrario, se rechaza que el conjunto R sigue una distribución uniforme. Ejemplo: Realizar la prueba chi-cuadrada al siguiente conjunto de R , que contiene 100 numeros. Considerar un nivel de confianza de 95%. 0.950 0.615 0.058 0.015 0.838 0.193 0.497 0.727 0.795 0.137 0.231 0.792 0.353 0.747 0.020 0.682 0.900 0.309 0.957 0.012 0.607 0.922 0.813 0.445 0.681 0.303 0.822 0.838 0.523 0.894 0.486 0.738 0.010 0.932 0.379 0.542 0.645 0.568 0.880 0.199 0.891 0.176 0.139 0.466 0.832 0.151 0.818 0.370 0.173 0.299 0.762 0.406 0.203 0.419 0.503 0.698 0.660 0.703 0.980 0.661 0.456 0.935 0.199 0.846 0.709 0.378 0.342 0.547 0.271 0.284 0.019 0.917 0.604 0.525 0.429 0.860 0.29 0.445 0.252 0.469 0.821 0.410 0.272 0.203 0.305 0.854 0.341 0.695 0.876 0.065 0.445 0.894 0.199 0.672 0.190 0.594 0.534 0.621 0.737 0.988 Para aplicar la prueba de una manera mas clara, se recomienda hacer una tabla como la siguiente: ( E i Oi ) 2 n E i Oi Intervalo E i s [0.00-0.10) 9 10 0.1 [0.10-0.20) 10 10 0.0 [0.20-0.30) 9 10 0.1 [0.30-0.40) 9 10 0.1 [0.40-0.50) 11 10 0.1 [0.50-0.60) 9 10 0.1 [0.60-0.70) 12 10 0.4 [0.70-0.80) 9 10 0.1 [0.80-0.90) 13 10 0.9 [0.90-1.00) 9 10 0.1 Con lo que el estadístico es: s ( E i Oi ) 2 2 0 0.1 0.0 0.1 0.1 0.1 0.1 0.4 0.1 0.9 0.1 2.0 E i i 1
El cual es menor que el estadístico para chi-cuadrada del 95% de aceptación. 2 0.05,9 16.9 , por lo que se puede concluir que no se puede rechazar que los números R siguen una distribución uniforme. Prueba Kolmogorov-Smirnov Con esta otra prueba también se busca determinar si los números del conjunto R se encuentran uniformemente distribuidos en el intervalo (0, 1). Se recomienda aplicar a conjuntos de números R pequeños, n<20 . El procedimiento consiste en lo siguiente: 1º Se ordenan de menor a mayor los números del conjunto R. R1 R2 R3 ... Rn + 2º Determinar los valores D , D- y D, con las siguientes ecuaciones: 1 D max Ri 1 i n n i 1 1 i n n D max D , D
D max Ri
3º Determinar el valor critico de D de acuerdo con los valores críticos de Kolmogorov-Smirnov para un grado de confianza α, y según el tamaño de la muestra n. 4º Si el valor D es mayor que el valor critico Dα,n, se concluye que los números del conjunto R no siguen una distribución uniforme; en caso contrario, se dice que no se ha detectado diferencia significativa entre la distribución de los números del conjunto R y la distribución uniforme. Ejemplo: Realizar la prueba de Kolmogorov-Smirnov, con un nivel de confianza del 90% para un conjunto R de 10 números, que se indica a continuación: R={0.456, 0.935, 0.199, 0.846, 0.709, 0.378, 0.342, 0.547 , 0.271, 0.284} El nivel de confianza del 90% implica α=10%. Ordenando de menor a mayor el conjunto de números R, se tiene: 0.199
0.271
0.284
0.342
0.378
0.456
0.547
0.709
0.846
Para determinar los valores de: , se debe realizar la siguiente tabla: i 1 2 3 4 5 6 7 8 9
0.935
10
i n
Ri
i n i
Ri
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.199
0.271
0.284
0.342
0.378
0.456
0.547
0.709
0.846
0.935
-0.099
-0.071
0.016
0.058
0.022
0.144
0.153
0.091
0.054
0.065
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
0.199
0.171
0.084
0.042
-0.022
-0.044
-0.053
0.009
0.046
0.035
1
n
Ri
0.10
i 1 n
n
D+
D-
D
10
0.153
0.199
0.199
De acuerdo con la tabla de valores para la prueba de Kolmogorov-Smirnov, el valor critico correspondiente a n=10 es D0.10,10 =0.368, que resulta mayor que el valor D=0.199; por lo tanto, se concluye que los números del conjunto R se distribuyen de manera uniforme.
2.2.3 Pruebas de independencia La Otra propiedad importante que deben satisfacer los números de un conjunto R es la independencia. A continuación se comentan algunas de las pruebas que tratan de corroborar si los números del conjunto R son independientes, es decir, si se puedan considerar aleatorios. Las pruebas que se van a considerar son: 3. Corridas arriba y abajo 4. Corridas encima y debajo de la media Estos procedimientos se basan en determinar una secuencia de números (S) que solo contienen unos y ceros de acuerdo con la comparación de los elementos del conjunto R. Prueba de corridas arriba y abajo Para este método, se hace la comparación entre los números r i y r i+1. Posteriormente se determina el numero de corridas observadas, C 0 (una corrida equivale a cada conjunto de ceros o unos consecutivos). Luego se calcula el valor esperado y la varianza del número de corridas y el estadístico Z 0, mediante las ecuaciones: (16n 29) ( 2n 1) 2 C C 0 3 90 C 0 C 0 Z 0
0
C
0
Si el estadístico Z 0 queda fuera del intervalo –z /2 ≤ z0 ≥ z /2, se concluye que los números del conjunto R no son independientes, de lo contrario, no se puede rechazar que el conjunto R es independiente. Para la construcción de la secuencia de unos y ceros, se coloca un cero si el número r i+1 es menor o igual al número anterior r i: en caso de ser mayor que el numero r i anterior, se coloca un uno. A continuación se da una tabla de números aleatorios para encontrar la secuencia de ceros y unos.
0.587301587 0.031746032 0.746031746 0.698412698
0.365079365 0.317460317 0.523809524 0.984126984
0.142857143 0.603174603 0.301587302 0.253968254
0.349206349 0.047619048 0.396825397 0.206349206
0.841269841 0.285714286 0.555555556 0.952380952
Considerando la secuencia de los 20 números del conjunto anterior (utilizando renglón por renglón), la secuencia de unos y ceros es: S={0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1} Debe aclararse que para este método la secuencia S tendrá solo n-1 números, es decir 19. Esto se debe a que el primer número (que es el 0. 587301587) no tiene un número anterior con el cual compararlo. Para contar las corridas, debe
recordarse que una corrida es una secuencia de puros unos o ceros consecutivos, a continuación se marcan las corridas correspondientes a la secuencia S. S={0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1} Contando el número de secuencias subrayadas, se obtiene que las corridas son 10. Ejemplo: Realizar la prueba de corridas arriba y abajo con un nivel de aceptación de 95% para el siguiente conjunto R de 40 números pseudo-aleatorios: 0.5828 0.4235 0.5155 0.3340
0.4329 0.2259 0.5798 0.7604
0.5298 0.6405 0.2091 0.3798
0.7833 0.6808 0.4611 0.5678
0.7942 0.0592 0.6029 0.0503
0.4154 0.3050 0.8744 0.0150
0.7680 0.9708 0.9901 0.7889
0.4387 0.4983 0.2140 0.6435
0.3200 0.9601 0.7266 0.4120
0.7446 0.2679 0.4399 0.9334
Realizando la secuencia de unos y ceros de acuerdo a la comparación con el anterior, que puede ser renglón por renglón o columna por columna, en este caso se sigue la primera. La secuencia S correspondiente usando cada renglón, se obtiene: S={ 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1} Donde el número de corridas es 26 y como la aceptación es del 95%. Los cálculos se presentan a continuación: (2(40) 1) 79 (16(40) 29) 2 C C 6.788 26 . 333 0 0 3 3 90 26 26.333 0.1278 Z 0 6.788 Que resulta estar dentro del intervalo del estadístico para el nivel de rechazo del 5% (=0.05), que corresponde a –1.96= z0.05/2 ≤ Z0 = -0.1278 ≤ z0.05/2 =1.96, por lo que no se puede rechazar los números del conjunto R son independientes. Prueba de corridas encima y debajo de la media En esta prueba se compara a cada elemento del conjunto de números R, con la media para una distribución uniforme (0.5). Posteriormente se determina el numero de corridas observadas, C0, así como la cantidad de ceros obtenidos, n0, así como la cantidad de unos, n1. C0 es el número de corridas, obtenidas de la misma manera que en la prueba anterior, además debe cumplirse que n0 + n1 = N. Luego se calcula el valor esperado y la varianza del número de corridas del conjunto y con las mismas el estadístico Z0 a partir de las ecuaciones siguientes: 2n0 n1 1 2n0n1 (2n0n1 n) 2 C
0
n
C
2
0
Z 0
n2 (n 1)
C 0 C 0 2
C
0
Al igual que en la anterior prueba, considerando un nivel de rechazo , si dicho estadístico queda fuera del intervalo –z /2 ≤ Z0 ≥ z /2, se concluye que los números del conjunto R no son independientes, de lo contrario, no se puede
rechazar que el conjunto R es independiente. La secuencia de unos y ceros se construye asignando un cero si el numero r i es menor o igual que 0.5 y un uno en caso contrario. A continuación se da una tabla de números pseudoaleatorios para encontrar la secuencia de ceros y unos. 0.587301587 0.031746032 0.746031746 0.698412698
0.365079365 0.317460317 0.523809524 0.984126984
0.142857143 0.603174603 0.301587302 0.253968254
0.349206349 0.047619048 0.396825397 0.206349206
0.841269841 0.285714286 0.555555556 0.952380952
Considerando la secuencia de los 20 números del conjunto anterior (utilizando renglón por renglón), la secuencia de unos y ceros es: S={1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1} El número de corridas se determina de la misma manera que en el anterior. En este caso son 11 corridas y se tienen 11 ceros y 9 unos. Ejemplo: Realizar la prueba de corridas encima y debajo de la media con un nivel de aceptación de 90% para el siguiente conjunto R de números aleatorios: q=round(v*10000)/10000 0.5828 0.4235 0.5155 0.3340
0.4329 0.2259 0.5798 0.7604
0.5298 0.6405 0.2091 0.3798
0.7833 0.6808 0.4611 0.5678
0.7942 0.0592 0.6029 0.0503
0.4154 0.3050 0.8744 0.0150
0.7680 0.9708 0.9901 0.7889
0.4387 0.4983 0.2140 0.6435
0.3200 0.9601 0.7266 0.4120
0.7446 0.2679 0.4399 0.9334
Realizando la secuencia de unos y ceros de acuerdo a la comparación con el anterior, que puede ser renglón por renglón o columna por columna, en este caso se sigue la primera. La secuencia S correspondiente usando cada renglón, se obtiene: S={ 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1} Donde el número de corridas es 27, el número de ceros es 19 y el número de unos es 21 y como la aceptación es del 90%. Los cálculos se presentan a continuación: 2(22)(18) 1 792 C 0.5 20.3 0 40 2 40 2
C
0
Z 0
2(22)(18)(2(22)(18) 40) 402 (40 1)
17 20.3
792(752) 1600(39)
595584 62400
9.54
3.03
0.98 3.09 9.54 Que resulta que están dentro del intervalo del estadístico para el nivel de rechazo del 10% (=0.1), que corresponde a –1.645= z0.1/2 ≤ Z0 = 0.98 ≤ z0.1/2 =1.645, por lo que no se puede rechazar los números del conjunto R son independientes. Algunas consecuencias de las propiedades de uniformidad e independencia son las siguientes: 1. Si el intervalo (0, 1) se divide en n clases, o subintervalos de igual magnitud, el número de observaciones esperado en cada intervalo es N/n, donde N es el número total de observaciones.
2. La probabilidad de encontrar un valor en un intervalo particular es independiente del valor anterior determinado. Para asegurar que dichas propiedades deseables sean alcanzadas, un numero de pruebas pueden ser desarrolladas (afortunadamente, las pruebas apropiadas son ya realizadas por la mayoría de los software de simulación comerciales).
2.3 Método de Monte Carlo. El método de Monte Carlo ha sido utilizado por siglos, sin embargo, durante la 2ª Guerra Mundial este fue utilizado para simular los aspectos probabilísticos de la difusión de neutrones (que fue la primera aplicación real). Se le conoce así por la capital de Mónaco (la capital mundial de los casinos), por el uso de los juegos de azar. 2.3.1 Características Los métodos que no son de Monte Carlo involucran ecuaciones ODE/PDE que describen al sistema. Los métodos de Monte Carlo son técnicas estocásticas como los juegos de azar. Están basados en el uso de números aleatorios y estadísticas de probabilidad para simular problemas. Algunas veces se puede llamar Método de Monte Carlo si se utilizan números aleatorios para examinar el problema que se está resolviendo. Por ejemplo, la resolución de ecuaciones que describen las interacciones de dos átomos. Esto sería factible sin necesidad de utilizar el método de Monte Carlo. Sin embargo, la resolución de las interacciones de miles de átomos que utilizan las mismas ecuaciones es imposible. Sin embargo, las soluciones son imprecisas y puede ser muy lenta si se desea una mayor precisión. Función de Densidad de Probabilidad En primer lugar, tendríamos que determinar la función de densidad de probabilidad (Probability Density Function o PDF). A continuación, realizar un muestreo aleatorio del PDF. Mantenemos registro de cada simulación realizada y de las cuentas correspondientes. Una función de densidad de probabilidad (o la función de distribución de probabilidad) es una función f definida en un intervalo (a, b) y que tiene las siguientes propiedades: Los Componentes de Simulación de Monte Carlo Funciones de Distribución de Probabilidad (pdf) - el sistema físico (o matemático) debe ser descrito por un conjunto de funciones PDF. Generador de Números Aleatorios - un conjunto de números aleatorios distribuidos uniformemente en el intervalo de la unidad debe estar disponible. Regla de Muestreo - una procedimiento para el muestreo de las funciones pdf especificadas, asumiendo que se debe dar la disponibilidad de números aleatorios en el intervalo unidad. Anotando (o Contando) - los resultados deben ser acumulados en recuentos totales o calificaciones de las cantidades de interés.
Estimación de Error - debe determinarse una estimación del error estadístico (varianza) como una función del número de pruebas y otras cantidades. Técnicas de Reducción de la Varianza - métodos para reducir la varianza en la solución estimada para reducir el tiempo de cómputo en la simulación de Monte Carlo. Paralelización y Vectorización - algoritmos que permita que los métodos de Monte Carlo puedan ser implementados de manera eficiente en arquitecturas de computadoras avanzadas.
2.3.2 Aplicaciones La importancia actual del método Montecarlo se basa en la existencia de problemas que tienen difícil solución por métodos exclusivamente analíticos o numéricos, pero que dependen de factores aleatorios o se pueden asociar a un modelo probabilístico artificial (obtención de áreas de figuras no comunes, resolución de integrales de difícil solución, minimización de funciones, etc.). Gracias al avance en diseño de los ordenadores, los cálculos Montecarlo que en otro tiempo hubieran sido inconcebibles, hoy en día se presentan como asequibles para la resolución de ciertos problemas. En estos métodos el error ~
1/√N, donde N es el número de pruebas y, por tanto, ganar una cifra decimal en la precisión implica aumentar N en 100 veces. La base es la generación de números pseudo-aleatorios de los que nos serviremos para calcular probabilidades. En el caso que presentamos hemos hecho uso de la función rand() incluida en Matlab como generador de números pseudoaleatorios. Las pruebas realizadas, algunas de las cuales se llevaron a cabo en el tema 2.2, verifican su calidad a la hora de calcular números pseudoaleatorios sin tendencia aparente a la repetición ordenada. Ejemplo Aproximación de Pi Podemos aproximar Pi (=3.14159) mediante la generación de dos números pseudo-aleatorios para las componentes x y y de un vector simulado. Entonces podemos obtener, mediante el uso de teorema de Pitágoras, si estas coordenadas se encuentran dentro o fuera del círculo de radio unitario. Contamos estos aciertos, y después de hacer esto muchas veces (100 o más), podemos obtener un valor estimado de Pi. La precisión de la estimación depende de la cantidad de "vectores". Un código de ejemplo podría ser (suponiendo que se fija el radio = 1): Código Matlab: aciertos=0; n=100; for i=1:n x(i)=rand(); y(i)=rand(); dist=sqrt(x(i)^2+y(i)^2); if dist<=1 aciertos=aciertos+1; end end
valorpi=4*aciertos/n; disp(valorpi);
Aproximación de una Area (integral de una función f(x) De la misma manera podemos aproximar el valor de una area mediante la generación de dos números pseudo-aleatorios para las componentes x y y de las coordenadas de un punto. Entonces se puede obtener, mediante el uso de teorema del concepto de pertenencia, si estas coordenadas se encuentran dentro o fuera del área buscada. Contamos estos aciertos, y después de hacer esto muchas veces (100 o más), podemos obtener un valor estimado del área pedida. La precisión de la estimación depende de la cantidad de "puntos". Un código de ejemplo podría ser (suponiendo que se utiliza una función f(x) = x2): Código: aciertos=0; n=100; for i=1:n x(i)=rand(); y(i)=rand(); if y(i)<=x(i)^2 aciertos=aciertos+1; end end valorintx2=aciertos/n; disp(valorintx2);
Requerimientos del Método de Monte Carlo Los métodos de Monte Carlo necesitan diferentes Generadores de Números Aleatorios (RNG por las siglas en inglés). – Por ejemplo, cuando se simula la dirección de salida de una partícula y las interacciones de la partícula con el medio, las siguientes serán propiedades deseables: – Los atributos de cada partícula debe ser independientes unos de otros. – El atributo de todas las partículas debe extenderse a lo largo de todo el espacio de atributos. Es decir, cuando nos acercamos a un número infinito de partículas, las partículas lanzadas en un espacio deben cubrir el espacio por completo. – Las propiedades necesarias del generador de números aleatorios son: – Cualquier subsecuencia de números aleatorios no debe ser correlacionado con cualquier otra subsecuencia de números aleatorios. Por ejemplo, cuando se simulan las partículas lanzadas, no debemos generar patrones geométricos. – Un número de repetición aleatoria debe ocurrir sólo después de una gran generación de números aleatorios.
– Los números aleatorios generados deben ser uniformes. Este punto y el primero de ellos son poco relacionada. Para lograr una mayor uniformidad, se deben establecer algunas correlaciones entre los números aleatorios. – El RNG debe ser eficiente. Debe ser vectorizable con bajo costo operativo. Los procesadores en sistemas paralelos, no deben ser obligados a comunicarse entre sí. 2.3.3 Solución de Problemas Los resultados obtenidos en el punto anterior proveen la motivación para aplicar el Método de Monte Carlo, ya que indican que con un número suficientemente alto de experimentos, es posible estimar el parámetro deseado incurriendo en pequeño error con alta probabilidad, y permiten cuantificar asintóticamente la relación entre dos valores importantes en la aproximación (error y probabilidad) a través de la distribución normal. Sin embargo, es preciso ser cauteloso en la aplicación práctica del método, ya que las implementaciones reales no verifican las hipótesis de algunos teoremas. Por un lado, las limitaciones computacionales imponen un límite superior a los valores de n que se pueden emplear (y cuando se emplean números pseudoaleatorios, la naturaleza cíclica de estos hace que no sea posible obtener una cantidad arbitraria de muestras independientes). Ejemplo Supongamos que tenemos un satélite, que para su funcionamiento depende de que al menos 2 paneles solares de los 5 que tiene disponibles estén en funcionamiento, y queremos calcular φ la vida útil esperada del satélite (el tiempo promedio de funcionamiento hasta que falla, usualmente conocido en la literatura como MTTF - Mean Time To Failure). Supongamos que cada panel solar tiene una vida útil que es aleatoria, y está uniformemente distribuida en el rango [1000 hs, 5000 hs] (valor promedio: 3000 hs). Para estimar por Monte Carlo el valor de φ, haremos n experimentos, cada uno de los cuales consistirá en sortear el tiempo de falla de cada uno de los paneles solares del satélite, y observar cual es el momento en el cuál han fallado 4 de los mismos, esta es la variable aleatoria cuya esperanza es el tiempo promedio de f uncionamiento del satélite. El valor promedio de las n observaciones nos proporciona una estimación de φ. Experimento No. 1 2 3 4 5 6 Prom.
Ejercicio 1
Panel 1 3027 4162 3655 2573 2977 3756 -
Panel 2 1738 4029 2896 2649 2724 4190
Tiempo de Falla de Panel 3 Panel 4 Panel 5 2376 4685 4546 4615 3455 3372 1378 4010 4144 2117 3956 1281 1355 2268 3262 1749 3398 2581
Satélite, X 4546 4162 4010 2649 2977 3756 Sn/n=3683
Supongamos que se desea aproximar el valor de la base de los logaritmos
Neperianos (e=2.71828183…), se utiliza una función:
f ( x) e x Cuya integral (área bajo la curva) se define como: 1 1
0
f ( x)dx e x
e 1 e 0 1 0
1 e
Que de acuerdo al método de Montecarlo corresponde a: 1 1 N f ( x)dx 1 ac 0 e N Y por lo tanto, la aproximación de e:
e
1
1
N ac N
2.4 Problemas 1. Encontrar 20 números pseudo-aleatorios por el método no congruencial indicado, de acuerdo con los siguientes parámetros, elabore la tabla correspondiente: a) Por el método de cuadrado de medios, con semilla de 6727. b) Por el método de producto por factor constante, con semilla 7531 y factor constante de 5379. c) Por el método de producto de medios, con semillas de 4219 y 3753. d) Por el método de producto por factor constante, con semilla 5571 y factor constante de 3753. e) Por el método de cuadrado de medios, con semilla de 1017. f) Por el método de producto de medios, con semillas de 4491 y 8245 2. Encontrar 20 números pseudo-aleatorios a partir de la ecuación generadora indicada y los parámetros incluidos, elabore la tabla correspondiente: a) X i 1 (13X i 11) mod 32, con semilla x0=19
( X i 21) mod 64, con semilla x0=13. con semilla x0=27. X i 1 (19 X i ) mod 32, con semillas x0=21 y x1=43. X i 1 (21 X i 15 X i 1 ) mod 32, con semilla 25 X i 1 ( X i 17) mod 32, X i ( X i 1 X i n ) mod 64, i n 1, n 2, n 3,... con 4 semillas x 1=51, x2=23,
b) X i 1 c) d) e) f)
x3=27 y x4=43. 3. Realizar la prueba de la media y la varianza, además de la indicada en cada punto, con el respectivo nivel de aceptación, al conjunto de números que se da a continuación: 0.0645
0.8065
0.6452
0.0968
0.1613
0.5161
0.4839
0.5806
0.3871
0.8710
0.4516
0.6774
0.2258
0.3226
0.0323
0.9032
0.1935
0.4194
0.7742
0.7419
a) La prueba de corridas arriba y abajo con un nivel de aceptación del 90%. b) La prueba de chi-cuadrada con un nivel de aceptación del 95%.
c) La prueba de corridas encima y debajo de la media con una aceptación del 95%. d) La prueba de Kolmogorov-Smirnov con un nivel de aceptación del 90%. e) La prueba de Póker con una aceptación del 95%. 4. Realizar la aproximación de la integral de la función indicada, utilizando el conjunto de números que se da en las tablas siguientes, la tabla A es x y la tabla B es y. a) f(x)=x3 b) f(x)=x2 c) f(x)=xx d) f(x)=e-x
Tabla A
Tabla B
0.2525
0.7325
0.9659
0.9576
0.9677
0.6129
0.2581
0.9355
0.3756
0.6556
0.2962
0.6997
0.4839
0.1290
0.8065
0.4516
0.1075
0.9811
0.7734
0.958
0.0000
0.6774
0.3226
1.0000
0.1556
0.2557
0.8147
0.7764
0.5484
0.1935
0.8710
0.5161
0.4211
0.5382
0.3736
0.2796
0.0645
0.7419
0.3871
0.0323
