Octópodo: junio 2008

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=m*X+b que mejor se ajusta a los datos (x,y) mediante un sencillo Algoritmo Genético (AG). Se compara el resultado obtenido con AG contra el obtenido con Mínimos Cuadrados.

Definición del problema

Se genera un conjunto de datos sintéticos x, y 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))])