%% PCR, Principal Component Regression, regresión por componentes principales, ejemplo Matlab
% © 2018 Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Vídeo-presentación en http://personales.upv.es/asala/YT/V/pcrm.html
% Este código ejecutó correctamente en Matlab R2018b.
% Tabla de Contenidos
% Preliminares
%     Generación de datos
%     Análisis preliminar, matrices VC
%     Mínimos cuadrados estandard
% PCR
%     Cálculo del predictor y dev. típica de sus parámetros
%     Estudio del error de predicción según número de componentes PCR elegidos
%% Preliminares
%% Generación de datos

N=600;
[X,Y]=formulasecreta(N);%ya sale media cero.
size(X)
maxcomponentes=size(X,2)
size(Y)

% Generemos datos de validación:

Nval=N;
[Xval,Yval]=formulasecreta(Nval); %ya sale media cero.

%% Análisis preliminar, matrices VC

SigmaX=X'*X/(N-1);
desvTipX=sqrt(diag(SigmaX))'
SigmaY=Y'*Y/(N-1);
desvTipY=sqrt(diag(SigmaY))'
VarianzaTotalY=trace(SigmaY)
Cov_XY=X'*Y/(N-1)

%% Mínimos cuadrados estandard

PredictorLS=SigmaX\Cov_XY
ErrorPred=Y-X*PredictorLS;
size(ErrorPred)
VarErrorPredLS=ErrorPred'*ErrorPred/(N-1);
VErrdiag=diag(VarErrorPredLS);
sqrt(VErrdiag)' %desv. típica del error en cada variable
VarTotErr=trace(VarErrorPredLS) %variación total del error
norm(ErrorPred,'fro')^2/(N-1) %interpretación determinista, coincide.
PorcentajeVzaResidual=VarTotErr/VarianzaTotalY
PorcentajeVzaExplicada=1-PorcentajeVzaResidual

%
%% PCR
%% Cálculo del predictor y dev. típica de sus parámetros
% calculemos el PCA de entrada X

[U,S,V]=svd(X,'econ');

% Transformemos a componentes principales:

Xpca=X*V; %U*S
size(Xpca)

% Comprobamos que la varianza es diagonal:

SigmaXpca=Xpca'*Xpca/(N-1)
desvtipXPCA=sqrt(diag(SigmaXpca))' 
diag(S/sqrt(N-1))' %coincide
Cov_XpcaY=Xpca'*Y/(N-1);
PredictorPCR=SigmaXpca\Cov_XpcaY
PredictorPCR_UnidadesORIG=V*PredictorPCR 

% Todo OK, coincide con PredictorLS, porque si gastamos todos los componentes es un mero cambio de variable.

PredictorLS

% Calculemos ahora la desviación típica de parámetros estimados del Predictor PCR...

inv_excitacion=diag(inv(Xpca'*Xpca)); %es matriz diagonal...
inv_excitacion'
VarParEstimPCR=kron(inv_excitacion,VErrdiag');
DesvTipParametrosEstimadosPCR=sqrt(VarParEstimPCR)

% Para hacernos una idea de la incertidumbre "relativa", presentemos el predictor dividido por su desv. típica; esto indica las veces la desv. típica necesarias para llegar a cero.

PCRrel=PredictorPCR./DesvTipParametrosEstimadosPCR

%% Estudio del error de predicción según número de componentes PCR elegidos
% Error de predicción según componentes elegidos con datos de entrenamiento:

for ncomp=1:maxcomponentes
    ErrorPred=Y-X*V(:,1:ncomp)*PredictorPCR(1:ncomp,:);
    VarErrorPred(ncomp)=trace(ErrorPred'*ErrorPred/(N-1));
end
VarErrorPred/VarianzaTotalY %fracción de varianza residual

% Lo mismo con datos validación:

VarianzaTotalYval=norm(Yval,'fro');
for ncomp=1:maxcomponentes
    ErrorPredV=Yval-Xval*V(:,1:ncomp)*PredictorPCR(1:ncomp,:);
    VarErrorPredV(ncomp)=trace(ErrorPredV'*ErrorPredV/(Nval-1));
end
VarErrorPredV/VarianzaTotalYval
bar(0:maxcomponentes,[VarianzaTotalY VarErrorPred; VarianzaTotalYval VarErrorPredV]'/diag([VarianzaTotalY,VarianzaTotalYval]))
grid on
legend('TrainingData','ValData'), title('Varianza residual')
xlabel('Num. componentes PCR')
bar(0:maxcomponentes,sqrt([VarianzaTotalY VarErrorPred; VarianzaTotalYval VarErrorPredV]'/diag([VarianzaTotalY,VarianzaTotalYval])))
grid on
xlabel('Num. componentes PCR'),legend('TrainingData','ValData'), title('Desv. Típica residual')

% Nota: Verificar sensibilidad resultados a ruido (repetir varias veces).
% Recordar que otros métodos (PLS) buscan maximizar covarianza y no varianza de información X y pueden ser mejores (más reducción de varianza error con mismo número de componentes).
%

function [X,Y]=formulasecreta(N)
ncomps=3;
    MatrizDeformaX=[5 -2 0;-5 1 0;0 0 5];
    X=randn(N,ncomps)*MatrizDeformaX;    
    ParamsIdeales=[2 5;1 -2;1 1]/2;
    ruido=randn(N,2)*6.5;   
    Y=ruido+X*ParamsIdeales;
    %y añadimos a X una columna que no tiene NADA QUE VER COn Y:
    X=[X randn(N,1)];
    X=X-mean(X);
    Y=Y-mean(Y);
end