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

Presentación: Introducción a la técnica Montecarlo con 2 ejemplos en Matlab

Montecarlo en matlab  

De acuerdo con la Wikipedia, el uso de los métodos de Montecarlo como herramienta de investigación, proviene del trabajo realizado en el desarrollo de la bomba atómica durante la segunda guerra mundial en el Laboratorio Nacional de Los Álamos en EE.UU.

Los métodos Montecarlo son una herramienta computacional que utilizan las propiedades deterministas de un sistema en conjunto con variables aleatorias para aproximar un resultado, es decir, se basa en un muestreo repetitivo aleatorio. En general se utiliza en problemas donde la evaluación analítica no es posible, o resulta muy complicada; en el estudio de sistemas con un gran número de grados de libertad acoplados, como son los fluidos, estructuras celulares, etc.

Más en general, los métodos de Montecarlo pueden ser útiles para modelar fenómenos con entradas de incertidumbre significantes, como en el cálculo de riesgo en los negocios.

¿Dónde más haz escuchado usar el método Montecarlo?

20 puntos extra: por mejorar interface y programa de captura de imagenes

Interface de video y analisis de una imagen capturada 

Este es un script  para obtener el video (con su archivo fig), la imagen congelada y estadística de esa imagen, funciona con la interface de una cámara web. Funciona en base a la captura de video, cuyo código se encuentra en nuestra práctica de grabación y manipulación de video por medio de Matlab.

En primer equipo que entregue este código con una función en la interface para variar y fijar la sensibilidad de la cámara web. Obtiene hasta 20 puntos (para cada uno de los miembros del equipo) en las evaluaciones de la clase.

El máximo de integrantes por equipo es 3.

La fecha límite para hacer el código es el 30 de noviembre.

Preguntas sobre este proyecto sólo serán contestadas en este post y en la clase presencial.


Links útilies:


Ayuda de Matlab Adq. Toolbox

Cómo simular en Matlab un patrón de difracción de una apertura circular

Apertura circular  y su patrón de difracción simulados con Matlab 

El siguiente guion es para usarse en Matlab y realizar una simulación de un patrón de difracción. Puede ser útil para motivarlos a estudiar más este tema:

%inicia script
%% script para simular un patrón de difración de una apertura circular
%% mediante la transfotmada discreta y rapida de Fourier
%
%% por Vicente Torres (18/8/2011), en base al script de Luis Mex

%limpieza de memoria
clear;
clc;
close all
%acondicionamiento
n=2^10;
M=zeros(n); % matriz base donde, que representa la pantalla, oscura
I=1:n;
x=I-n/2;
y=n/2-I;
[x,y]=meshgrid(x,y);
R=10;
a=(x.^2+y.^2<=R.^2); % acondicionando el área de la apertura,
% estos seran indices de la matriz
M(a)=1; % la apertura deja pasar la luz, que se representa con 1
figure
subplot(1,2,1)
imagesc(M),colormap( [0 0 0;1 1 1]),title('\bf abertura circular'),axis image;
D=fft2(M); % Calculo del patron de difraccion
D1=fftshift(D);
subplot(1,2,2)
imagesc(abs(D1)),colormap(hot),axis image, title('\bf patrón de difracción');
%fin del script

En Matlab central pueden encontrar otros ejemplos de patrones de difracción, basado en la ayuda de Matlab.

Preguntas para pensar:

La forma del patrón de difracción en campo lejano  parece concordar con el experimento, como se puede comprobar que su perfil de intensidad corresponde adecuadamente con la teoría o con el experimento. 

Ejercicios:
1) En base a este ejemplo obtenga el patrón de difracción de una rendija simple.
2) Obtenga la simulación de un patrón de difracción  de una rendija cuadrada

Cómo hacer un diagrama de barras estadístico con Matlab

Existen varias posibilidades para representar diagramas de barras. Por ejemplo, los datos se escriben como un vector y son introducidos a Matlab:

>> x = [10 2 3 4 18 20 15];

Entonces, podemos usar los comandos bar, barh, bar3 y bar3h para generar los gráficos. Usamos como apoyo el comando subplot para ver en una sola imagen todos estos resultados.

>> subplot(2,2,1), bar(x), title('Barras verticales')

>> subplot(2,2,2), barh(x), title('Barras horizontales')

>> subplot(2,2,3), bar3(x), title('Barras verticales 3D')

>> subplot(2,2,4), bar3h(x), title('Barras horizontales 3D')

El resultado lo podemos ver en la siguiente figura:

Diagrama de barras

Por supuesto,  las graficas 3D se pueden modificar empleando el comando rotate3d.

Además, los datos pueden estar agrupados. Como ejemplo, consideremos los siguientes datos

x = [1 2 3; 4 3 6; 10 9 8; 4 2 7; 12 10 7 ];

Utilizando los mismos comandos que antes, con modificaciones como: bar3(x,'group') y bar3(x,'steack'). Obtenemos la siguiente figura.

Diagrama de barras apiladas y agrupados 3D

Preguntas para pensar.

En que casos prácticos son útiles matrices de MxN mostradas como diagrama de barras.

Ejercicios:

Empleando la función rand genere un vector de 1000 datos,  agrupe los valores entre 0-10, 10-20, 20-30 y así.  ahora haga el diagrama de barraras con la agrupación, por lo que solo aparecernos 10 barras.