UNIDAD III: SIMULACION MONTE CARLO
El método Montecarlo es un método numérico que permite resolver problemas físicos y matemáticos mediante la simulación de variables aleatorias. Lo vamos a considerar aquí desde un punto de vista didáctico para resolver un problema del que conocemos tanto su solución analítica como numérica. El método Montecarlo fue bautizado así por su clara analogía con los juegos de ruleta de los casinos, el más célebre de los cuales es el de Montecarlo, casino cuya construcción fue propuesta en 1856 por el príncipe Carlos III de Mónaco, siendo inaugurado en 1861.
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ística artificial (resolución de integrales de muchas variables, minimización de funciones, etc.). Gracias al avance en diseño de los ordenadores, cálculosMontecarlo 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 aleatorios de los que nos serviremos para calcular probabilidades. Conseguir un buen generador de estos números así como un conjunto estadístico adecuado sobre el que trabajar son las primeras dificultades con la nos vamos a encontrar a la hora de utilizar este método. En el caso que presentamos hemos hecho uso de la funciónrandom() incluida en la clase Math que la máquina virtual Java trae por defecto como generador. Las pruebas realizadas, algunas de las cuales se propondrán como ejercicio, verifican su calidad a la hora de calcular números aleatorios sin tendencia aparente a la repetición ordenada.
Para resolver la ecuación elíptica de nuestro problema usando el método de Montecarlo, se ha dividido el recinto bidimensional en una malla cuadrada de puntos. Todos los situados en su frontera se consideran inicializados a un valor de temperatura conocido. Suponemos en principio una partícula situada en uno de los puntos y que tiene la posibilidad de moverse libremente por todos los que constituyen la malla. La única condición que imponemos es que en un solo salto, su movimiento se limite a los 4 nodos vecinos, los situados su izquierda, derecha, arriba o abajo. La probabilidad de elegir una cualquiera de las 4 direcciones posibles es la misma. Dejando a la partícula viajar por toda la red sin más restricciones contamos el número de veces que, partiendo de un mismo punto de coordenadas (i,j) sale por cada uno de los que constituyen la frontera, momento en el cual suponemos que ha terminado su viaje. Considerando un número elevado de pruebas podemos calcular la probabilidad de que, partiendo de un mismo punto, salga por cada uno de los puntos del contorno después de recorrer una trayectoria aleatoria. Los detalles de camino seguido desde el inicio hasta el final del viaje no nos importan, tan solo nos vamos a fijar en el número de veces que sale del recinto por cada uno de los puntos posibles.
La temperatura a la que se encuentra el punto desde donde ha partido la partícula es la suma, extendida a todos los puntos frontera (if,jf), de la temperatura de dichos puntos (determinada por las condiciones de contorno) y por la probabilidad de que estando en (i,j) salga por (if,jf). Ver ec.(6).
Si tomamos una malla pequeña de 10x10 (salvo consideraciones de simetría) hay que calcular probabilidades para 102puntos. Una precisión razonable requerimos que para cada uno de ellos hay que calcular ~106 trayectorias aleatorias. Con sólo estas estimaciones podemos aventurar que el tiempo de computación requerido para solucionar la ecuación de Laplace en una malla pequeña va a ser superior al necesario en cualquiera de los otros métodos propuestos.
Bajo el nombre genérico de "simulaciones de Monte Carlo" suelen englobarse todas las simulaciones que emplean números aleatorios para resolver problemas estáticos (estocásticos o deterministas), es decir, problemas en los cuales el tiempo no juega un papel relevante.
Supongamos por ejemplo que queremos evaluar la integral:
donde g(x) es una función real monovaluada. La resolucion de este problema determinista mediante simulacion de Monte Carlo es la siguiente:
Sea Y una variable aleatoria definida de la forma:
donde X es una variable aleatoria uniforme en el intervalo [a,b], U(a,b).
El valor esperado de Y, E(Y), es:
donde la función:
es la función densidad de probabilidad de la variable aleatoria uniforme en el intervalo [a,b], U(a,b).
Asi pues. el problema de evaluar la integral se ha reducido la estima del valor esperado E(Y). En particular, E(Y) puede estimarse calculando la media de cierto número, n, de muestras:
donde X1,...,Xn son variables aleatorias U(a,b) IID.
Puede demostrarse que E(E(Y)(n))=I, es decir, que E(Y)(n) es un estimador no sesgado de I, y Var(E(Y)(n))=(Var(Y))/n. Como Var(Y) es un número fijo, E(Y)(n) puede hacerse arbitrariamente próximo a I para un n suficientemente grande.
Ejemplo: Supongamos que se evaluara la integral:
La integral posse primitiva analitica (-cos(x)), y el valor de la integral definida entre 0 y
p es 2.
En la práctica, si el integrando tiene "buen comportamiento" (ni la funcion ni sus derivadas presentan discontinuidades) suelen ser preferibles otros métodos numéricos de integración al de Monte Carlo.
El siguiente programa realiza la simulación:
|
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define PI 3.141593
struct intervalo{
float inicio
}
float Uniforme_interv_A_B(interv)
struct intervalo interv;
{ |
/* Generacion de U(A,B) */
float v_unif;
v_unif=rand()/32727.0;
if (v_unif>1.0) v_unif=1.0;
return(interv.inicio+(interv.fin-interv.inicio)*v_unif);
}
float Func_integrando(x)
float x;
{
return(sin(x));
}
void Entrada_datos(num_muestras,interv_integ)
unsigned int *num_muestras;
struct intervalo *interv_integ;
{
clrscr();
/* Numero de muestras */
printf("tNumero de muestras (1-65535): ");
scanf("%d",num_muestras);
/* Intervalo de integracion */
interv_integ->inicio=0;
interv_integ->fin=PI;
/* Inicializacion rutina generacion numeros aleatorios */
randomize();
return;
}
void Informe(num_muestras,estima_integ)
unsigned int num_muestras;
float estima_integ;
{
printf("tI N F O R M Ent-------------nn");
printf("Numero de muestras: %dn",num_muestras);
printf("Estima de la integral: %fn",estima_integ);
}
main()
{
unsigned int num_muestras,
contad=0;
float func_acumulada=0,
estima_integ;
struct intervalo interv_integ;
Entrada_datos(&num_muestras,&interv_integ);
while (contad<num_muestras) {
func_acumulada+=Func_integrando(Uniforme_interv_A_B(interv_integ));
contad++;
}
estima_integ=(interv_integ.fin-interv_integ.inicio)*
func_acumulada/num_muestras;
Informe(num_muestras,estima_integ);
return;
}
Simulando varias veces, para diferentes valores del número de muestras, obtenemos los siguientes resultados.
Simulacion 1:
| n |
10 |
50 |
100 |
500 |
1000 |
| E(Y)(n) |
1.980 |
1.832 |
2.021 |
1.958 |
2.001 |
Simulacion 2:
| n |
10 |
50 |
100 |
500 |
1000 |
| E(Y)(n) |
2.302 |
1.845 |
1.799 |
1.971 |
1.993 |
Simulacion 3:
| n |
10 |
50 |
100 |
500 |
1000 |
| E(Y)(n) |
2.286 |
1.926 |
1.801 |
1.999 |
1.998 |
Simulacion 4:
| n |
10 |
50 |
100 |
500 |
1000 |
| E(Y)(n) |
2.371 |
1.867 |
1.926 |
1.943 |
1.937 |
Simulacion 5:
| n |
10 |
50 |
100 |
500 |
100 |
| E(Y)(n) |
2.300 |
1.968 |
1.916 |
1.975 |
2.004 |
Deben establecerse criterios para determinar el tamaño de la muestra n, que permita obtener una determinada precisión.
Similarmente, si queremos evaluar la integral:
podemos realizar el cambio de variable:
con lo cual, reescribimos la integral:
