Curso: M´ etodos de Monte Carlo. Unidad 2, Sesi´ on 3: Estimaci´ on de vol´ umenes Departamento de Investigaci´ on Operativa Instituto de Computaci´ on, Facultad de Ingenier´ıa Universidad de la Rep´ ublica, Montevideo, Uruguay dictado semestre 1 - 2010

Contenido: 1. Introducci´ on y definiciones. 2. Estimaci´ on de vol´ umenes. 3. Ejemplo. 4. Ejercicio.

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

1

Monte Carlo para estimar vol´ umenes e integrales Esta unidad introduce los problemas b´asicos que surgen al aplicar el m´etodo de Monte Carlo para resolver el conjunto de problemas encontrados con mayor frecuencia en el ´area del c´alculo num´erico. En su forma m´as simple, el problema consiste en evaluar el volumen de una regi´on acotada en un espacio euclidiano multi-dimensional. Una formulaci´ on m´as general consiste en la evaluacion de la integral de una funci´ on en una regi´ on de este estilo. El m´etodo de Monte Carlo con frecuencia es una manera competitiva de resolver este problema (m´as a´ un, muchas veces es la u ´nica forma pr´actica para lograrlo). Las ventajas del m´etodo surgen cuando la forma de la regi´on de inter´es hace imposible resolver el problema por m´etodos anal´ıticos, y, en el caso general de la integracion de funciones, cuando las propiedades de suavidad Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

2

y de variaci´ on del integrando son desconocidas o impiden la aplicaci´on de t´ecnicas de evaluaci´ on num´erica cl´asicas.

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

3

Definiciones Esta unidad se centra en los problemas relacionados a la obtenci´on de estimadores puntuales y de estimadores de intervalo bajo forma de intervalos de confianza. • Estimador puntual: denota una expresi´ on matem´atica que transforma un conjunto de datos (derivados de una muestra obtenida experimentalmente) en un u ´nico n´ umero, conocido como estimaci´ on puntual, que sirve como aproximaci´ on para el valor exacto (pero desconocido) del volumen (o m´as generalmente de la integral) de inter´es. • Estimador de intervalo o Intervalo de confianza: denota una expresi´ on matem´atica que transforma un conjunto de datos (derivados de una muestra obtenida experimentalmente) para calcular un intervalo num´erico que, con una probabilidad especificada, contiene el valor exacto del volumen. Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

4

• En la pr´actica, en general se emplea un intervalo de confianza para complementar y servir como indicaci´ on del error potencial cometido por un estimador puntual. • Un punto central en todo m´etodo de Monte Carlo es el n´ umero de sorteos independientes o replicaciones que deben realizarse para garantizar una cierta cota del error. Este n´ umero es llamado tama˜ no de muestra, y depende del error aceptable en cada estudio, de las caracter´ısticas del problema (y del tiempo computacional disponible). • Muchas veces es posible (eventualmente empleando algo de informaci´on previamente conocida) calcular antes del muestreo, un tama˜ no de muestra de peor caso que garantiza que, si se emplea al menos este n´ umero de experimentos, el estimador puntual resultante tendr´a un error de acuerdo a las especificaciones deseadas. En general el tama˜ no de muestra de peor caso resulta en un esfuerzo computacional mayor que el estrictamente necesario, pero conocierlo con anticipaci´ on puede ser u ´til en varios contextos. Por un lado, si no Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

5

hay limitaciones en tiempo de c´alculo, el usuario puede emplear este tama˜ no de muestra para asegurarse la calidad del resultado. Si por el contrario el tiempo de c´alculo es una limitaci´ on importante, es posible emplear el conocimiento previo del tama˜ no de muestra de peor caso para varias alternativas como forma de elegir aquella m´as conveniente, y tambi´en para estudiar como mejorar la eficiencia computacional del m´etodo. • En algunos casos tambi´en es posible calcular a priori tama˜ no de muestra de mejor caso, que dan una cota inferior por debajo de la cu´al es seguro que no puede alcanzarse las especificaciones de error requeridas. La relaci´ on entre el tama˜ no de muestra de mejor y peor caso sirven como indicadores de que tan pesimista es esta segunda medida.

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

6

