Caminando hacia ninguna parte

Qué encontrarás en esta entrada?

• Algoritmo para crear en MATLAB una caminata al azar en dimensión menor o igual a 3.
• Gráficas de las trayectorias.
• Un poco de teoría matemático-estadística de fondo. 

No es la primera vez que oigo en clase la expresión "caminata al azar" (o "random walk"). Y no es una referencia a lo perdido que estemos en el mundo y cómo damos bandazos de un lado a otro sin encontrar nuestro sitio :p, sino a un problema que consiste en que una partícula puede moverse en intervalos discretos de tiempos, dando un salto a cada intervalo, hacia ciertas direcciones ponderadas con sendas probabilidades.

Este, a priori, curioso ingenio teórico resulta estar en correspondencia con ciertas situaciones físicas reales, como es el caso del "movimiento browniano", del cual ya hablamos en esta página.

Gráfica realizada por camaz en 3D

El caso es que se presta bastante a ser programado, así que anoche, en un rato, a falta de algo mejor que hacer, hice un sencillo script en MATLAB que lo ilustra. El script básicamente dibuja la trayectoria de una partícula que se rija por la dinámica anteriormente descrita en tres dimensiones. Las probabilidades por defecto son iguales en todas las direcciones principales del espacio (x⁺, x⁻, y⁺, y⁻, z⁺ y z⁻), pero se pueden cambiar introduciéndolas en los argumentos del script, por lo que es trivial pasar a un caso bidimensional (P(z⁺)=P(z⁻)=0), o unidimensional (P(y⁺)=P(y⁻)=P(z⁺)=P(z⁻)=0).

Caminata al azar en 1D con 10 saltos (equiprobable izquierda/derecha)

Caminata al azar en 2D con 10.000 saltos (equiprobable izquierda/derecha/arriba/abajo). Es curioso lo que me recuerdan estas gráficas a mapas, yo en este dibujo veo África, Asia y Europa, incluyendo al Reino Unido xD!

Teoría matemática para 1D (equiprobable)

En una dimensión, en cada salto, tenemos la opción de ir hacia la izquierda o la derecha. La probabilidad de dar n⁺ saltos a la derecha viene dada entonces por la distribución binomial:

La posición de la partícula será la diferencia entre el número de saltos multiplicada por la escala. La posición más probable está relacionada con la media del número de saltos hacia uno de los dos lados:

Implementando esto en el programa, podemos comparar lo que pasa empíricamente con la estimación teórica, y los resultados no son del todo malos.

Cuando la probabilidad hacia ambos lados es idéntica, el valor esperado es cero (si se puede mover de igual manera hacia la izquierda que hacia la derecha, no hay forma de elegir ninguna de las dos direcciones). En esos casos se demuestra que el recorrido libre medio es del orden de la raíz cuadrada de n, lo cuál también está implementado en el script.

Código

El código que he utilizado es el sigueinte:

function []=camaz(N,p1,p2,p3,p4,p5,p6)
%
% (CAM)inata al (AZ)ar
%
% Fecha: 09/03/11
% Autor: Astaroth (O.R.G.)
% Web: http://astarothsworld.blogspot.com/
%
% Este script dibuja una caminata al azar en tres (o menos) dimensiones y realiza
%unas simples estimaciones estadísticas.
%
% Modo de uso:
% ------------
%
% camaz(N,p(x+),p(x-),p(y+),p(y-),p(z+),p(z-))
%
% Donde N es el número de saltos y p(x_i+-) es la probabilidad de que la partícula
%salte en la dirección positiva (o negativa) de la coordenada x_i.
%
% Tenga en cuenta que la suma de las probabilidades ha de ser 1:
%
% p(x+)+p(x-)+p(y+)+p(y-)+p(z+)+p(z-)=1
%
% La última probabilidad es innecesaria, puesto que se calcula usando la expresión
%anterior.
%
% Si no introduce ninguna probabilidad, se representará por defecto el caso 3D
%equiprobable.
%
if ( nargin == 1 )

p1=1/6;

p2=1/6;

p3=1/6;

p4=1/6;

p5=1/6;

p6=1/6;

end
if ( p1+p2+p3+p4+p5 > 1 )

disp('Las probabilidades que ha introducido son incorrectas!')

disp('No se puede continuar con el programa.')

else

if ( p1+p2+p3+p4+p5+p6 > 1 )
disp('Las probabilidades que ha introducido son incorrectas!')
disp('Se puede continuar recalculando la probabilidad para el eje z negativo.')
disp('Nueva probabilidad p(z⁻):')
p6=1-(p1+p2+p3+p4+p5);
disp(p6)
end

%Inicialización de valores:

close all

x1=0;

y1=0;

z1=0;

plot3(x1,y1,z1,'g*')

grid

hold on

for i = 1:N
x0=x1;
y0=y1;
z0=z1;
r=rand;
if ( r < p1 )
x1=x0+1;
else if ( r < p1+p2 )
x1=x0-1;
else if ( r < p1+p2+p3 )
y1=y0+1;
else if ( r < p1+p2+p3+p4 )
y1=y0-1;
else if ( r < p1+p2+p3+p4+p5 )
z1=z0+1;
else
z1=z0-1;
end
end
end
end
end v=[x1,y1,z1]-[x0,y0,z0]; plot3([x0,x1],[y0,y1],[z0,z1])
end plot3(x1,y1,z1,'r*') d=(x1^2+y1^2+z1^2)^(1/2); xlabel('Eje x') ylabel('Eje y') zlabel('Eje z') title('Caminata al azar') disp('Distancia recorrida (en valor absoluto):') disp(d) if ( p3+p4+p5+p6 == 0 )
disp (' ') if ( p1 ~= p2 )
disp('Estimación teórica de la posición basada en la distribución binomial para 1D:')
xmed=N*(2*p1-1);
disp(xmed)
else
disp('Estimación teórica de la distancia recorrida basada en la distribución binomial para 1D:')
xeff=sqrt(N);
disp(xeff)
end
end

disp(' ') disp('Distancia media por salto:') disp(d/N)

end

Como siempre, podéis encontrar este programa en la sección de descargas de la página, o pinchar directamente aquí para descargarlo.