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.

Código Matlab para acondicionar y filtrar señales de electrocardiogramas.

Sergey Chernenko nos muestra, explica y regala el código Matlab para detectar una característica importante de las señales provenientes de un electrocardiógrafo: el pico R de un electrocardiograma.

Para estudiar este excelente ejemplo de procesamiento de señales en ciencias medicas, lo mejor es bajar el archivo zip, con el código original.

Hay que estudiarlo para comentarlo en la clase.

Como nos comenta el mismo Sergey en su sitio Web, la tarea básica en el procesamiento de electro-cardiogramas (ECG) es la detección de los picos R. Sin embargo, debido a la respiración del paciente, la detección es complicada; pues se presenta en los ECGs artefactos que complican la detección de los picos R. Por ejemplo, las distancias irregulares entre los picos, la forma irregular en los picos, la presencia de componentes de baja frecuencia en los ECGs.

Es por ello que se deben contar con varias etapas de procesamiento que reduzcan la influencia de estos factores. Sergey nos muestra un caso donde una señal ECG original se le eliminan las frecuencias bajas por medio de FFT (Fast-Fourier-Transform). En la segunda etapa, se usa un filtro pasa banda, el cual localiza los máximos significativos en la señal. En la etapa final, se eliminan los datos más alejados del promedio. Este proceso permite la detección de los picos R.

Código para obtener la primera derivada por coeficientes de una función

Como se puede consultar en el libro: transformada de Fourier para peatones, podemos definir la derivada Yk como la “primera diferencia central” de la función F:

La cual tiene una función de transferencia, que no presenta problemas de corrimientos de fase:

Sin embargo, presenta un efecto de filtro muy pronunciado a frecuencias altas.

El siguiente código matlab aplicado a la función sin^2(x) con un paso de 0.1 funciona muy bien. En las figuras de salida del script se observan la función original, la derivada teórica y la derivada obtenida. Observe que presenta un comportamiento muy similar las dos derivadas. Efectivamente, cuando comparamos la derivada obtenida sobre la deriva teórica obtenemos en la gráfica una recta identidad. Por todo lo anterior, afirmamos que es un buen algoritmo para operar la derivada.

%%Script que calcúla numéricamente la primera derivada de la función:

%%sin(x)^2.

clear all

close all

clc

%Definiendo intervalos

x=0:0.01:1;

%función original:

funcion=(sin(x)).^2;

%preparando parámetros:

s=length(funcion); %longitud de la funcion

Y=zeros(1,s); % espacio de memoria para la derivada numerica

% Algoritmo de primera derivada

for k=2:s-1

Y(k)=(funcion(k+1)-funcion(k-1))/(2*0.01);

end

% Graficación de funciones teoricas y su compraración con la estimación

% numerica de la primer derivada

figure

subplot(2,1,1)

hold on

plot(x,funcion,'r') % función original

plot(x,2*sin(x).*cos(x), '-b') % derivada teorica de la funcion original

plot(x,Y,'ok') % derivada numerica obtenida

axis([0,1.04,0,1.1]);

xlabel('Vector <>')

ylabel('Funciones utilizadas')

text(0.02,0.8,'Funciones Sin^2(x) (rojo), derivada teorica (azul), derivada numerica (negro)')

box on

hold off 

subplot(2,1,2)

plot(2*sin(x).*cos(x),Y,'-b')% comparacion entre funciones;

% derivada teorica contra la derivada numerica

xlabel('Función 2*sin(x).*cos(x), derivada teorica')

ylabel('Derivada numerica obtenida')

axis([0,1.04,0,1.1]);

grid

%%% fin del script 

Ejercicios:

Matlab tiene su comando/función para hacer la derrivada: diff, su algoritmo de trabajo es:


1) ¿Cómo se compara este algoritmo con el propuesto en el post?

2) Muestre ejemplos del uso de este comando

3) ¿Cuales son las limitaciones de esta forma de definir la derivada?

Enlaces relacionados:

Capítulo 2, Derivación numérica del curso de métodos numéricos de Agustín A. Rosas.

Comparación del desempeño de un filtro convolución y un filtro pasa bajas de promedios ponderados en una misma función/señal

Ahora mostramos dos ejemplos de filtros. La funcion a trabajar es:

La cual tiene un pico extendido en t=5 y dos pequeños picos de estructura en t=2. Pues bien la entrada a un instrumento  cuenta con  filtro descrito por:

Efectivamente, la señal filtrada por  convolución será la transforma de Fourier inversa de la multiplicación de las transformadas de Fourier de la señal y del filtro; es decir:

A continuación aplicaremos el filtro en código Matlab.

Filtro 1: Convolución.

Para comparar la señal y la salida filtrada, a ambas las normalizamos. Observamos que hay un corrimiento en el tiempo. Este efecto en los indices y en la amplitud son muy comunes cuando se usa la transformada de Fourier.

