%% Script de la Practica8
function Practica8_Prediccion

%% Situate en el directorio donde se encuentra corpus.mat

%%% limpiamos el entorno de matlab
clear
close

%% cargamos los datos desde corpus.mat
load corpus.mat

%% preparamos un conjunto de entrenamiento de 2 clases
%% XT2clases: variables predictivas del corpus de entrenamiento con 2 clases
%% YT2clases: etiequeta de salida del corpus de entrenamiento
XT1=Xpitraining(Ytraining==1,:);
XT2=Xpitraining(Ytraining==2,:);
YT1=Ytraining(Ytraining==1);
YT2=-1*ones(size(Ytraining(Ytraining==2)));
XT2clases=[XT1;XT2];
YT2clases=[YT1;YT2];

%% preparamos un conjunto de test de 2 clases
%% Xt2clases: variables predictivas del corpus de test con 2 clases
%% Yt2clases: etiequeta de salida del corpus de test
Xt1=Xpitest(Ytest==1,:);
Xt2=Xpitest(Ytest==2,:);
Yt1=Ytest(Ytest==1);
Yt2=-1*ones(size(Ytest(Ytest==2)));
Xt2clases=[Xt1;Xt2];
Yt2clases=[Yt1;Yt2];

%%
%% ejercicio 1: discriminante lineal, estimacion por minimos cuadrados
%%

XT=XT2clases;
Xt=Xt2clases;
YT=YT2clases;
Yt=Yt2clases;

%% ampliamos la matriz X para el termino independiente
XT= [XT,ones(size(XT,1),1)];
Xt= [Xt,ones(size(Xt,1),1)];

%% estimacion
%% W = ((XT'*XT)^-1)*XT'*YT
W = (inv(XT'*XT))*XT'*YT

%% regla de clasificacion basada en la funcion discriminante lineal
gt = Xt*W %valor de la funcion discriminante para cada caso de test
yt = gt>0 %Si gt(i) > 0, entonces yt(i) 1, sino -1

yt(yt==0)=-1

a = yt == Yt

%% calculo del error de test
err = 1 - (sum(a) / size(a,1)) % casos incorrectamente clasificados dividido numero de casos de test

%%
%% ejercicio 2: clasificador gausiano
%%

%% clasificacion binaria LDA y QDA

XT=XT2clases;
YT=YT2clases;
Xt=Xt2clases;
Yt=Yt2clases;

%% estimacion y clasificacion lineal de las clases 1 y 2

[class,err,POSTERIOR,logp,coeff] = classify( Xt, XT, YT, 'linear');

%% estudia la variable POSTERIOR
POSTERIOR
%% estudia la variable coeff para averiguar la frontera de decision
coeff
%% calcula el error de test
err

%% estimación y clasificación cuadrática de las clases 1 y 2
[class,errT,POSTERIOR,logp,coeff] = classify( Xt, XT, YT, 'quadratic');

%% estimacion y clasificacion lineal de las clases 1, 2 y 3
class
%% estimacion del error de test
TODO
%% Cálculo de la matriz de confusión
confus(Ytest,class)

%% Clasificacion de 3 clases basada en la extraccion de caracteristicas
%% mediante PCA

%% tipificar la union de conjuntos de entrenamiento y test
Xtip=tipificar([Xpitraining;Xpitest]);
XTtip = Xtip(1:size(Xpitraining,1),:);
Xttip = Xtip(size(Xpitraining,1)+1:size(Xpitraining,1)+size(Xpitest,1),:);


%% estudiamos kolmogorov-smirnov, miramos k-s de 1 var en mitad de la señal
[h,p,k,c] = kstest(Xtip(:,floor(size(Xtip,2)/2)),[],0.05,0);

%% Proyeccion PCA del conjunto de entrenamiento tipificado
TODO = princomp(XTtip);
%% Proyeccion del conjunto de test tipificado mediante loadings
TODO
%% clasificacion de 3 clases con los 10 primeros scores
[class,err,POSTERIOR,logp,coeff] = classify( TODO



%% ejercicio 3: sequentialfs: forward y backward-selection, 5 repeticiones

for dir=['fb']
for i=1:5
    c = cvpartition(Ytraining,'k',10);
    opts = statset('display','iter','TolFun',1e-3);
    fun = @(XT,yT,Xt,yt)...
          (sum(yt ~= classify(Xt,XT,yT,'linear')));
    [fs,history] = sequentialfs(fun,Xpitraining,Ytraining,'cv',c,'options',opts, 'direction',dir);
    Xpitrainingfs = Xpitraining(:,fs);
    Xpitestfs = Xpitest(:,fs);
    [class,err,POSTERIOR,logp,coeff] = classify(Xpitestfs, Xpitrainingfs, Ytraining, 'linear');
    confus(Ytest,class);
    'sw lda 3c'
    err=sum(Ytest~=class)/size(Ytest,1)

end
end


%% ejercicio 4: svm


XT=XT2clases;
YT=YT2clases;
Xt=Xt2clases;
Yt=Yt2clases;

%clasificacion con kernel polinomico de grado 3
C=2;
svmStruct = svmtrain( TODO %X,t,kernel_,width,polypower,C
classt = svmclassify( TODO %svmStruct,Xt)
errT=sum(YT~=classt)/size(YT,1);
class = svmclassify(svmStruct, Xt);
'svm 2c'
err= TODO %calcula el error
confus(Yt+2,class+2);


%PCA[1:2] + SVM

%con PCA de vector completo claificamos 2 clases,
%en PC tenemos los loadings ordenados de mayor a menor

%tipificar
Xtip=tipificar([XT;Xt]);

XTtip = Xtip(1:size(XT,1),:);
Xttip = Xtip(size(XT,1)+1:size(XT,1)+size(Xt,1),:);

%miramos k-s de 1 var en mitad de la señal
[h,p,k,c] = kstest(Xtip(:,floor(size(Xtip,2)/2)),[],0.05,0);

[PC,scoresT,variance] = princomp(XTtip);
scorest=Xttip*PC;

svmStruct = svmtrain( TODO
classT = svmclassify( TODO
errT=sum(YT~=classT)/size(YT,1);
class = svmclassify(svmStruct, scorest(:,1:2));
'pca svm 2c'
err=sum(Yt~=class)/size(Yt,1)
confus(Yt+2,class+2);

%plot de vectores soporte, fronteras y casos de test
figure
scoresT12=scoresT(:,1:2);
scorest12=scorest(:,1:2);
Step=0.5;
mn = min(scoresT12);
mx = max(scoresT12);
[x1,x2]=meshgrid(floor(mn(1)):Step:ceil(mx(1)),floor(mn(2)):Step:ceil(mx(2)));
[n11,n12]=size(x1);
[n21,n22]=size(x2);
XG=[reshape(x1,n11*n12,1) reshape(x2,n21*n22,1)];
alpha_index=svmStruct{3};
KG = kernel_func(scoresT(alpha_index,1:2),XG,'poly',1,3);
alpha=svmStruct{1};
w_0 = svmStruct{2};
f = (YT(alpha_index).*alpha)'*KG + w_0;
plot(scoresT12(alpha_index,1),scoresT12(alpha_index,2),'go');
hold
plot(scorest12(1:52,1),scorest12(1:52,2),'.')
plot(scorest12(53:75,1),scorest12(53:75,2),'.r')
contour(x1,x2,reshape(f,[n11,n12]),[0,0]);
hold off


end