Estimaci´ on de vol´ umenes • Sea R una regi´ on de volumen desconocido λ(R), en el hipercubo unitario m-dimensional J m = [0, 1]m = [0, 1] × · · · × [0, 1]. • Si R fuera de tama˜ no arbitrario (pero acotada), supondremos que es posible transformarla para que quede incluida en J m. • Supongamos que R est´a definida por un conjunto de desigualdades y de relaciones impl´ıcitas entre las variables espaciales 0 ≤ xi ≤ 1, i = 1, . . . , m, y que la forma de R es tal que el c´alculo exacto de λ(R) es intratable. • Si disponemos de un procedimiento que genera una secuencia de puntos (j) (j) Xm,n = {x(j) = (x1 , . . . , xm ∈ J m; j = 1, . . . , n}, es posible ¯ aproximar λ(R) a trav´es de λ(R) calculada por el algoritmo que sigue: Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

7

Procedimiento Volumen (region R, entero n) Entrada: region R, n tama˜ no de la muestra ¯ Salida: λ(R) 1. S = 0; /* Inicializaci´ on */ 2. For j = 1, . . . , n do 2.1 Generar x(j) de Xm,n; 2.2 If x(j) ∈ R then φ(x(j)) = 1 else φ(x(j)) = 0; 2.3 S = S + φ(x(j)); /* Acumular*/ ¯ 3. λ(R) = S/n;

¯ • La exactitud de λ(R) depender´a directamente de las propiedades de R y de Xm,n. Por ejemplo, una forma sencilla es generar n puntos de evaluaci´ on en la malla m − dimensional: Xm,n = {x = (x1 , . . . , xm ) : xi = (zi + 1/2)/k; zi = 0, 1, . . . , k − 1; i = 1, . . . , m, n = k

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

m

)

8

1


















































































































































0



































































x_1







0




x_2




x_2




1


1




• En este caso, cada uno de los n puntos en Xm,n es el centro de un hipercubo de m dimensiones y volumen 1/k m = 1/n, y de esta forma S/n, el volumen total de cada uno de los S m-cubos con centros en R, ser´a una aproximaci´ on de λ(R).

x_1

1

Figure 1: Ejemplo de una superficie a estimar, y la estimaci´on a partir de la malla 2-dimensional. Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

9

• La figura 1 presenta un ejemplo en 2 dimensiones (en el que el volumen es en realidad una superficie), con la malla correspondiente y la estimaci´ on que surge de la misma (figura a la derecha). Es claro que se cometen errores de dos tipos, por un lado 2-cubos (cuadrados) que estan parcialmente fuera de la figura, por otro partes de la figura que no ¯ son cubiertas por ning´ un 2-cubo. El error total, λ(R) − λ(R), es la diferencia de estos dos errores. • En el caso 2-dimensional, sea s(R) el largo de la frontera de R. Entonces el error absoluto tiene la siguiente cota superior: ¯ |λ(R) − λ(R)| ≤ s(R)/k. Este error surge de considerar que los errores pueden agruparse en un rect´angulo que tiene, en el peor caso, largo s(R), y ancho 1/k (el ancho de cada 2-cubo de la malla), y de usar el ´area de este rect´angulo como cota del error en el peor caso. • En general, si s(R) es la “superficie” (en n − 1 dimensiones) de R, la Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

10

misma f´ ormula sigue siendo v´alida: ¯ |λ(R) − λ(R)| ≤ s(R)/k. Dado que n = k m, podemos expresar esta misma f´ormula as´ı: ¯ |λ(R) − λ(R)| ≤ s(R)/n1/m, lo que demuestra que la tasa de convergencia del error disminuye al aumentar la dimensi´ on m del problema. En particular, para garantizar un error absoluto de a lo sumo  ∈ (0, 1), es necesario tomar una malla (de dimensi´ on m) con al menos n(m, ) = d[s(R)/]me puntos (notaci´ on: dxe = ceiling(x), menor entero mayor o igual a x). • Por lo tanto, si tuvi´eramos una secuencia de regiones de dimensi´on creciente e igual volumen, para obtener igual precisi´on  necesitar´ıamos un n´ umero de puntos que crece exponencialmente con la dimensi´on m. Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

11

Estimaci´ on de volumen por Monte Carlo • La forma m´as sencilla de aplicar Monte Carlo para resolver este problema consiste en usar el mismo algoritmo anterior, con la u ´nica modificaci´ on consistente en sustituir el empleo de la secuencia determin´ıstica Xm,n = (x(1), . . . , x(n)) por una muestra aleatoria (j) (j) independiente (X(1), . . . , X(n)) donde cada X(j) = (X1 , . . . , Xm ) se sortea de acuerdo a la distribuci´ on uniforme en J m. • Si X = X1, . . . , Xm es uno de estos puntos, su p.d.f. tiene la forma siguiente (donde x = (x1, . . . , xm)):  fX (x) =

1 si 0 ≤ xi ≤ 1 para todo i = 1, . . . , m, 0 si no.

Esto corresponde a tener distribuciones U (0, 1) independientes para cada una de las componentes del vector X; por lo tanto, para sortear un Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

12

valor del mismo, alcanza con sortear de manera uniforme cada uno de los valores Xi. • El siguiente seudoc´ odigo describe este m´etodo de Monte Carlo. Procedimiento MonteCarlo-Volumen (region R, entero n) Entrada: region R, n tama˜ no de la muestra, ¯ ¯ Salida: λ(R) estimaci´ on del volumen, V [λ(R)] estimaci´on de la varianza 1. S = 0; /* Inicializaci´ on */ 2. For j = 1, . . . , n do 2.1 Sortear X(j) con distribuci´ on uniforme en J m; 2.2 If X(j) ∈ R then φ(X(j)) = 1 else φ(X(j)) = 0; 2.3 S = S + φ(X(j)); /* Acumular*/ ¯ 3. λ(R) = S/n; /*Estimador puntual de λ(R)*/ ¯ 4. V [λ(R)] = (S/n)(1 − S/n)/(n − 1); /*Estimador puntual de la
¯ ¯ varianza Var λ(R) del estimador puntual λ(R)*/

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

13

• Se puede ver que en el seudoc´ odigo precedente utilizamos una f´ormula para estimar la varianza distinta de la discutida en la Unidad 02. Esta nueva f´ ormula se basa en propiedades de v.a. de Bernoulli (v.a. discretas que toman valores 0 y 1 exclusivamente). Para el problema de estimar vol´ umenes (y en general para los problemas donde se acumulan v.a. Bernoulli), ambas f´ ormulas son equivalentes. • Derivamos la f´ ormula empleada ahora. Comenzamos observando que φ(X(j)) es una v.a. que toma dos valores, 0 o 1; la probabilidad de que φ(X(j)) = 1 es la probabilidad de que X(j) ∈ R, que es igual a λ(R), y la probabilidad de que valga 0 es su complemento, 1 − λ(R). • Por lo tanto, φ(X(j)) es una v.a. e Bernouilli con par´ametro λ(R), y su varianza se puede calcular facilmente, resultando igual a λ(R)(1 − λ(R)). • Usando la independencia entre los X(j), y propiedades de la varianza de una combinaci´ on lineal de v.a. independientes, sabemos que el Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

14

¯ estimador λ(R) = S/n tendr´a varianza n   X
1 1 (j) ¯ Var λ(R) = 2 Var (S) = 2 Var φ(X ) = n n j=1

=

1 1 nλ(R)(1 − λ(R)) = λ(R)(1 − λ(R)). 2 n n

• Aunque no conocemos el valor exacto de λ(R)(1 − λ(R)) , podemos estimarlo. Para esto, observamos que
1
1 2 E ((S/n)(1 − S/n)) = E (S/n) − S /n ) = E (S) − 2 E S = n n 2

