29 junio, 2008

Ajustar una linea recta usando algoritmos genéticos en MATLAB

Encuentra los parámetros $m$ y $b$ de la ecuación de la linea recta $y=mx+b$ que mejor se ajusta a los datos $\{(x_i,y_i)\}_{i=1}^{n}$ mediante un sencillo Algoritmo Genético (AG). Se compara el resultado obtenido con AG contra el obtenido con Mínimos Cuadrados.

Contenido

Definición del problema

Se genera un conjunto de datos sintéticos $\{(x_i,y_i)\}_{i=1}^{n}$ donde la variable dependiente $y$ es un vector de números aleatorios en el intervalo de 0 a 1 por lo que se espera que la mejor solución sea cercana a $m=0$, $b=0.5$

x=linspace(1,100); % dominio de 1 a 100
y=rand(1,100);     % números aleatorios en el intervalo de 0 a 1

Parámetros iniciales

Se configua el Algoritmo definiendo parámetros análogos en la Genética

generaciones = 800; % iteraciones del algoritmo
n = 200;            % tamaño de la población de soluciones
probabilidad = 0.6; % probabilidad de que cada gen mute
mutacion = [1 1];   % máximo valor que puede mutar

Población inicial

Se genera aleatoriamente una población inicial de posibles soluciones. La primer columna representa el gen que define la pendiente $m$ de la recta. La segunda columna representa el gen que define la ordenada al origen $b$. Cada uno de los $n$ individuos (renglones) representa una posible solución.

poblacion = rand(n,2);

Iteraciones

La población de soluciones evoluciona mediante la cruza y mutación

s = n;
hw = waitbar(0,'Evolucionando...');
for i = 1:generaciones

Orden aleatorio de parejas

    poblacion(:,3) = rand(s,1);
    poblacion = sortrows(poblacion,+3);

Cruza

    poblacion=[poblacion; poblacion(1:2:end-1,1) poblacion(2:2:end,2) zeros(fix(s/2),1); poblacion(2:2:end,1) poblacion(1:2:end-1,2) zeros(fix(s/2),1)];

Mutacion

    s = size(poblacion,1);
    poblacion=[poblacion; poblacion+[(rand(s,3)<probabilidad).*(rand(s,3)*2-1).*repmat([mutacion 0],s,1)]];

Cálculo de la aptitud (o viabilidad)

La tercer columna de poblacion representa la aptitud (viabilidad) de la solución, menor es mejor.

    s = size(poblacion,1);
    m = repmat(poblacion(:,1),1,100);
    b = repmat(poblacion(:,2),1,100);
    X = repmat(x,s,1);
    Y = repmat(y,s,1);
    Y2=m.*X+b;
    poblacion(:,3)=sum(abs(Y-Y2),2);

Seleccion

    poblacion = unique(poblacion,'rows');
    poblacion = sortrows(poblacion,+3);
    s = min([size(poblacion,1),n]);
    poblacion = poblacion(1:s,:);
    waitbar(i/generaciones,hw)
end
close(hw)

Resultados

Se compara el resultado obtenido con AG contra el obtenido con mínimos cuadrados

error_ag = poblacion(1,3);
p = polyfit(x,y,1);
y3 = polyval(p,x);
error_mc = sum(abs(y3-y));
plot(x,y,'.',x,polyval(poblacion(1,1:2),x),x,y3,':')
axis([0 100 -0.5 1.5])
legend('Datos Sintéticos','Algoritmos Genéticos (AG)','Mínimos Cuadrados (MC)')
text(1,-.167,['AG: Error = ' num2str(error_ag) ', m = ' num2str(poblacion(1,1)) ', b = ' num2str(poblacion(1,2))])
text(1,-.333,['MC: Error = ' num2str(error_mc) ', m = ' num2str(p(1)) ', b = ' num2str(p(2))])