%% Componentes principales con introducción explícita de regresores no lineales: ejemplo Matlab.
% © 2020, Antonio Sala Piqueras, Universitat Politècnica de València.
% Presentación en vídeo en http://personales.upv.es/asala/YT/V/pcanlx.html
% Este código ejecutó correctamente en Matlab R2018b.
% Tabla de Contenidos
% 1.-Motivación y esbozo de teoría
% 2.-Ejemplo
%     PCA lineal
%     PCA no lineal
%         Reconstrucción de la curva (variedad) en coordenadas originales
% 3.- Conclusiones
%%
%% 1.-Motivación y esbozo de teoría
% Si se forma un vector de regresores (características, "feature vector") que pueda incluir no-linealidades, \phi(x), el PCA puede descubrir que los datos se agrupan en variedades ("manifolds") con curvatura, descritas por una combinación lineal de las características no lineales.
% El PCA lineal sólo puede determinar subespacios (hiperplanos) donde se agrupen los datos.
% Por lo tanto, dados (x1, x2, ... xn), en vez de formar X=[x1 x2 .... xn]  para PCA lineal, formamos X=[\phi(x_1) \ \phi(x_2)\ \dots\ \phi(x_n)] como matriz de datos. La descomposición en valores singulares de Xproveerá de información sobre si algunas combinaciones lineales de estas características tienen una desviación típica muy pequeña, que podamos, por tanto, suponer que son un "modelo no lineal" identificado, generalizando la idea tras el PCA o Total Least Squares lineal.
% El PCA no-lineal puede considerar unas características de más dimension que los datos originales para modelar mejor: \phi(x)=(x, \|x\|, sin(x_1) log(x_1-x_2) \dots).
%
% Nota: aunque los "ruidos de medida" en los elementos de x podrían considerarse no correlados, en las características no lineales aparecerán repetidas veces (introduciendo correlación) y una hipotética "distribución normal" del ruido de medida sería deformada en \phi(x). Por tanto, la interpretación estadística (TLS) ya no es formalmente válida, pero sigue siendo útil en la práctica.
%% 2.-Ejemplo

N=250;
r=1.5+randn(N,1)*0.075;
th=rand(N,1)*2*pi;
x=r.*cos(th);
y=r.*sin(th);
plot(x,y,'o'), grid on, title('Datos Originales'), axis equal

%% PCA lineal

X0=[x y]
medias=mean(X0)
desvtip=std(X0)
X=(X0-medias)*inv(diag(desvtip));

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

% No detecta dimensionalidad reducida, todos los valores son similares.
%% PCA no lineal
% Sospechamos que podría existir alguna relación cuadrática entre los datos. Definimos un vector de características:

phi=@(x1,x2) [x1 x2 x1.^2 x2.^2 x1.*x2];

% y formamos la matriz de datos extendida.

X0=phi(x,y)
medias=mean(X0)
desvtip=std(X0)
X=(X0-medias)*inv(diag(desvtip)) %escalado media cero, desv.típ. unidad.

% X contiene el vector de características normalizado para que todos sus componentes tengan media cero y varianza unidad.
% \phi_{normalizado}(x)=(\phi(x)-mean(\phi(x)))./std(\phi(x))

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

% Hemos detectado que hay un componente (el último) de mucha menor desviación típica que el resto. Por tanto, los datos podrían estar próximos a una variedad de dimensión 1 (una curva en el plano) expresable como "[polinomio de grado 2] = 0".
%
%% Reconstrucción de la curva (variedad) en coordenadas originales
% La ecuación de la curva, en concreto, sería \phi_{normalizada}(x)*V(:,5)=0.
% Para escribirlo en las coordenadas de los datos, hay que deshacer la normalización:
% 0=\frac{\phi(x)-media}{desvtip}*V(:,5)=\phi(x)*\frac{V(:,5)}{desvtip}-media*\frac{V(:,5)}{desvtip}=\phi(x)*w-m

w=inv(diag(desvtip))*V(:,5)
m=medias*w
phitrans=@(x,y) phi(x,y)*w-m
plot(x,y,'x') %datos
hold on
fimplicit(@(x,y) phi(x,y)*w-m)  %modelo estimado
fimplicit(@(x,y) x.^2+y.^2-1.5^2) %somos DIOS, y sabemos de verdad qué originó los datos.
hold off, grid on, axis equal
legend('datos','PCA polin.','exacto')

%% 3.- Conclusiones
% Si se conoce el tipo de no-linealidad subyacente en un conjunto de datos, puede incorporarse al PCA para obtener mejores resultados en identificación y posterior detección de cambios.