Adicionalmente, una vista cercana a la estructura de los dos picos pequeños muestra que la salida del filtro cuenta con estos picos más pequeños y casi sin desernir. Con todo, esta forma de presentar los datos muestra a detalle que el filtro es bueno sin ser excelente para eliminar la estructura de la curva principal .

A continuación, mostramos el código Matlab utilizado para hacer estas curvas.

 %%% Inicio código 1: Filtro convolución

 clc

clear

close all

x = 0:0.0001:10;

y = 4*((sin(x-2)).^2).*exp(-50*((x-2)).^2) + 0.5*( ((x-5).^2 +0.25).^(-.25) );

r = ones(size(x));

for k= 1:length(x)

    if x(k) <= 0.15

        r(k)= .50;

    elseif x(k) <= 0.30

        r(k)= 0.25;

    elseif x(k) <= 0.60

        r(k)= 0.10;

    else

        r(k)= 0;

    end

end

plot(x,r,  'LineWidth',3)

axis([-.2,3,-0.2,0.6])

Y= ifft( fft(y).*fft(r));

figure

hold on

plot(x, y/max(y), '--b','LineWidth',3)

plot(x,Y/max(Y), 'k', 'LineWidth',3)

axis([0,10,0.28,1.05])

hold off

%%% Fin código 1: Filtro convolución

Filtro 2: Pasa bajas con promedio ponderado

Ahora bien, utilizaremos el filtro pasa bajas (PB), que ya mostramos en otro post, en este caso usamos 6,0000 iteraciones. Dado que la estructura en la función puede representar una oscilación rápida aplicaremos el filtro PB.   Observamos en la salida que tanto el cambio en la amplitud como el corrimiento en los indices de los vectores de las curvas es casi minimo. Además de que se ve más disminuida la estructura en la función.

Este filtro PB presenta mejores resultados que el de convolución; esto es por el tipo de filtro r(t) utilizado. Sin embargo, el filtro PB  tarda mucho más en desarrollarse. Puedes  usar los comandos de Matlab TIC yTOC antes y despues de que inicien las rutinas para observar la velocidad de ejecución de cada función o escript.

Podemos asegurar que dependiendo que se busquemos filtrar, hay un filtro que nos funcionara mejor que otros, pero hay que saberlo encontrar :)

A continuacion el segundo código Matlab:

%%% Inicio código 2: filtro pasa bajas con promedio ponderado

clc

clear

close all

tic % inicia el contador de tiempo

t = 0:0.001:5*pi;

f = 4*((sin(t-2)).^2).*exp((-50)*(t-2).^2) +0.5*((t-5).^2+0.25).^(-0.25); %señal de trabajo

hold on

plot(t, f, 'r') % para visualizar la señal con ruido

for corrida = 1:60000 %iteracion del filtro

for k = 2:(length(t)) -1 ;

f(k) = 0.25*(f(k-1) + 2*f(k) + f(k+1)); %filtro

end

end

plot(t, f, '--b') %para visualizar la señal filtrada

hold off

toc %finaliza el contador del tiempo

%%% Fin código 2: filtro pasa bajas con promedio ponderado

Ejercicios: 

1) Muestre el grado de atenuación de las protuberancias después de pasar por el filtro.

2) Diseñe otro filtro de convolución que atenue mejor a los dos picos que forman protuberancias en la curva original

3) Calcule la salida del filtro de convolución con r(w)*I(w) en lugar de  I(w)*r(w). Describa sus observaciones.

Ejemplos de código para formar señales discretas en Matlab

El objeto más básico en Matlab es una matriz numérica con la posibilidad de almacenar números complejos. Por supuesto, los datos obtenidos en el estudio de señales y sistemas son muy bien representados en forma de matrices. En este post usaremos Matlab para la generación de señales elementales: cuadrada, triangular, entre otras.

El ToolBox de procesamiento de señales de Matlab posee una gran variedad de funciones para la generación de señales, estas señales requieren de una representación vectorial de la variable tiempo, de manera continua o discreta. Para realizar una simulación de un intervalo continuo, se usa un vector de valores discretos con un intervalo de muestreo muy pequeño.

Como vimos en post anteriores, el siguiente comando genera un vector llamado t de valores que representan la variable tiempo, con un intervalo de muestreo de 1 ms entre 0 y 1 segundo.

t = 0:0.001:1;

Después de creado el vector que representa la variable tiempo, es posible iniciar el desarrollo de alguna señal de interés.

En Matlab una señal discreta en el tiempo se representa exactamente, porque los valores de la señal son representados como los elementos de un vector. Sin embargo las señales de tiempo continuo en Matlab son tan solo aproximaciones. La aproximación consiste de un vector cuyos elementos son muestras de la verdadera señal de tiempo continuo. Cuando se usa esta técnica para la representación de señales continuas es importante escoger el intervalo de muestreo lo suficientemente pequeño para asegurar que las muestras capturan todos los detalles de la señal.

EJEMPLOS DE SEÑALES EN MATLAB 

La generación de señales periódicas tales como ondas cuadradas y triangulares es una actividad muy fácil de realizar en MATLAB.