2

1 1 = E (S)− 2 (Var (S)+(E (S))2) = λ(R)−λ(R)(1−λ(R)/)n−λ(R)2 = n n (n − 1)λ(R) − (n − 1)λ(R)2 (n − 1)λ(R)(1 − λ(R)) = . = n n Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

15

Por lo tanto, 



λ(R)(1 − λ(R)) (S/n)(1 − S/n) =E , n (n − 1) y entonces (S/n)(1 − S/n)/(n − 1) es un estimador insesgado de la varianza.

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

16

Ejemplo Problema: Supongamos que tenemos un c´ırculo con centro en el punto (x1, x2) = (0.3, 0.2) y radio 0.4, y queremos calcular la superficie de la intersecci´ on de este c´ırculo con el cuadrante positivo (es decir, x1 ≥ 0, x2 ≥ 0). 1

r=0.4

(0.3, 0.2) 0

1

Figure 2: Ejemplo de una superficie a estimar. Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

17

Si aplicamos Monte Carlo, sortearemos una cantidad n de puntos independientes y uniformemente distribuidos en (0, 1) × (0, 1), y usaremos como estimador del ´area a calcular la cantidad de puntos que cayeron dentro de la misma sobre el total de puntos sorteados. La siguiente tabla muestra el resultado de n = 10 sorteos independientes. Dado que de 10 sorteos s´ olo 3 puntos cayeron en la figura R cuya ´area se ¯ desea calcular, nuestra estimaci´ on es λ(R) = 0.30. La estimaci´on de la ¯ ¯ varianza es λ(R)(1 − λ(R))/(n − 1) = (0.3)(0.7)/9 = 0.0233, y la estimaci´ on de la desviaci´ on est´andar es la ra´ız cuadrada de este valor, aproximadamente 0.153.

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

