TEMA 4 NÚMEROS ALEATORIOS Y MÉTODOS DE MONTE CARLO

1. Métodos Numéricos y de Simulación TEMA 4NÚMEROS ALEATORIOS YMÉTODOS DE MONTE CARLO Rocío del Río Fernández Dpto. Electrónica y Electromagnetismo

2. Métodos Aleatorios Contenidos • Funciones de MATLAB relacionadas con la aleatoriedad:rand, randn, random y randperm • Generadores de números pseudo-aleatorios:Concepto, tests de caracterización y repetitividad • Métodos de Monte Carlo y aplicaciones: • Paseo aleatorio • Desintegración radioactiva • Integración numérica 1 sesión teórica + 1 sesión práctica

3. Saber utilizar las funciones rand y randn de MATLAB para la generación de secuencias aleatorias. Comprender la naturaleza pseudo-aleatoria de cualquier secuencia generada de manera determinista por un ordenador. Ser capaz de aplicar tests simples para determinar la uniformidad y aleatoriedad de secuencias pseudo-aleatorias. Comprender la utilidad de los métodos de Monte Carlo en la simulación de procesos físicos estocásticos. Ser capaz de evaluar integrales multi-dimensionales mediante el muestreo del valor medio. Métodos Aleatorios Objetivos

4. Métodos Aleatorios Bibliografía • R.H. Landau, M.J. Páez, C.C. Bordeianu: Computational Physics: Problem Solving with Computers, 2/E. Wiley, 2007.http://physics.orst.edu/~rubin/ Transparencias, videos, applets, … • S.R. Otto, J.P. Denier: An Introduction to Programming and Numerical Methods in MATLAB. Springer, 2005. • A. O’Hare: Numerical Methods for Physicists. 2005.http://www-teaching.physics.ox.ac.uk/computing/NumericalMethods/NMfP.pdf

5. Genera números reales aleatorios con distribución uniforme en el intervalo [0,1]. rand(M,N) genera una matriz MxN Métodos Aleatorios Funciones de MATLAB Función rand EJEMPLO_01.m>> … >> rand(1,100) >> …

6. Métodos Aleatorios Funciones de MATLAB >> rand(1,100) La función de distribución es tanto más visible cuanto mayor es el número de muestras generadas. >> rand(1,1e5)

7. Números reales aleatorios con distribución uniforme en el intervalo [0,1]. >> rand(1,1e5) Métodos Aleatorios Funciones de MATLAB Desplazamiento de Distribuciones Uniformes Números reales aleatorios con distribución uniforme en el intervalo [a,b]. >> a+(b-a)*rand(1,1e5) EJEMPLO_02.m a = -0.5 b = 3

8. Genera números reales aleatorios con distribución normal (gaussiana) de media 0 y desviación estándar 1. randn(M,N) genera una matriz MxN Métodos Aleatorios Funciones de MATLAB Función randn EJEMPLO_03.m>> … >> randn(1,100) >> …

9. Métodos Aleatorios Funciones de MATLAB >> randn(1,100) Media = -0.0177 Std = 1.0150 La función de distribución es tanto más visible cuanto mayor es el número de muestras generadas. >> randn(1,1e5) Media = 0.0021 Std = 0.9979

10. Números reales aleatorios con distribución normal de media 0 y desviación estándar 1. >> randn(1,1e5) Métodos Aleatorios Funciones de MATLAB Desplazamiento de Distribuciones Normales Números reales aleatorios con distribución normal de media  desviación estándar . >> mu+sigma*randn(1,1e5) EJEMPLO_04.m mu = 5 sigma = 0.2

11. Genera números reales aleatorios de acuerdo a distribuciones especificadas.random(NAME,A,B,C,M,N) Matriz MxN con distribución NAME y parámetros A, B, C Métodos Aleatorios Funciones de MATLAB Función random >> help random … NAME can be: 'beta' or 'Beta', 'bino' or 'Binomial', 'chi2' or 'Chisquare', 'exp' or 'Exponential', 'ev' or 'Extreme Value', 'f' or 'F', 'gam' or 'Gamma', 'gev' or 'Generalized Extreme Value', 'gp' or 'Generalized Pareto', 'geo' or 'Geometric', 'hyge' or 'Hypergeometric', 'logn' or 'Lognormal', 'nbin' or 'Negative Binomial', 'nbin' or 'Negative Binomial', 'ncf' or 'Noncentral F', 'nct' or 'Noncentral t', 'ncx2' or 'Noncentral Chi-square', 'norm' or 'Normal', 'poiss' or 'Poisson', 'rayl' or 'Rayleigh', 't' or 'T', 'unif' or 'Uniform', 'unid' or 'Discrete Uniform', 'wbl' or 'Weibull'.

12. Métodos Aleatorios Funciones de MATLAB Función randperm • Realiza una permutación aleatoria de números enteros. randperm(N) Permutación aleatoria de los enteros 1 a N. • Útil en la aleatorización de listas. EJEMPLO_05.m

13. Métodos Aleatorios Funciones de MATLAB Otros Ejemplos de Uso • Generación de N números naturales aleatorios con distribución uniforme en el intervalo [1,n].ceil(n*rand(N)+eps) EJEMPLO_07.m EJEMPLO_06.m

14. Métodos Aleatorios Generadores pseudo-aleatorios ¿Aleatoriedad Determinista? • Los ordenadores son deterministas  a las mismas entradas a un programa, las mismas salidas (salvo error). • ¿Cómo pueden entonces generar números aleatorios? • ¿Cómo de bien pueden hacerlo? • ¿Para qué nos sirve esto? • Simulación de eventos aleatorios: • movimiento térmico • juegos de azar • desintegración radioactiva • solución estadística de ecuaciones • …

15. Métodos Aleatorios Generadores pseudo-aleatorios Pseudo-Aleatoriedad • Los ordenadores son deterministas  a las mismas entradas a un programa, las mismas salidas (salvo error). • ¿Cómo pueden entonces generar números aleatorios? • Estrictamente hablando, no son capaces de generar secuencias aleatorias (sin correlación entre los números). • Si conocemos r1, r2, …, rm, siempre es posible determinar rm+1. • Realmente se generan secuencias pseudo-aleatorias. • ¿Cómo de bien pueden hacerlo? • Muy bien: • Baja correlación • Periodo de repetitividad alto • … siempre que el generador sea bueno.

16. Métodos Aleatorios Generadores pseudo-aleatorios Generadores Congruentes Lineales • Es el método más común para generar secuencias de números pseudo-aleatorios con distribución uniforme. • Genera números enteros en el intervalo [0,M-1] a partir del resto del cociente (rem = remainder) de números grandes: • r1 = semilla (seed), proporcionada por el usuario (o la máquina) • M = número muy grande (periodicidad máxima de la secuencia) • a = número grande; c = magia negra

17. Métodos Aleatorios Generadores pseudo-aleatorios Ejemplo ri = (ari-1 + c)modM, con a = 57, c = 1, M = 256, r1 = 10 • Secuencia pseudo-aleatoria en [0, M-1]: • r1 = 10 • r2 = (5710 + 1)mod256 = rem(571/256) = 59 • r3 = (5759 + 1)mod256 = rem(3364/256) = 36 • … • r257 = 10  Longitud de la secuencia (periodicidad M = 256) • Secuencia pseudo-aleatoria en [0, 1]  r/M • 0.0391 0.2305 0.1406 0.0195 0.1172 0.6836 0.9688 0.2227 …pero lógicamente con la misma periodicidad en la secuencia (M)

18. Métodos Aleatorios Generadores pseudo-aleatorios • Una prueba simple para chequear una secuencia pseudo-aleatoria es utilizar el córtex visual para reconocer los patrones de correlación. EJEMPLO_08.m (N = 2000) LCG con a = 57, c = 1, M = 256periodo 256 MATLAB randAlgoritmo Mersenne Twister, periodo (2^19937-1)/2

19. Métodos Aleatorios Generadores pseudo-aleatorios Tests de Aleatoriedad y Uniformidad • Mirar el listado de números generados. • Dibujar la gráfica de scattering r(i) frente i. • Evaluar el momento k-ésimo de la secuencia: • Si se cumple la ecuación, la distribución es uniforme. Si las desviaciones varían como 1/sqrt(N), la distribución es además aleatoria. • Evaluar la correlación con vecinos próximos: • Si se cumple la ecuación, los números no están correlacionados. Si las desviaciones varían como 1/sqrt(N), la distribución es además aleatoria.

20. Métodos Aleatorios Generadores pseudo-aleatorios Tests de Aleatoriedad y Uniformidad (cont.) • Dibujar la gráfica de scattering r(2i+1) frente r(2i). • Correr el cálculo o simulación con r1,r2,r3,… y con (1-r1),(1-r2),(1-r3),… Los resultados no deben diferir más allá de la estadística. • Correr el cálculo o simulación con una secuencia de números verdaderamente aleatorios tabulados y comparar con el generador de números pseudo-aleatorios. • Realizar el test 3 para k = 1, 3, 7 y N = 1e2, 1e4, 1e5. • Realizar el test 4 a la serie levemente correlacionada r1,(1-r1),r2,(1-r2),r3, (1-r3),… para N = 1e2, 1e4, 1e5.

