Transcript Métodos Aleatorios - IHMC Public Cmaps

Métodos Numéricos y de Simulación
TEMA 4
NÚMEROS ALEATORIOS Y
MÉTODOS DE MONTE CARLO
Rocío del Río Fernández
Dpto. Electrónica y Electromagnetismo
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Objetivos
 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.
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Funciones de MATLAB
Función rand
 Genera números reales aleatorios con distribución uniforme en el
intervalo [0,1].
rand(M,N)  genera una matriz MxN
EJEMPLO_01.m
>> …
>> rand(1,100)
>> …
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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)
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Funciones de MATLAB
Desplazamiento de Distribuciones Uniformes
Números reales aleatorios con
distribución uniforme en el
intervalo [0,1].
Números reales aleatorios con
distribución uniforme en el
intervalo [a,b].
>> rand(1,1e5)
>> a+(b-a)*rand(1,1e5)
EJEMPLO_02.m
a = -0.5
b = 3
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Funciones de MATLAB
Función randn
 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
EJEMPLO_03.m
>> …
>> randn(1,100)
>> …
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Funciones de MATLAB
Desplazamiento de Distribuciones Normales
Números reales aleatorios con
distribución normal de media 0
y desviación estándar 1.
Números reales aleatorios con
distribución normal de media 
desviación estándar .
>> randn(1,1e5)
>> mu+sigma*randn(1,1e5)
EJEMPLO_04.m
mu = 5
sigma = 0.2
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Funciones de MATLAB
Función random
 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
>> help random
…
NAME can be:
'beta' or
'bino' or
'chi2' or
'exp'
or
'ev'
or
'f'
or
'gam'
or
'gev'
or
'gp'
or
'geo'
or
'hyge' or
'logn' or
'nbin' or
'Beta',
'Binomial',
'Chisquare',
'Exponential',
'Extreme Value',
'F',
'Gamma',
'Generalized Extreme Value',
'Generalized Pareto',
'Geometric',
'Hypergeometric',
'Lognormal',
'Negative Binomial',
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
'nbin'
'ncf'
'nct'
'ncx2'
'norm'
'poiss'
'rayl'
't'
'unif'
'unid'
'wbl'
or
or
or
or
or
or
or
or
or
or
or
'Negative Binomial',
'Noncentral F',
'Noncentral t',
'Noncentral Chi-square',
'Normal',
'Poisson',
'Rayleigh',
'T',
'Uniform',
'Discrete Uniform',
'Weibull'.
Métodos Numéricos y de Simulación – Curso 2010/11
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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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_06.m
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
EJEMPLO_07.m
Métodos Numéricos y de Simulación – Curso 2010/11
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
• …
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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.
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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.
 ari1  c 
ri  (ari1  c)modM  rem

 M 
 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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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)
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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 = 256
MATLAB rand
periodo 256
Algoritmo Mersenne Twister, periodo (2^19937-1)/2
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Generadores pseudo-aleatorios
Tests de Aleatoriedad y Uniformidad
1. Mirar el listado de números generados.
2. Dibujar la gráfica de scattering r(i) frente i.
3. Evaluar el momento k-ésimo de la secuencia:


1
1 N k
1
x   x i   xkP(x)dx   1 / N 
0
N i1
k 1
k
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.
4. Evaluar la correlación con vecinos próximos:


1
1
1 N
1
C(k)   x ix ik   dx xyP(x , y)dy   1 / N 
0
0
N i1
4
k  1,2,...
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.
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Generadores pseudo-aleatorios
Tests de Aleatoriedad y Uniformidad (cont.)
5. Dibujar la gráfica de scattering r(2i+1) frente r(2i).
6. 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.
7. 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.
8. Realizar el test 3 para k = 1, 3, 7 y N = 1e2, 1e4, 1e5.
9. Realizar el test 4 a la serie levemente correlacionada r1,(1-r1),r2,(1-r2),r3,
(1-r3),… para N = 1e2, 1e4, 1e5.
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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.
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
Métodos Aleatorios
Métodos de Monte Carlo
Paseo Aleatorio - Teoría

¿A qué distancia del origen llegamos después de N pasos?
R2  (x1  x2  ...  xN )2  (y1  y2  ...  yN )2
 x12  ...  xN2  2x1x2  2x1x3  2x2x1  ...  (x  y)
•
Si las direcciones son aleatorias, para N alto, los términos
cruzados se anulan:
R2  x 12  ...  x N2  y 12  ...  y N2  N r 2
R  Nrrms 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.
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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), N/t ↓ con t
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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:
1.
2.
3.
4.
5.
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.
Npond
Npond  Nbox

Apond
Abox
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)

Apond 
Npond
Npond  Nbox
Box
Abox
Métodos Numéricos y de Simulación – Curso 2010/11
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.
1.
2.
3.
4.
Imagine un estanque circular encerrado en un cuadrado
de lado 2.
Sabemos que  dA  pr 2  p
Generamos una secuencia de números aleatorios {ri} en
[-1,1].
Pond
r=1
(0,0)
Asignamos (xi, yi) = (r2i-1, r2i) para i = 1, …,N.
xi2
yi2
5.
6.
Si +
Usar
7.
Npond  Nbox
Aumente N hasta obtener p con 3 cifras decimales.
Box
< 1, Npond = Npond+1. En otro caso, Nbox = Nbox+1.
Apond 
Npond
Abox  p
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
¿N=1e5, 1e6?
¿Tiempo de cómputo?
Pruebe EJEMPLO_12b.m
Métodos Numéricos y de Simulación – Curso 2010/11
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:
b
I   f(x)dx  (b  a) f
a

Se usa una secuencia xi de N números
aleatorios uniformes en [a,b] para
muestrear <f>:
1 N
f   f(x i )
N i1

El error en el valor de la integral disminuye como 1/sqrt(N).
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
Métodos Numéricos y de Simulación – Curso 2010/11
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:
b
d
f
a
c
e
I   dx dy f(x , y, z)dz  (b  a)(d  c)(f  e) f

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.m
Integración en 3D
1
1
1
0
0
0
EJEMPLO_14.m
Integración en 10D
I   dx1  dx2  dx3 (x1  x2  x 3 )
Rocío del Río – Dpto. Electrónica y Electromagnetismo (US)
2
1
1
0
0
I   dx1 ... dx10 (x1  ...  x10 )2
Métodos Numéricos y de Simulación – Curso 2010/11