montecarlos
A. Casares M
4 de mayo de 2018
Ejercicios con Métodos de Montecarlo:
Son ejemplos simples con el propósito de ilustrar su forma de aplicación.
Referencia 1: Montecarlo Methods, Jonathan Pengelly, February 26, 2002.
Referencia 2: Wikipedia: Buffon’s Needle.
Referencia 3: www.rpubs.com/acasares/246542.
1. Presentación:
Todos los ejemplos están orientados a estimar el valor de \(\pi\).
Cada ejemplo se ha desarrollado en un guión MATLAB independiente. Se expone sólo la parte fundamental de cada guión.
El resultado de cada caso se presenta por lo general en un par de gráficos enlazados, que muestran, a más de la estimación calculada de \(\pi\), la evolución de los errores relativos y el tiempo de ejecución para cada orden de tamaño de la muestra de números aleatorios.
2. Ejemplo 1: \(\pi\) como área de un círculo.
Evalúa el área A de un círculo de radio r dado: A = pi*r^2.
Implementa la técnica de “Crude Montecarlo”, en el 1er cuadrante, forma que podría considerarse “estratificada” (la estratificación consistiría en usar sólo el 1er cuadrante y extrapolar luego el valor de área obtenido a toda la figura).
También estima el error relativo en la respuesta obtenida e informa el tiempo de ejecución.
% Datos:
r = 2; % radio del círculo (puede cambiarse).
fc = 0.017453292519943; % Factor de conversión a radianes
np = 7;
% Aplicación para varios tamaños de muestra:
for in = 1:np
N = 10^in; % Tamaño de la muestra
tic
sf = 0;
for i = 1:N
% Genera un punto del círculo en el 1er cuadrante:
x = r*rand; fx = sqrt(r^2-x^2);
sf = sf + fx; % Suma las ordenadas
end
end
% Estimación de pi en este método crudo, según la fórmula (1):
% El factor 4 se debe a que son 4 cuadrantes
epi = 4*sf/(r*N);
% Tiempo de cálculo:
toc
errel(in)=(epi-pi)/pi;
elapstime(in)=toc;  |
|---|
| Figura 1 |
Conclusión de ejemplo 1:
El método converge en forma lenta. Aún con 10 millones de puntos sólo ha alcanzado una precisión de 4 cifras significativas correctas. A partir de 10000 puntos la reducción del error relativo ya no es apreciable - aunque parece lineal - pero en cambio el tiempo de ejecución comienza a crecer notablemente desde los 100000 puntos en adelante.
3. Ejemplo 2: \(\pi\) como área de un círculo inscrito en un cuadrado.
Evalúa el área de un círculo inscrito en un cuadrado, mediante la técnica de “Aceptación-Rechazo”, en el 1er cuadrante.
También en este caso el procedimiento podría considerarse “estratificado”, al igual que en el ejemplo 1 y por la misma razón.
Estima el error relativo en la respuesta obtenida e informa el tiempo de ejecución.
% Datos:
r = 2; % radio del círculo (puede cambiarse)
fc = 0.017453292519943; % Factor de conversión a radianes
np = 7;
for in = 1:np
N = 10^in; % Tamaño de la muestra
% Acceptation-Rejection Monte Carlo:
tic
h = 0;
for i = 1:N
% Genera las coordenadas en el 1er cuadrante dentro del cuadrado:
x = r*rand; y = r*rand;
d = sqrt(x^2+y^2);
% Chequea si cae o no dentro del círculo:
if d <=r % Sí cae:
h = h+1;
end
end
% Estimación de pi:
% Como el área C del círculo es pi*r^2, la del cuadrado A = 4*r^2
% y la relación C/A debe aproximarse por (4*h)/(4*N), podemos estimar:
epi = 4*h/N;
toc
errel(in)=(epi-pi)/pi;
elapstime(in)=toc;
end  |
|---|
| Figura 2 |
Conclusión de ejemplo 2:
El método converge en forma más lenta y no uniforme. Aún con 10 millones de puntos sólo ha alcanzado una precisión de 3 cifras significativas correctas. La reducción del error relativo no es uniforme, pero en cambio el tiempo de ejecución comienza a crecer notablemente desde los 100000 puntos en adelante.
4. Ejemplo 3: \(\pi\) como área de un círculo circunscrito a un cuadrado.
Evalúa el área de un círculo circunscrito a un cuadrado, mediante la técnica de “Aceptación-Rechazo”.
Estima el error relativo en la respuesta obtenida e informa el tiempo de ejecución.
% Datos:
r = 2; % radio del círculo (puede cambiarse)
lm = r*sqrt(2)/2;
fc = 0.017453292519943; % Factor de conversión a radianes
np = 7;
for in = 1:np
N = 10^in; % Tamaño de la muestra
tic
h = 0;
for i = 1:N
% Genera las coordenadas polares de un punto dentro del círculo:
rp = r*rand();
theta = rand()*360*fc; % En radianes
% Chequea si cae o no dentro del cuadrado:
x = rp*cos(theta); y = rp*sin(theta);
if abs(x) <= lm && abs(y) <= lm % Sí cae:
h = h+1;
end
end
% Estimación de pi:
% Como el área C del círculo es pi*r^2, la del cuadrado A = 2*r^2
% y la relación C/A debe aproximarse a N/h, podemos estimar:
epi = 2*N/h;
toc
errel(in)=(epi-pi)/pi;
elapstime(in)=toc;
end |
|---|
| Figura 3 |
Conclusión de ejemplo 3:
El método no converge al valor de \(\pi\) (queda como ejercicio determinar qué está fallando y corregirlo). El módulo del error relativo se estabiliza alrededor de 20% para las tres muestras mayores. El tiempo de ejecución comienza a crecer notablemente desde la muestra con cien mil puntos.
Puede hallar una explicación sobre el comportamiento equivocado de este procedimiento en la NOTA final de este documento.
4. Ejemplo 4: \(\pi\) como probabilidad al arrojar repetidamente la aguja de Buffon.
Evalúa la probabilidad de que la aguja de Buffon corte las líneas de un entablado vertical.
Para la teoría, consulte la referencia 2 para el caso de aguja corta.
Estima el error relativo en la respuesta obtenida e informa el tiempo de ejecución.
% Montecarlo4
% Guión para mostrar el método de Montecarlo en un ejemplo clásico:
% la aguja de Buffon, utilizada para evaluar $\pi$.
% Referencia 1: tutorial "Montecarlo Methods" de Jonathan Pengelly (2002)
% Implementa la técnica denominada "Aceptación-Rechazo", que resulta
% la apropiada para este ejemplo.
np = 7; errel = []; elapstime = [];
fprintf('Estimación de pi usando la aguja de Buffon\n')
fprintf('==========================================\n')
for in = 1:np
% Datos:
N = 10^in; % Tamaño de la muestra
L = 0.8; % Longitud de la aguja
t = 1; % Separación de la malla (vertical)
if L > t
fprintf('En este programa la aguja no debe medir más de t = %3.2f\n',t)
return
end
fc = 0.017453292519943; % Factor de conversión a radianes
diary montecarlo.txt
% Sea T un entablado formado por tablas verticales contiguas.
% Sean T0, T1, T2 tres "tablas" contiguas de T, de igual anchura = t.
% Pongamos el origen de coordenadas en algún punto del borde izquierdo de T1.
% Entonces, la abscisa izquierda de T1 = t1 = 0,y las abscisas izquierdas
% de T0 y T2 son, respectivamente, t0 = -t y t2 = t, mientras el borde
% derecho de T2 tiene abscisa t3 = 2t.
% Se generará aleatoriamente un extremo de la aguja, de abscisa Px, tal que
% 0 <= Px <= t3. Luego se generará igualmente el ángulo theta entre la aguja
% y el eje horizontal x, tal que 0 <= theta <= 360.
% Con esto se puede calcular la abscisa Rx del otro extremo de la aguja.
% Sea Z1 la tabla en que quedó Px y Z2 aquélla en la que está Rx.
% Entonces, la aguja habrá cruzado un borde si y sólo si Z1 y Z2 son diferentes.
% Según la referencia 2, $\pi$ es aproximadamente 2L.N/(t.h).
% Técnica de Aceptación-Rechazo:
tic
h = 0;
for i = 1:N
% Genera la abscisa de un extremo de la aguja:
% Es mucho más rápido que usar random('unif',0,2*t):
Px = rand*2*t;
% Genera el ángulo entre la dirección de la aguja y la horizontal
theta = fc*rand*360;
%fprintf('Abscisa extrema: %6.4f, Angle = %7.4f\n',Px,theta)
% Calcula el otro extremo de la aguja:
Rx = Px+L*cos(theta);
% Encuentra las zonas (-1,0 o 1) de ambos extremos:
Z1 = floor(Px/t); Z2 = floor(Rx/t);
% Chequea si cruza o no una vertical de la malla:
if Z1 ~= Z2 % Sí la cruza:
h = h+1;
end
end
% Estimación de pi:
epi = 2*L*N/(t*h);
toc
errel(in)=(epi-pi)/pi;
elapstime(in)=toc;
end
subplot(211), plot(1:np,abs(errel))
title('Error relativo')
grid on
subplot(212), plot(1:np,elapstime)
strepi = sprintf('%10.8f',epi);
title('Tiempo de ejecución (sec)')
grid on
xlabel('log_{10}N')
suptitle([{'Cálculo de pi usando la aguja de Buffon:'},...
{['\pi =',strepi]}]) |
|---|
| Figura 4 |
Conclusión de ejemplo 4:
El método converge al valor de \(\pi\) en forma más rápida que en los casos en que se han usado relaciones de áreas. Los módulos de errores relativos van disminuyendo nítidamente y se llega a una precisión de 5 dígitos. El tiempo de ejecución comienza a crecer notablemente desde la muestra con cien mil puntos.
6. Ejemplo 5: \(\pi\) como semiperímetro de un polígono regular de N lados encerrado entre dos cotas.
Esta no es una aplicación de Montecarlo, sino una estimación de \(\pi\) usando intervalos. Se debe a Arquímedes, pero se han ajustado las fórmulas para evitar los problemas de inestabilidad numérica que se presentan al utilizar punto flotante.
La teoría básica y las fórmulas pueden consultarse en la referencia 3.
EL valor de \(\pi\) está encerrado en un intervalo [S1,S2], donde S1 es el semiperímetro de un polígono regular inscrito y S2 el de uno circunscrito a un círculo de radio 1.
Conforme el número de lados de estos polígonos crece (Arquímedes duplica el número de lados en cada iteración), el radio del intervalo decrece, y así la aproximación a \(\pi\) mejora.
% Programa para calcular pi por Arquímedes (primer ejemplo conocido de
% analisis de intervalos).
% Algoritmo estable, que evita la sustracción cancelativa.
% La cota inferior se calcula por el semiperímetro del polígono inscrito.
% La cota superior sale del semiperímetro del circunscrito.
diary arquimedes
fprintf(' CALCULO DE PI: FORMULA ESTABLE\n')
fprintf(' =================================\n\n')
l=1;n=6; % Comienza desde el hexágono inscrito de lado 1.
si=3;sc = 6/sqrt(3);
fprintf(' n lado semip.inscr semip.circuns.\n');
fprintf('------------------------------------------------------------------------------\n')
fprintf('%10ld%15.6e%25.15e%26.15e\n',n,l,si,sc);
tic
it = 0;
while (sc - si)/si > eps
it = it+1;
l = l/sqrt(2 + sqrt(4 - l*l)); % Fórmula estable para el nuevo lado */
si = n * l; % semiperímetro inscrito
sc = 2 * si /sqrt(4 - l*l); % semiperímetro circunscrito
n = 2 * n; % número de lados actual
% Salida en precisión doble:
fprintf('%10ld%15.6e%25.15e%26.15e\n',n,l,si,sc);
relerr(it) = (sc-si)/(2*pi);
elapstime(it)=toc;
end
toc
epi = (si+sc)/2; strepi = sprintf('%21.15e',epi);
fprintf('Valor de pi estimado:%21.15e\n',epi)
fprintf('|Error relativo|: %10.4e\n',relerr(it))
diary off
figure
plot(elapstime,abs(relerr))
grid on
ylabel('| Error relativo | ')
xlabel('Tiempo empleado (sec.)')
title([{'Cálculo de pi usando intervalos (Arquímedes)'},...)
{['\pi =',strepi]}])
Resultados de la ejecución:
CALCULO DE PI: FORMULA ESTABLE
=================================
n lado semip.inscr semip.circuns.
------------------------------------------------------------------------------
6 1.000000e+00 3.000000000000000e+00 3.464101615137755e+00
12 5.176381e-01 3.105828541230249e+00 3.215390309173472e+00
24 2.610524e-01 3.132628613281238e+00 3.159659942097500e+00
48 1.308063e-01 3.139350203046867e+00 3.146086215131434e+00
96 6.543817e-02 3.141031950890509e+00 3.142714599645368e+00
192 3.272346e-02 3.141452472285462e+00 3.141873049979823e+00
384 1.636228e-02 3.141557607911858e+00 3.141662747056849e+00
768 8.181208e-03 3.141583892148318e+00 3.141610176604689e+00
1536 4.090613e-03 3.141590463228050e+00 3.141597034321526e+00
3072 2.045307e-03 3.141592105999271e+00 3.141593748771351e+00
6144 1.022654e-03 3.141592516692156e+00 3.141592927385096e+00
12288 5.113269e-04 3.141592619365383e+00 3.141592722038613e+00
24576 2.556635e-04 3.141592645033690e+00 3.141592670701997e+00
49152 1.278317e-04 3.141592651450766e+00 3.141592657867843e+00
98304 6.391587e-05 3.141592653055036e+00 3.141592654659305e+00
196608 3.195793e-05 3.141592653456103e+00 3.141592653857171e+00
393216 1.597897e-05 3.141592653556370e+00 3.141592653656637e+00
786432 7.989483e-06 3.141592653581437e+00 3.141592653606504e+00
1572864 3.994742e-06 3.141592653587703e+00 3.141592653593970e+00
3145728 1.997371e-06 3.141592653589270e+00 3.141592653590837e+00
6291456 9.986854e-07 3.141592653589662e+00 3.141592653590053e+00
12582912 4.993427e-07 3.141592653589759e+00 3.141592653589858e+00
25165824 2.496714e-07 3.141592653589784e+00 3.141592653589809e+00
50331648 1.248357e-07 3.141592653589791e+00 3.141592653589797e+00
100663296 6.241784e-08 3.141592653589792e+00 3.141592653589794e+00
201326592 3.120892e-08 3.141592653589793e+00 3.141592653589794e+00
Elapsed time is 0.016447 seconds.
Valor de pi estimado:3.141592653589793e+00
|Error relativo|: 7.0679e-17 |
|---|
| Figura 5 |
Conclusión de ejemplo 5:
La convergencia es continua: el módulo del error relativo decrece uniforme y acusadamente desde el comienzo.
El orden de convergencia es mucho mayor que en los casos anteriores, como se puede ver tanto en el comportamiento del error relativo como en la precisión alcanzada de 16 dígitos (la mayor alcanzable con la precisión doble del sistema de punto flotante, utilizada por MATLAB) en un tiempo menor que el de los casos anteriores.
NOTA. Sobre la no convergencia a \(\pi\) del ejemplo 3:
El problema se debe a que en el método de Montecarlo la distribución de los puntos en la superficie que se va a estimar debe ser uniforme, y en la forma en que el guión está programado, en la cual cada punto se genera por sus dos coordenadas polares, si bien el radio vector y el ángulo tienen sendas distribuciones uniformes, los puntos en sí no la tienen.
Una explicación intuitiva de la razón de este comportamiento puede ser que se está generando, no una malla rectangular de puntos, sino una malla radial, en la cual, para cada radio, el número de ángulos que se consideran (y por ende, el número de puntos generados) es siempre el mismo (por la distribución uniforme del ángulo). Y como la circunferencia va creciendo proporcionalmente al radio, mientras mayor sea éste, más distanciados estarán los puntos sobre ella (como sucede en el segmento circular), mientras que, en las vecindades del centro, donde los radios son menores, su densidad de puntos será más grande.
Para rectificar este comportamiento, el radio polar no debe calcularse con la fórmula que se usó en el ejemplo 3: (1) rp = r.rand(), sino con (2): rp = r.sqrt(rand()). Es decir, usando, no directamente cada número aleatorio uniforme, sino su raíz cuadrada.
En el gráfico siguiente se nota claramente la diferencia entre las distribuciones de puntos obtenidas con las fórmulas (1) y (2), para una muestra de diez mil puntos:
 |
|---|
| Figura 6 |
Con ese cambio en el guión del ejemplo 3 se obtiene el resultado resumido en la figura 7:
 |
|---|
| Figura 7 |
En conclusión, se ha corregido el comportamiento anterior, y ahora el método converge a \(\pi\) en forma similar a los ejemplos 1 y 2. Con 10 millones de puntos alcanza una precisión de 4 cifras significativas correctas. La reducción del error relativo es uniforme desde la muestra de 10000 puntos en adelante, y el tiempo de ejecución comienza a crecer notablemente desde los 100000 puntos.
Se hace hincapié entonces, a partir de este comportamiento, que en este tipo de aplicaciones del método de Montecarlo es muy importante asegurarse de que la distribución de los puntos generados sea uniforme sobre la superficie que se desea estimar.