Ejemplo del uso de la técnica Montecarlo en una función y su cambio de variable

Este es un ejemplo concreto del uso de técnica Montecarlo con una función y su cambio de variable.  la Integral original es:



El cambio de variable es:
 

De modo que la integral después de realizar el cambio de variable es:



Cuando graficamos la función original y la del cambio de variable, ver la Fig 1, encontramos que el cambio de variable permite obtener una función mucho más suave en el intervalo de integración.

Figura 1. Gráficas de funciones: original y con cambio de variable. La función original presenta más estructura, mientras que el cambio de variable muestra una función casi constante en el intervalo de trabajo

Ahora bien, cuando aplicamos la técnica Montecarlo varias veces a las dos funciones, observamos que las funciones convergen a un valor de integral, ver Fig. 2. Es decir, muestran estabilidad para presentar una respuesta única en un intervalo de confianza. Adicionalmente, encontramos que la función con cambio de variable converge mucho más rápido que la función original, lo cual evidencia la importancia de realizar un adecuado cambio de variable.

Fig. 2.  Graficas de los cálculos de integral por método Montecarlo (con 100 interacciones) de la función original y la de cambio de variable. Note que alcanza más rápido la estabilidad la integral con el cambio de variable.

Valor de la integral original 0.4520, con un error porcentual asociado 0.0143%. En contraste, la integral con el cambio de variable presenta un valor promedio de -0.2948 con un error porcentual = 0.0077% (Mucho menor es el error porque el proceso alcanza rápidamente la estabilidad).

% incia rutina, begins script
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function montecarloevalua
%Rutina para graficar una integral por el metodo de MonteCarlo,
%calcular el tiempo de ejecución, y contar el con la estabilidad
clc % parte donde se limpia la pantalla y la memoria
close all
clear
a=exp(-pi^2); % límite inferior de la integral
b=1; % límite superior de la integral
x1=linspace(a,b,1000); % vector donde ser realiza la integral
y1=1/2*(sin(exp(-x1.^2)))./exp(-x1.^2); % función a integrar
plot(x1,y1) % gráfica de la función original a integrar
figure
avg=[0 0]; %donde se va a contener y
%inicia método montecarlo
tic % inicia contador del reloj
iteracion=100;% numero de veces que se repite la interacción MonteCarlo,
for i=1:10 % Número de graficas para saber la estabilidad del cálculo
for N = 1:iteracion; %iteracion dentro del montecarolo
x1 = a + (b-a)*(rand(1000,1));
term= sum([(1/2*(sin(exp(-x1.^2)))./exp(-x1.^2)) (1/2*(sin(exp(-x1.^2)))./exp(-x1.^2)).^2])/1000;
avg= ((N-1)*avg + term)/N;
z(N)=(b-a)*avg(1);
end
%fin estricto del proceso MonteCarlo.
prom(i) = (b-a)*avg(1); % último punto calculado del proceso MonteCarlo
N = 1000*(1:iteracion); % vector para graficar
V = [z'];
plot(N,V);
hold on;
xlabel('N')
end
toc % finaliza el contador del reloj
promedio_int = mean(prom) %promedio de las integrales calculadas
error = std(prom); %error asociado al promedio de la integral
error_porcentual = 100*error/promedio_int % error porcentual
% termina la rutina, end script
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%inicia método montecarlo tic % inicia contador del reloj iteracion=100;% numero de veces que se repite la interacción MonteCarlo,
for i=1:10 % Número de graficas para saber la estabilidad del cálculo
for N = 1:iteracion; %iteracion dentro del montecarolo x1 = a + (b-a)*(rand(1000,1)); term= sum([(1/2*(sin(exp(-x1.^2)))./exp(-x1.^2)) (1/2*(sin(exp(-x1.^2)))./exp(-x1.^2)).^2])/1000; avg= ((N-1)*avg + term)/N; z(N)=(b-a)*avg(1); end %fin estricto del proceso MonteCarlo. prom(i) = (b-a)*avg(1); % último punto calculado del proceso MonteCarlo N = 1000*(1:iteracion); % vector para graficar
V = [z']; plot(N,V); hold on; xlabel('N') end toc % finaliza el contador del reloj promedio_int = mean(prom) %promedio de las integrales calculadas error = std(prom); %error asociado al promedio de la integral error_porcentual = 100*error/promedio_int % error porcentual % termina la rutina, end script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Sustituyendo en la rutina las siguientes líneas

a=exp(-exp(-2*pi^2));
b=1/exp(1);
y1=-1/4*(sin(x1)./((x1.^2).*sqrt(-log(x1))));

Ejercicios:

1) Proponga otras funciones y su cambio de variable para repetir este ejercicio
2) Cómo se puede hacer más rápido esta algoritmo
3) Modifique el algoritmo para contar con una función que se alimente con la función a integrar y los limites de integración