21. Métodos Aleatorios Generadores pseudo-aleatorios Repetitividad en Secuencias Pseudo-Aleatorias • A veces resulta útil poder correr un cálculo o simulación con exactamente la misma secuencia pseudo-aleatoria: • Para el debugging de scripts • Para chequear si el método computacional es correcto • Para que una secuencia pseudo-aleatoria sea reproducible debemos utilizar siempre la misma semilla. • Para cambiar la secuencia, basta usar una semilla distinta: • Cambio automático del estado tras una ejecución de la función en MATLAB • Uso del reloj del sistema como semilla EJEMPLO_09.m

22. Métodos Aleatorios Métodos de Monte Carlo Monte Carlo Familia de métodos que involucran el uso de un número alto de muestras aleatorias para solucionar problemas matemáticos o físicos mediante técnicas estadísticas. • Simulación de procesos físicos aleatorios: • Movimiento térmico • Desintegración radioactiva • Método matemático: • Integración numérica por rechazo • Integración numérica por valor medio

23. Métodos Aleatorios Métodos de Monte Carlo Paseo Aleatorio • Procesos físicos en los que una partícula parece moverse aleatoriamente: • Movimiento browniano • Transporte de electrones en metales • Problema: ¿Cuántas colisiones (de media) debe sufrir una partícula para recorrer una distancia radial R? • Modelo: Caminar N pasos de longitud r en direcciones aleatorias.

24. Métodos Aleatorios Métodos de Monte Carlo Paseo Aleatorio - Teoría • ¿A qué distancia del origen llegamos después de N pasos? • Si las direcciones son aleatorias, para N alto, los términos cruzados se anulan: • Si damos un paseo de N pasos en direcciones aleatorias cubriendo una distancia total de Nrrms, acabamos el paseo (de media) a una distancia radial sqrt(N)rrms del punto de comienzo.

25. Métodos Aleatorios Métodos de Monte Carlo Paseo Aleatorio - Simulación EJEMPLO_10.m Media de 31 paseos de 1000 pasos 1 paseo de 1000 pasos

26. Métodos Aleatorios Métodos de Monte Carlo Desintegración Radioactiva • Proceso natural en el que una partícula (átomo, núcleo), sin ningún estímulo externo, se desintegra en otras partículas. • El instante en el que una partícula se desintegra es aleatorio. Es independiente de cuánto lleva existiendo la partícula o de qué le ocurre a otras a su alrededor. • La probabilidad P(t) de desintegración por unidad de tiempo y por partícula es constante: • P(t) = probab desinteg / t / partícula = -l • Al ir diminuyendo el número de partículas N, lógicamente también lo hará el número de desintegraciones: • N(t), DN/Dt ↓ con t

27. Métodos Aleatorios Métodos de Monte Carlo Desintegración Radioactiva - Simulación EJEMPLO_11.m La simulación del proceso demuestra: • Ley de desintegración exponencial: • N(t)  e-lt si N alto • Proceso estocástico si N bajo

28. Box Métodos Aleatorios Métodos de Monte Carlo Integración Numérica Tirar piedras en un estanque con forma irregular como método para determinar su área: Señalice una caja que encierre completamente el estanque y quite todas las piedras que contenga. Mida el área de esa caja (Abox). Coja un buen montón de piedras y láncelas al aire en direcciones aleatorias. Cuente el número de piedras que caen en el estanque (Npond) y el número de piedras que caen al suelo dentro de la caja (Nbox). Asumiendo que tiró las piedras uniforme y aleatoriamente, Npond debe ser proporcional a Apond.

29. Pond r = 1 (0,0) Box Métodos Aleatorios Métodos de Monte Carlo Integración Numérica por Rechazo EJEMPLO_12.m Determinar el número p utilizando la técnica de muestreo anterior. Imagine un estanque circular encerrado en un cuadrado de lado 2. Sabemos que Generamos una secuencia de números aleatorios {ri} en [-1,1]. Asignamos (xi, yi) = (r2i-1, r2i) para i = 1, …,N. Si xi2 + yi2 < 1, Npond = Npond+1. En otro caso, Nbox = Nbox+1. Usar Aumente N hasta obtener p con 3 cifras decimales. • ¿N=1e5, 1e6? • ¿Tiempo de cómputo? • Pruebe EJEMPLO_12b.m

30. Métodos Aleatorios Métodos de Monte Carlo Integración Numérica por Valor Medio • La técnica estándar de Monte Carlo para integración está basada en el Teorema de Valor Medio: • Se usa una secuencia xi de N números aleatorios uniformes en [a,b] para muestrear <f>: • El error en el valor de la integral disminuye como 1/sqrt(N).

31. Métodos Aleatorios Métodos de Monte Carlo Integrales Multi-Dimensionales • La integración por valor medio puede generalizarse fácilmente a espacios multi-dimensionales: • Existen otros métodos de integración mejores para integrales simples y dobles, pero cuando la dimensionalidad es alta las técnicas de Monte Carlo son las mejores. EJEMPLO_13.mIntegración en 3D EJEMPLO_14.mIntegración en 10D