Algoritmos de simulación y aprendizaje de maquina en MatLab
Los algoritmos de simulación son parte del aprendizaje de maquina y son usados en diversos campos de estudio en los cuales la predicción de determinadas variables permite la toma de decisiones, a continuación les comparto la implementación en Matlab de los algoritmos mas usados.
Algoritmo de gradiente descendente: para la implementación del algoritmo se usó la fórmula de gradiente descendente (1) el cual actualiza los pesos en cada iteración, también se usó la estimación del error (2) para detener el algoritmo cuando su error fuera inferior a 1e-10.
function w = gradiente(X, Y, eta, num)
%validacion de argumentos
if nargin < 3, eta = 0.001; end
%Numero de filas y columnas
[N,D] = size(X);
%extandarizacion
X = zscore(X);
%matriz extendida
Xext = [X,ones(N,1)];
%maximo de iteraciones
Maxiter = num;
%vector de coeficientes
W =ones(1,D+1);
%vector de errores
E =zeros(1,Maxiter);
for iter = 1:Maxiter
%gradiente
W = W - eta*((1/N)*((W*Xext')' - Y)'*Xext);
%error
E(iter + 1) = (1/(2*N))*sum(((W*Xext')' - Y).^2);
%detener a iteracion si el error es minimo
if abs(E(iter + 1) - E(iter)) < 1e-10
break;
end
end
%grafica del error
E = E(1:iter);
plot(E);
%-----------------
% function Y = milogsitic( x )
% Y = 1./(1+exp(x));
% end
Algoritmo de gradiente descendente para clasificacion: Si se desean evaluar un modelo de clasificación se debe implementar la función sigmoide, a continuación una implementacion generando dos conjuntos de muestras artificiales:
function g = sigmoide( x )
g = exp(x)./(1+exp(x));
end
eta = 0.001; %tasa de aprendizaje
X1 = 1.*randn(100,2)+2; %media 2 y varianza 1
X2 = 2.*randn(100,2)+5; %media 5 varianza 2
X = [X1 ; X2];
Y = [ones(100,1);zeros(100,1)];
[N,D] = size(X);
%estandarización --- no es necesaria
%X = zscore(X);
%matriz extendida
Xext = [X,ones(N,1)];
%maximo de iteraciones
Maxiter = 10000;
%vector de coeficientes
W =ones(1,D+1);
%vector de errores
E =zeros(1,Maxiter);
for iter = 1:Maxiter
%version especial con llamado a la funcion sigmoide
W = W - eta*((1/N)*(sigmoide(W*Xext')' - Y)'*Xext);
%error
E(iter + 1) = (1/(N))*sum(abs(((W*Xext')' - Y)));
%detener a iteracion si el error es minimo
if abs(E(iter + 1) - E(iter)) < 1e-10
break;
end
end
plot(X1(:,1),X1(:,2), '*b');
hold on
plot(X2(:,1),X2(:,2), '*r');
%dibujar la frontera
x1 = linspace(-2,10,100);
x2 = -(W(1)/W(2))*x1 - (W(3)/W(2))
plot(x1,x2,'k');
Algoritmo para clasificación basado en funciones discriminantes Gaussianas: Encontrarmos las funciones de densidad de probabilidad (fdp) de las clases y se realiza la clasificación de una muestra con base en la probabilidad de que esa muestra pertenezca a una u otra clase.
%generativo con funcione gausianas
X1 = randn(100,2) + 2;
X2 = 2*randn(100,2) + 5;
%modelo original
% plot(X1(:,1),X1(:,2), '*b');
% hold on
% plot(X2(:,1),X2(:,2), '*r');
X = [X1 ; X2];
Y = [ones(100,1);zeros(100,1)];
%face de entrenamiento
modelo = train(X,Y);
Z = zeros(200,3);
for i=1:200
%funciones de probabilidad para la clase 1
tem = (X(i,:) - modelo(1).media);
p1 = 1/(((2*pi)*sqrtm(det(modelo(1).cova)))) * exp(-0.5 * tem * inv(modelo(1).cova) * tem');
%funciones de probabilidad para la clase 2
tem = (X(i,:) - modelo(2).media);
p2= 1/(((2*pi)*sqrtm(det(modelo(2).cova)))) * exp(-0.5 * tem * inv(modelo(2).cova) * tem');
if (p1 >= p2)
%pertenece a la clase 1
Z(i,:) = [1 X(i,:)];
else
%pertence a la clase 2
Z(i,:) = [0 X(i,:)];
end
end
%clasificaion!!!
plot(Z(Z(:,1)==0,3), Z(Z(:,1)==0,2), '*b');
hold on
plot(Z(Z(:,1)==1,3), Z(Z(:,1)==1,2), '*r');
Algoritmo SMOTE y SubMuestro Inteligente: Este algoritmo nos permite realizar crear nuevas muestras sinteticamente utilizando el algoritmo de K-vecinos, tambien se usa sub muestreo inteligente que permite reducir la clase mayoritaria cuando la base de datos esta desbalanceada
%--------------------------------SMOTE-----------------------------------
%agrego libreria para llamar el algoritmo de submuestreo inteligente
addpath('CSNN');
%cargo los datos que deseo utilizar
datos = load ('glass12.dat');
datos2 = datos;
%ahora seleccionas el numero de muestras minoritarias
datos = datos(datos(:,end) == 1, :); %ya se que son los que tienen 1 (76)
%separamos x y y
% lo valores de entrada x
x = datos(:, 1:end-1);
% los valores de salida en este caso la clase a que pertenece (-1 y 1)
%llamo el algoritmo, le enviamos como parametros, las clases minoritarias
%T, el porcetaje N% que hace referencia a la cantidad de smote que se
%quiere realizar y el numero de vecinos a evaluar K. El algoritmo me
%devolvera nuevas muestras de la clase minoritaria creadas sinteticamente,
%el numero de muestras sera (N/100) * T
syn = smote(x, 200, 5); %152 MUESTRAS NUEVAS
%muestras de la clase minoritaria combinadas con las nuevas muestras;
x = [x;syn];
%armo de nuevo el conjunto de datos
x = [x ones(size(x,1),1)];
datosNew = [x; datos2(datos2(:,end) == -1, :)];
%-------------------------UNDER-SAMPLING----------------------------------
x = datosNew(:, 1:end-1);
% los valores de salida en este caso la clase a que pertenece (-1 y 1)
y = datosNew(:, end);
%para el uso de submuestreo inteligente se uso una libreria ya construida,
%para ello se usa la funcion underSamplig del CSNN, se le entrega como
%parametro, el conjunto de datos x, su salida y, las clases, la relacion de
%salida, y un vector que indica si los datos son numericos o discretos
%(0 para numericos), entrega como salida el conjunto de los nuevos datos (x,y)
[newTrain,newTrainLabel]=UnderSampling(x',y',[1 -1],[1 2],[0,0,0,0,0,0,0,0,0]);
datos = [newTrain' newTrainLabel'] %nuevos datos
%-------------------------IMPLEMENTACION SMOTE--------------------------
function synthetic = smote(T, N, k)
numattrs = size(T,2); %numero de atributos
nt = size(T,1); %cantidad de muestras
if(N < 100)
%combino las filas aleatoreamente
T = T(randperm(size(T,1)),:);
nt = floor((N/100)*nt); %redondeo el numero
N = 100;
end
N = floor(N/100);
%calculo los k vecinos y retorno su indice en la muestra original
cercaindex = kvecinos(k, T);
synthetic = []; %nueva muestra
for i =1:nt
%claculo las nuevas muestras, envio como parametros el % de smote,
%el indice de la muestra, sus k vecinos, el numero de atrivutos y
%las muestras
synthetic = [synthetic; populate(N,i,cercaindex(i,:), numattrs, T)];
end
end
%-------------------------POPULATE----------------------------------
function synthetic = populate(N, i, nnarray, numattrs, sample)
newindex = 1;
synthetic = []; %tamaño de la nueva muestra
%mientras N sea diferente de cero
while (N ~= 0)
%selecciono un indice del array aleatoriamente
nn = randi(size(nnarray,2));
for attr = 1:numattrs
%diferencia entre un atributo seleccionado aleatoriamente de
%los k vecinos mas cercanos y el atributo en cuestion
dif = sample(nnarray(nn),attr) -sample(i,attr);
% numero aleatorio entre 0 -1
gap = rand(1);
%final mente se calcula la nueva muestra
synthetic(newindex, attr) = sample(i,attr) +(gap*dif);
end
%avanzo en el for
newindex = newindex +1;
N = N -1;
end
end
%-------------------------K-VECINOS----------------------------------
function cercaindex = kvecinos(k, x)
num = size(x,1) ; %numero de filas
cercaindex = zeros(num,k);
for i=1:num
%para calcular la distancia se puede hacer manualmente de la
%siguiente manera:
%----------------------------------------------------------------
%N1 = size(x, 1);
%N2 = size(x(i,:), 1);
%L = x(i,:)*x';
%|xi-x*|^2 = sum(Xtrain.^2,2)*ones(1,N1)+ones(N2,1)*sum(X.^2,2)'-2*L
%D = (sum(x(i,:).^2,2)*ones(1,N1)+ones(N2,1)*sum(x.^2,2)'-2*L).^(1/2)
%----------------------------------------------------------------
%sin embargo investigando se encontro que existe una funcion en
%matlab que permite calcular la distancia de manera mas simple
%usando dist:
D = dist(x(i,:), x');
%utilizamos la funcion sort para encontrar los k vecinos mas
%cercanos al vector, luego obtenemos su indice (o) y los extraemos
[c, o] = sort(D);
cercaindex(i,:) = o(:,2:k+1);
end
end
Algortimo K-Vecinos para clasificacion: Se usó el conjunto conocido como base de datos Iris que corresponde a un problema de clasificación de 3 tipos de orquídeas con base en información sobre los pétalos, del repositorio Machine Learning Repository. Estos datos fueron cargados y evaluados usando un modelo para la clasificación de las muestras empleando el método de los k-vecinos con k variando de 1 a 7. Se realizó un cambio en la base de datos cambiando la variable de salida para facilitar su manejo, así: Iris-setosa = 1, Iris-versicolor = 2, Iris-virginica = 3. Se llevaron a cabo mediciones del error de clasificación, así como su respectiva matriz de confusión. A continuación la implementación del algoritmo:
%se cargan los datos, se cambiaron las clases por números Iris-setosa = 1
%Iris-versicolor = 2 Iris-virginica = 3
datos = load('iris.data');
%seleccionamos x y
x = datos(:,1:end-1);
y = datos(:,end);
%separación de datos de entrenamiento y de prueba
%selecciono muestras aleatorias sin repetir 70%
[datosTrain, o] = datasample(datos,floor((size(datos, 1)*70)/100),'Replace',false);
xTrain = datosTrain(:,1:end-1);
yTrain = datosTrain(:,end);
%escojo las muestras que no han sido seleccionadas 30%
c = 1:size(x,1);
i = setdiff(c,o);
xTest = x(i,1:end);
yTest = y(i,end);
%las diferentes clases para construir la matriz de confusión
My = unique(y);
eTem = 0;
kTem = 0;
for k=1:1:7
[e,m]= kvecinosCla(xTrain,yTrain,xTest,yTest,k, My);
if(k == 1)
eTem = e;
mTem = m;
kTem =k;
else
if(e < eTem)
eTem = e;
mTem = m;
kTem =k;
end
end
end
disp(strcat('La matriz con mejor resultado fue en k = ', num2str(kTem)));
disp(mTem);
function [Error,MC] = kvecinosCla(xTrain,yTrain,xTest,yTest,k, My)
%normalización del vector de entrenamiento, se obtiene la media y la
%desviacion
[xTraiN,media,desvia] = zscore(xTrain);
medias = repmat(media,size(xTest,1),1);
xTestN =(xTest- medias)./repmat(desvia,size(xTest,1),1);
num = size(xTestN,1) ; %numero de filas
yResul = zeros(num, 1);
for i=1:num
%distancia
D = dist(xTestN(i,:), xTraiN');
[c, o] = sort(D);
%saco los k vecinos
vecinos = o(:,1:k);
%calculo la media y ese sera la prediccion
yResul(i) = mode(yTrain(vecinos));
end
%error
dif = (yResul(:,end) ~= yTest(:,end));
Error = sum(dif)/num;
disp(strcat('Error de clasificacion= ', sprintf('%0.5e',Error)));
%Matriz de confusión;
disp('Matriz de confucion: ');
[MC, ord] = confusionmat(yTest,yResul);
MC = [ord';MC];
MC = [[0; ord] MC];
disp(MC);
end
Algortimo K-Vecinos con estandarizacion de datos (Regresion):
%se cargan los datos
datos = xlsread('Concrete_Data.xls');
%seleccionamos x y
x = datos(:,1:end-1);
y = datos(:,end);
%separación de datos de entrenamiento y de prueba
%selecciono muestras aleatorias sin repetir 70%
[datosTrain, o] = datasample(datos,floor((size(datos, 1)*70)/100),'Replace',false);
xTrain = datosTrain(:,1:end-1);
yTrain = datosTrain(:,end);
%escojo las muestras que no han sido seleccionadas 30%
i = setdiff(c,o);
xTest = x(i,1:end);
yTest = y(i,end);
for k=1:1:7
kvecinos(xTrain,yTrain,xTest,yTest,k);
end
function kvecinos(xTrain,yTrain,xTest,yTest,k)
%normalización del vector de entrenamiento, se obtiene la media y la
%desviacion
[xTraiN,media,desvia] = zscore(xTrain);
medias = repmat(media,size(xTest,1),1);
xTestN =(xTest- medias)./repmat(desvia,size(xTest,1),1);
num = size(xTestN,1) ; %número de filas
yResul = zeros(num, 1);
for i=1:num
%distancia
D = dist(xTestN(i,:), xTraiN');
[c, o] = sort(D);
%saco los k vecinos
vecinos = o(:,1:k);
%calculo la media y ese será la predicción
yResul(i) = mean(yTrain(vecinos));
end
%error cuadrático medio
dif = (yResul - yTest).^2;
Error = sum(dif)/num;
disp(strcat('Error Cuadratico Medio = ', sprintf('%0.5e',Error)));
mdl = LinearModel.fit(xTest,yResul,'linear');
%coeficiente de determinacion ordinario
disp(strcat('Coe?ciente de Determinacion = ', sprintf('%0.5e',mdl.Rsquared.Ordinary) ));
%grafica de salida real vs predicción del sistema
figure
scatter(yTest, yResul);
title(strcat('Y Real vs Y de prediccion k = ',num2str( k )));
xlabel('Real');
ylabel('Train');
end
Algoritmo de ventana de Parzen con estimador de Nadaraya-Watson: Modelo de predicción basado en el método de ventana de Parzen utilizando el estimador Nadaraya-Watson para diferentes valores de h (0.01<h<10), con medición del error cuadrático medio y el coeficiente de terminación para cada h. se realiza la mejor predicción al tener el menor MSE. A continuación la implementación del algoritmo
%se cargan los datos
datos = xlsread('Concrete_Data.xls');
%seleccionamos x y
x = datos(:,1:end-1);
y = datos(:,end);
%separación de datos de entrenamiento y de prueba
%selecciono muestras aleatorias sin repetir 70%
[datosTrain, o] = datasample(datos,floor((size(datos, 1)*70)/100),'Replace',false);
xTrain = datosTrain(:,1:end-1);
yTrain = datosTrain(:,end);
%escojo las muestras que no han sido seleccionadas 30%
c = 1:size(x,1);
i = setdiff(c,o);
xTest = x(i,1:end);
yTest = y(i,end);
ventana_parzen_Nadara(xTrain,yTrain,xTest,yTest, 0.01);
ventana_parzen_Nadara(xTrain,yTrain,xTest,yTest, 0.05);
ventana_parzen_Nadara(xTrain,yTrain,xTest,yTest, 0.1);
ventana_parzen_Nadara(xTrain,yTrain,xTest,yTest, 1.5);
ventana_parzen_Nadara(xTrain,yTrain,xTest,yTest, 3);
ventana_parzen_Nadara(xTrain,yTrain,xTest,yTest, 6);
ventana_parzen_Nadara(xTrain,yTrain,xTest,yTest, 9);
function ventana_parzen_Nadara(xTrain,yTrain,xTest,yTest, h)
%normalización del vector de entrenamiento, se obtiene la media y la
%desviación
[xTraiN,media,desvia] = zscore(xTrain);
medias = repmat(media,size(xTest,1),1);
xTestN =(xTest- medias)./repmat(desvia,size(xTest,1),1);
num = size(xTestN,1) ; %numero de filas
yResul = zeros(num, 1);
for i=1:num
%distancia
D = dist(xTestN(i,:), xTraiN');
numD = size(D,2);
%kernel
k = exp(-(0.5)*(D/h).^2)/2;
%Nadaraya
yTem = zeros(numD, 1);
tem = D./h;
yTem = (sum(k.*tem*yTrain))./(sum(k.*tem));
yResul(i) = mean(yTem);
end
%error cuadrático medio
dif = (yResul - yTest).^2;
Error = nansum(dif)/num;
disp(strcat('Error Cuadratico Medio = ', sprintf('%0.5e',Error)));
mdl = LinearModel.fit(xTest,yResul,'linear');
%coeficiente de determinación ordinario
disp(strcat('Coe?ciente de Determinacion = ', sprintf('%0.5e',mdl.Rsquared.Ordinary) ));
%grafica de salida real vs predicción del sistema
figure
scatter(yTest, yResul);
title(strcat('Y Real vs Y de prediccion h = ',num2str( h )));
xlabel('Real');
ylabel('Test');
end
Algoritmo de ventana de Parzen para clasificación: Se construyó un modelo para la clasicación de la base de
datos Iris empleando el método de ventana de Parzen evaluando diferentes
valores de h.
%se cargan los datos, se cambiaron las clases por numeros Iris-setosa = 1
%Iris-versicolor = 2 Iris-virginica = 3
datos = load('iris.data');
%seleccionamos x y
x = datos(:,1:end-1);
y = datos(:,end);
%separacion de datos de entrenamiento y de prueba
%selecciono muestras aleatorias sin repetir 70%
[datosTrain, o] = datasample(datos,floor((size(datos, 1)*70)/100),'Replace',false);
xTrain = datosTrain(:,1:end);
%escojo las muestras que no han sido seleccionadas 30%
c = 1:size(x,1);
i = setdiff(c,o);
xTest = x(i,1:end);
yReal = y(i,end);
calcularVetana(xTrain,xTest,yReal,0.01);
calcularVetana(xTrain,xTest,yReal,0.05);
calcularVetana(xTrain,xTest,yReal,0.1);
calcularVetana(xTrain,xTest,yReal,1.5);
calcularVetana(xTrain,xTest,yReal,3);
calcularVetana(xTrain,xTest,yReal,6);
calcularVetana(xTrain,xTest,yReal,9);
function calcularVetana(xTrain,xTest,yReal,h)
num = size(xTest,1);
yTest = zeros(num, 1);
for i=1:num
p1 = ventana_parzen_Cla(xTrain(xTrain(:,end) == 1,1:end-1),xTest(i,1:end), h);
p2 = ventana_parzen_Cla(xTrain(xTrain(:,end) == 2,1:end-1),xTest(i,1:end), h);
p3 = ventana_parzen_Cla(xTrain(xTrain(:,end) == 3,1:end-1),xTest(i,1:end), h);
[c, o] = max([p1 p2 p3]);
yTest(i) = o;
end
dif = (yTest(:,end) ~= yReal(:,end));
Error = sum(dif)/num;
disp(strcat('Error de clasificacion= ', sprintf('%0.5e',Error)));
function suma = ventana_parzen_Cla(xTrain,xTest, h)
num = size(xTest,1) ; %numero de filas
suma = 0;
for i=1:num
%distancia
D = dist(xTest(i,:), xTrain');
k = exp(-(0.5)*(D/h).^2);
suma = sum(k);
end
suma = suma/num;
end
Excelente aporte
ResponderEliminar