%% Mejor predicción lineal. Ejemplo 2 variables.
% © 2018, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo en personales.upv.es/asala/YT/V/prelim.html
%
% Este código ejecutó sin errores en Matlab R2018b.
% Tabla de Contenidos
% Generación y visualización de los datos
% Análisis de varianzas-covarianzas
% Predicción de y dado x
% Predicción de x dado y
% Gráfica de resultados
%
%% Generación y visualización de los datos

Nmuestras=600;
[x,y,ysecreto]=formulasecreta(Nmuestras);

xmedio=mean(x); ymedio=mean(y);
x=x-xmedio;y=y-ymedio; %normalizamos a media 0
ysecreto=ysecreto-ymedio;%este es "desconocido"
plot(x,y,'o'), grid on

%% Análisis de varianzas-covarianzas

varx=cov(x)
stdx=sqrt(varx)
vary=cov(y)
stdy=sqrt(vary)
varxy=cov([x y])
varx=varxy(1,1);
vary=varxy(2,2);
covxy=varxy(1,2);
correlacion=covxy/sqrt(varx*vary)

%% Predicción de y dado x

bestpredygivenx_coef=covxy/varx
%coincide con mínimos cuadrados:
pinv(x)*y
vary_posteriori=vary-covxy*inv(varx)*covxy

% Comparemos la desviación típica sin información (a priori) y la desviación típica a posteriori

stdypriori=sqrt(vary)
stdyposteriori=sqrt(vary_posteriori)

porcent_NO_explicado=stdyposteriori/stdypriori
sqrt(1-correlacion^2)

%% Predicción de x dado y

bestpredxgiveny_coef=covxy/vary
%coincide con mínimos cuadrados:
pinv(y)*x
varx_posteriori=varx-covxy*inv(vary)*covxy

% Comparemos la desviación típica sin información (a priori) y la desviación típica a posteriori

stdxpriori=sqrt(varx)
stdxposteriori=sqrt(varx_posteriori)

porcent_NO_explicadox=stdxposteriori/stdxpriori
sqrt(1-correlacion^2)

%% Gráfica de resultados

plot(x,y,'o')
hold on
plot(x,bestpredygivenx_coef*x,'dy')
plot(x,bestpredygivenx_coef*x+2*stdyposteriori,'.y')
plot(x,bestpredygivenx_coef*x-2*stdyposteriori,'.y')
plot(bestpredxgiveny_coef*y,y,'xm')
plot(bestpredxgiveny_coef*y+2*stdxposteriori,y,'.m')
plot(bestpredxgiveny_coef*y-2*stdxposteriori,y,'.m')
plot(x,ysecreto,'.k')
hold off
grid on
axis tight
legend('Data','y dado x','','','x dado y','','','secreto','Location','best')

%

function [xc,y,ysecreto]=formulasecreta(Nmuestras)
    xc=max(-7,min(7,randn(Nmuestras,1)*2+(randn(Nmuestras,1).^(3))*.5));
    n=(randn(Nmuestras,1)).^2*2.5+randn(Nmuestras,1)*3.75;
    ysecreto=3*xc+sin(xc).*(1+1.2*xc)+2.5;
    y=(ysecreto-2.5)+n;
end