Monte Carlo:
experimentos simulados
Tema 2
Itziar Aretxaga
Números aleatorios: distribución uniforme
0 ≤ x ≤ 1 P(x)dx=dx
Es el generador básico de números aleatorios.
Todos los lenguajes de programación cuentan con
uno de estos generadores:
.
Ejem. FORTRAN:
iseed=1.
.
x=ran(iseed)
P(x)
1
0
0
1
x
normalmente están basados en generadores lineales congruenciales: .
.
Ij+1=aIj+b (mod m)
que sufren de:
1. Dados k números aleatorios, estos se distribuyen en un plano k-1.
2. Los bits de menor orden están más asociados que los de mayor orden
Es por esto que se recomienda utilizar rutinas portátiles especializadas. En
Numerical Recipes, por ejemplo:
. ran0: tiene un periodo de repetición 2x109 y sufre los fallos de los
generadores lineales congruenciales estándar.
. ran1: es como ran0 pero sin los problemas de dimensionalidad y
correlación. Es uno de los generadores recomendados.
. ran2: tiene un periodo de 2x1018, pero es más lento.
...
(Press et al., “Numerical Recipes”)
Números aleatorios
(Cortesía de Leonardo Sandoval)
Colección DIEHARD de tests de aleatoriedad (George Marsaglia,
ftp://stat.fsu.edu/diehard/index.html)
Otros tests de aleatoriedad y generadores alternativos de números
aleatorios (http://burtleburtle.net/bob/rand/testsfor.html)
Números aleatorios: método de transformación
Supongamos que queremos generar números aleatorios que sigan una
densidad de probabilidad P(y), que tiene asociada una probabilidad
acumulada F(y)
(Fig. © “Numerical Recipes”)
Esta distribución se puede asociar a la distribución uniforme P(x)
P ( y ) dy  P ( x ) dx  x 

y
P ( y ) d y   F ( y )
a
si F(y) es invertible, entonces el número aleatorio y=F−1(x). Por lo tanto,
se generan números aleatorios x bajo una distribución uniforme, y se
transforman en números aleatorios y bajo la distribución P(y).
Ejemplo: distribución exponencial P(y)=e−y
e−ydy=dx
x=1−e−y
0≤ x≤1
(Press et al., “Numerical Recipes”)
y=−ln(1−x)
x
1
0≤y≤∞
0
y
Números aleatorios: gaussianas
Sean x,y dos números aleatorios generados por distribuciones normales.
Si son independientes, su distribución sobre un plano será
 (x2  y2)
P ( x, y ) 
exp  

2
2


1
P ( d , ) 
o en coordenadas polares (R,θ) y haciendo d=R2
 ( x, y )
 ( d , )
P ( x, y ) 
1 1
2 2
exp(  d / 2 )
lo que es equivalente al producto de una distribución exponencial de vida
media 2, y una distribución uniforme definida en el intervalo [0,2π].
Ésta es la base de la transformación de Box-Müller:
Sean dos números aleatorios u1, u2 derivados de una distribución uniforme.
Se realizan las transformaciones
R   2 ln u 1
x  R cos  
 2 ln u 1
cos(2  u 2 )
  2 u 2
y  R sin  
 2 ln u 1
sin(2  u 2 )
2
que nos llevan a dos números aleatorios x,y cuya probabilidad sigue una
distribución gaussiana.
Puesto que las transformaciones dependen de funciones trigonométricas, no
son muy eficientes para el cálculo computacional.
(Press et al., “Numerical Recipes”)
Números aleatorios: gaussianas
(−1,1)
(1,1)
R
θ v2
v1
Para hacer el algoritmo de Box-Müller más rápido
se definen las variables
v1 =2u1−1
v2 =2u2−1
Se generan números hasta que (v1,v2) se
encuentre dentro del círculo de radio R=1.
cos  
(−1,−1)
v1
R
(1,−1)
sin  
v2
R
  2 ln d 
x  v 1

d




v1
2
2
1/ 2
2
1/ 2
(v 1  v 2 )
v2
2
(v 1  v 2 )
1/ 2
  2 ln d 
y  v2

d


1/ 2
para d ≤ 1
.
Estas transformaciones modificadas
son más eficientes en el cálculo.
(Press et al., “Numerical Recipes”)
Números aleatorios: método del rechazo
Supongamos que queremos generar números aleatorios que sigan una
densidad de probabilidad P(y), cuya probabilidad acumulada F(y), o bien
no sea analítica o no sea invertible.
(Fig. © “Numerical Recipes”)
Se busca una función envolvente f(y) que tenga una integral I(y)=∫ f(y)dy
finita e invertible. Se genera un número aleatorio de distribución uniforme
x en el intervalo (0,A), entonces y=F−1(x) es un número aleatorio de la
distribución envolvente f(y). Si ahora se genera un segundo número
aleatorio de distribución uniforme x2 en el intervalo (0,f(y)), entonces y es
aleatorio de la distribución P(y) si x2 ≤ P(y).
Ejemplo: distribución poissoniana
y se expande el área como si fuera
una distribución continua
f ( y) 
 e
j
P ( j) 
j!
c0
1  ( y  y0 ) / a0
y  a 0 tan(  x )  y 0
(Press et al., “Numerical Recipes”)

2
2
(Fig. © “Numerical Recipes”)
Monte Carlo: método de la fuerza bruta
Simulación de los resultados de un experimento utilizando una
computadora y un generador de números aleatorios. Este tipo de análisis
se utiliza cuando un cálculo es difícil de realizar por otros métodos
numéricos o algebraicos, o cuando somos demasiado vagos como para
pensar en cómo solucionarlo por métodos más y
elegantes.
Ejemplo clásico: área debajo de una curva.
Dada un área A fácil de medir, que contiene una
curva f(x) difícil de integrar, se puede calcular el
área debajo de la curva mediante la generación
N veces de dos números aleatorios (x,y) que
representen las coordenadas. Se cuentan los puntos
por encima y por debajo de la curva.
f(x)
x
o

f ( x ) dx  A
n de puntos bajo la curva
o
n de puntos
Este argumento se puede aplicar también a volúmenes, por ejemplo para calcular el
volumen comprendido por los censos de pincel.
El error del cálculo es proporcional a 1 /
N
Monte Carlo: método de la fuerza bruta
Ejemplo: cálculo del ángulo sólido observable por temporada, descontando el
espacio ocupado por el plano galáctico.
(O. Vega, 2001, MSc, INAOE)
Monte Carlo: método de la fuerza bruta
Ejemplo: curvas de luz de cúmulos
estelares de edades 10  60 Myr
(Aretxaga & Terlevich 1994, MNRAS,
269, 462):
Cada flecha señala el tiempo en el
que explota una SN, derivado de
una distribución uniforme en un
intervalo grande de tiempo (103
años). Este método se aproxima a
una distribución poissoniana de
media igual a la tasa de
explosiones señalada en cada
recuadro.
Por medio de un cálculo así, se
pueden medir fácilmente
propiedades tales como la
desviación cuadrática media de las
curvas de luz, etc.
Monte Carlo: método de la fuerza bruta
Ejemplo: curvas de luz de cúmulos
estelares de edades 10  60 Myr
(Aretxaga & Terlevich 1994, MNRAS,
269, 462):
Cada flecha señala el tiempo en el
que explota una SN, derivado de
una distribución uniforme en un
intervalo grande de tiempo (103
años). Este método se aproxima a
una distribución poissoniana de
media igual a la tasa de
explosiones señalada en cada
recuadro.
Por medio de un cálculo así, se
pueden medir fácilmente
propiedades tales como la
desviación cuadrática media de las
curvas de luz, etc.
Monte Carlo: método de la fuerza bruta
Ejemplo: curva de luz real de NGC4151
(en el recuadro de arriba) y simulación
utilizando la superposición de modelos
simples de SN (en el recuadro de
abajo).
El Monte Carlo siempre resulta más
convincente si se incluye además una
simulación de las condiciones de
observación, tales como errores
fotométricos (que se pueden aproximar
por distribuciones gaussianas), o el
muestreo temporal de la curva de luz.
(Aretxaga & Terlevich 1994)
Métodos de remuestreo:
Monte Carlos rápidos y sucios
Se suelen utilizar para calcular errores, estimar sesgos, etc.
♦ Filosofía: construir poblaciones hipotéticas derivadas del mismo conjunto
de observaciones, cada una de las cuales se analiza para determinar cómo
los indicadores estadísticos cambian con fluctuaciones aleatorias de los
datos.
♦ Ventajas: no necesita ninguna suposición sobre la distribución
subyaciente, y si existen efectos de selección en la muestra de datos,
estos también se consideran.
Remuestreos habituales (Feigelson & Babu, “Astrostatistics”):
• mitad de la muestra (Mahalonovis 1946): Selecciona aleatoriamente la
mitad de la muestra, n veces.
• “jacknife” (navaja): selecciona las n muestras que contienen n−1 puntos
diferentes de la muestra original.
• “bootstrap” (agujetas): genera un gran número de muestras con n
puntos seleccionados al azar , cada una de las cuales puede contener
puntos repetidos u omitidos de la muestra original.
.
“To pull oneself up by one´s bootstraps”
Métodos de remuestreo: bootstrap
Ejemplo: cálculo del error en la
medida de la mediana del índice
xi  x que define la intensidad
de la variabilidad de una muestra
de QSOs (Hook et al 1994, MNRAS,
268, 305)
Métodos de remuestreo: bootstrap
Ejemplo: cálculo del error en los
parámetros derivados del ajuste de
un modelo por la minimización de
χ2, y elipsoides de confianza
asociados al cálculo.
(Numerical Recipes, Press et al.)
Descargar

Simulaciones