%% Simulación e Interpolación (Kriging) con procesos gaussianos no estacionarios
% © 2020, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/wienkrig.html
% Este código ejecutó correctamente con Matlab R2020b
% Objetivo: comprobar cómo, aparte de cambiar la función de covarianza, la estimación en procesos no estacionarios basados en el kernel de covarianza es idéntica al caso estacionario.
% Comprobar las relaciones de la covarianza del proceso de Wiener (ruido tras un integrador) o el proceso de Wiener integrado (ruido tras doble integrador) con otras técnicas de interpolación.
% Table of Contents
% Modelo estocástico "a priori"
% Estimación en puntos intermedios: ridge regression (kriging)
% Conclusiones
% Apéndice: funciones auxiliares
%%
%% Modelo estocástico "a priori"
% Supondremos que tenemos una función que pasa por cero "seguro", y luego pasa por otros puntos, habiendo algunas funciones más probables que otras según la medida inducida por un Kernel de covarianza.
% Nota: podríamos poner ciertas condiciones de "función media no nula" o de "incertidumbre en t=0", omitido por brevedad.

Ts=2e-2;
Xtest=0:Ts:7.5;Ntest=length(Xtest);
ktt=KernelMatr(Xtest,Xtest);
Q=chol(ktt+1e-12*eye(Ntest))'; %para hacer realizaciones ktt=Q*Q'

%
% Hago realizaciones

proc=Q*randn(Ntest,1);
plot(Xtest,proc), grid on
v=diff(proc)/Ts;
plot(v), title("Velocidad"), grid on
plot(diff(v)/Ts,'.'), title("Aceleración"), grid on

%

priorvar=diag(ktt); %varianza a priori
plot(Xtest,sqrt(priorvar)), grid on, title("desviación típica a priori del proceso gaussiano")
surf(Xtest,Xtest,ktt), axis tight, title("Covarianza entre muestras"), shading interp

%% Estimación en puntos intermedios: ridge regression (kriging)
% Supongamos que tenemos estas muestras medidas (datos de entrenamiento):

X=[0.4 1 2 3 4 5.5]';N=length(X);
Y=[0.3 -1.5 -1.5 3 4 5.5]';

% Calculamos el Kernel de los datos de entrenamiento medidos

K=KernelMatr(X,X)+.01^2*eye(N); %suponemos datos corrompidos por ruido desv. típ. 0.01 de medida

% Luego, la covarianza entre puntos a estimar y puntos de entrenamiento

Kt=KernelMatr(Xtest,X);

% La fórmula de la mejor predicción lineal obtiene la media

pred_mean=Kt*(K\Y);

% Y la matriz de varianzas-covarianzas del error de predicción (si los datos han sido generados por el modelo propuesto):

predvar=ktt-Kt*(K\Kt');
preddt=sqrt(diag(predvar)); %desviación típica
plot(Xtest,pred_mean,'LineWidth',3), hold on
plot(X,Y,'xk','MarkerSize',12,'LineWidth',3)
plot(Xtest,[pred_mean-2*preddt pred_mean+2*preddt],'-.','Color',[.7 .8 .5],'LineWidth',2)
grid on, legend("estimado","medido","interv. conf.",'Location','best')
hold off, axis tight
plot(diff(pred_mean)/Ts), grid on, title("Velocidad")
plot(diff(diff(pred_mean))/Ts^2), grid on, title("Aceleración")
plot(diff(diff(diff(pred_mean)))/Ts^3), grid on, title("Jerk")

%% Conclusiones
% Que el proceso subyacente sea o no estacionario no cambia en nada las fórmulas de ridge-regresión/interpolación/Kriging asociadas.
% El Kernel originado por el proceso de Wiener mantiene constantes al extrapolar, hace interpolación lineal entre muestras.
% El Kernel originado por la integral del proceso de Wiener mantiene rectas al extrapolar, y tiene relación con las splines cúbicas de interpolación dado que hace un ajuste cúbico a trozos (y''' constante \rightarrowy'' recta\rightarrowy' cuadrática\rightarrowy es un polinomio cúbico).
% En el caso uni-dimensional (funciones del "tiempo"), la versión en "espacio de estados" de esta regresión es el Filtro de Kalman (extrapolación causal "a la derecha de la última muestra") o el suavizado Rauch-Tung-Striebel (equivalente a los desarrollos aquí).
%
%% Apéndice: funciones auxiliares

function K=KernelMatr(X1,X2)
    n1=length(X1);    n2=length(X2);
    K=zeros(n1,n2);
    for i=1:n1
        for j=1:n2
            K(i,j)=kappa3(X1(i),X2(j)); 
        end
    end
end


%stationary, exponential kernel
function kk=kappa1a(x1,x2)
    kk=exp(-abs(x2-x1));
end

%stationary, Gaussian kernel
function kk=kappa1b(x1,x2)
    kk=exp(-abs(x2-x1)^2/2);
end

%Wiener Process, preserves constants, NON stationary
%cond. inciales W(0)=0 en media y varianza.
%entrada ruido PSD=1
function kk=kappa2(x1,x2)
    kk=min(x1,x2);  %válido si x1>=0, x2>=0
end

%esto es "Kalman" (RTS) del doble integrador
%Integrated Wiener Process, preserves lines, NON stationary
%cond. inciales W(0)=0, W'(0)=0 en media y varianza.
%entrada ruido PSD=1
function kk=kappa3(x1,x2)
    s=min(x1,x2); t=max(x1,x2);
    kk=s^2/2*(t-s/3);%válido si x1>=0, x2>=0
end