18

x1 0.12 0.85 0.25 0.59 0.45 0.44 0.25 0.09 0.87 0.13

x2 0.94 0.92 0.73 0.43 0.53 0.91 0.20 0.57 0.37 0.73

Pertenece a R 0 0 0 1 1 0 1 0 0 0

Table 1: Una simulaci´ on detallada con n = 10 experimentos.

La figura siguiente muestra la ubicaci´ on de los 10 puntos sorteados. Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

19

1

0

1

Figure 3: Ejemplo de sorteo de 10 puntos uniformemente distribuidos en (0, 1) × (0, 1).

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

20

Una implementaci´ on C de Monte Carlo para resolver este problema se encuentra disponible en la siguiente direcci´ on: http://www.fing.edu.uy/inco/cursos/mmc/codigoC/circulo.c

Dicho c´ odigo genera puntos independientes y uniformemente distribuidos en (0, 1) × (0, 1), y controla si su distancia al centro del c´ırculo es menor o mayor que el radio del mismo (lo que equivale a preguntar si est´an en su interior o en su exterior).

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

21

Preguntas para auto-estudio • ¿Cu´al es la diferencia entre un estimador puntual y un estimador de intervalo? • Si tenemos una secuencia de regiones de igual volumen pero dimensi´on m creciente, y queremos estimar este volumen con una cierta precisi´on prefijada : ¿cu´al es el comportamiento computacional de un m´etodo de aproximaci´ on cl´asico y de un m´etodo Monte Carlo en funci´on de m? • Dar el seudoc´ odigo del M´etodo de Monte Carlo para calcular la medida de una regi´ on arbitraria incluida dentro de J m = [0, 1]m. • Explicar c´ omo se calcula la varianza en el caso de experimentos basados en sumas de variables aleatorias de Bernoulli.

Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

22

Ejercicio Problema: se desea estimar el volumen de una region R de [0, 1]4 definida por todos los puntos de la hiper-esfera de centro (0.3, 0.2, 0.7, 0.8) y radio 0.50, que adem´as cumplan las restricciones siguientes: 5x1 + 10x2 ≤ 4; x3 + x4 ≤ 1; 2x1 − 3x2 ≥ 0. Ejercicio 3.1: • Parte a: implementar un programa que reciba como par´ametro la cantidad de replicaciones n a realizar, y emplee Monte Carlo para calcular (e imprimir) la estimaci´ on del volumen de R, y la desviaci´on est´andar de este estimador. Incluir c´ odigo para calcular el tiempo de c´alculo empleado por el programa. Utilizar el programa con n = 106 para estimar el volumen de R. • (en la sesi´ on 5 se continuar´a este ejercicio). Fecha entrega: ver calendario en pagina web. Curso “M´etodos de Monte Carlo” - Facultad de Ingenier´ıa, Universidad de la Rep´ ublica (Uruguay)

23