1) SEÑAL CUADRADA

Consideremos primero la generación de una onda cuadrada de amplitud A, frecuencia fundamental w (medida en radianes por segundo) y ciclo útil rho. Recordemos que el ciclo útil es la fracción de cada periodo en donde la señal es positiva.

Para generar dicha señal se pueden escribir lo siguiente en la linea de comandos:

>> A = 1;

>> w = 10 * pi;

>> rho = 0.5;

>> t = 0:0.001:1;

>> sq = A*square(w*t+rho);

>> plot(t,sq);

En la segunda línea de este ejemplo, pi es una función interna de Matlab que calcula el número más cercano a la constante PI en formato de coma flotante. El último comando es usado para vizualizar la señal generada. El comando plot dibuja líneas conectando los valores sucesivos de la señal y así da la apariencia de una señal en tiempo continuo.

2) SEÑAL TRINGULAR 

Consideremos ahora la generación de una onda triangular de amplitud A, frecuencia fundamental w y ancho Wdt . El periodo de la onda triangular será T con el máximo valor de la señal ocurriendo en t = WT . El comando básico para generar esta señal es:

A * sawtooth(w * t + Wdt)

El resultado se puede observar en la gráfica a la izquierda

3) SEÑAL ESCALÓN

En Matlab, el comando ones(M, N) genera una matriz de unos de tamaño MxN, y el comando zeros(M, N) es una matriz de ceros del mismo tamaño. Se puede hacer uso de estas dos matrices para generar dos señales comúnmente usadas: la señal escalón y la señal impulso.

Una señal paso de amplitud uno, puede ser generada con el siguiente comando.

U = [zeros(1, 10), ones(1, 11)];

Para la versión continua creamos un vector que represente el tiempo el cual tenga muestras de un intervalo separados por valores muy pequeños. Los comandos y los resultados se muestran a continuación:

>> u=[zeros(1,10),ones(1,11);

>> t=-1:0.1:1;

>> plot(t,u)

Como se menciono anteriormente, una señal generada en Matlab es inherentemente de naturaleza discreta. Para visualizar una señal en tiempo discreto se puede hacer uso del comando stem. Específicamente stem(n, x), bosqueja los datos contenidos en el vector x como una señal de tiempo discreto con los valores de tiempo definidos por el vector n. Los vectores n y x deben tener dimensiones compatibles, es decir deben tener el mismo número de elementos. Así, para este caso para obtener la representación de esta señal en tiempo discreto creamos un vector-tiempo el cual debe tener valores separados por una unidad.

>> u=[zeros(1,10), ones(1,11)];

>> n=-10:10;

>> stem(n,u)

Recuerde que para poder usar las funciones plot y stem, es requisito que los vectores (t y u) ó (n y u) tengan iguales dimensiones. Por esta razón el vector u se forma como una composición de diez ceros y 11 unos, debido a que los arreglos t y n, tienen dimensión 21 dado que incluyen un elemento central el cual es el número cero. Para probar este hecho, se puede hacer uso de la función Matlab llamada size que devuelve como resultado un vector con las dimensiones de la matriz que se le pasa como parámetro.

4) SEÑAL IMPULSO:

La versión discreta de la señal impulso se puede también generar con ayuda de las funciones zeros y ones, realizando una composición como sigue:

>> delta = [ zeros( 1 ,10 ), 1 , zeros( 1 ,10 ) ];

>> n = -10:10;

>> stem(n,delta);

5) SEÑAL RAMPA 

Para generar la señal rampa, tan solo es necesario recordar que esta función puede ser creada, como la composición de una recta Y(x) = x a partir de cero y de la recta Y(x) = 0 para valores de x menores de cero, así la versión discreta se muestra a continuación:

>> t1=0:0.1:10;

>> rampa1=t1;

>> rampa=[zeros(1,101),rampa1];

>> t2=-10:0.1:0;

>> t=[t2,t1];

>> plot(t,rampa)

Ejercicios:

1) Desarrollar un conjunto de comandos Matlab para aproximar las siguientes señales periódicas en tiempo continuo, dibujando 5 ciclos de cada una:

a) Onda Cuadrada, de amplitud 5 Volts, frecuencia fundamental 20 Hz y ciclo útil del 60%.

b) Señal diente de sierra, amplitud 5 Volts y frecuencia fundamental 20Hz

Práctica 2: Señales harmónicas moduladas y su filtrado

Ya pueden consultar su práctica 2, para lo cual les recomiendo que leen sus libros sobre transformada de Fourier.

Por cierto, esta presentación sobre función de transferencia les será útil en los siguientes semanas.

Funciones De Transferencia

View more presentations from Grupo 4 Señales y Sistema

Por cierto, también deben consultar el artículo: How can Trigonometry help to understand Lock-in Amplifier operation? De Ernesto Marín, R. Ivanov. Pues la técnica lock-in también es un método de filtrado de señales, muy aplicado cuando la razón señal/ruido es muy baja.