%% Identificación estática 1 salida, 2 entradas
% © 2021, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Una presentación/screencast de los contenidos de este archivo aparece en: http://personales.upv.es/asala/videos/id2.html
%
% Vamos a identificar un modelo y=Cu, con y vector 1x1, C matriz 1x2, u vector 2x1.

C=[2 -1]; %este vector es "desconocido", queremos "adivinarlo"

% Conocemos, u, y... varias muestras de ellos.

entradasexperim1=[ 1 -1 4 -5 6 7 9 -3 -8; ...
           2  4 -3 3 1 0 5 2 3]'
       

% *Cada muestra está en una FILA de la matrix entradasexperim1...

info1=svd(entradasexperim1)

entradas2=[1 1.1 1.05 0.1 -1;...
         2 2.1 2.15  0.2 -1.9]';
info2=svd(entradas2)
[U,S,V]=svd(entradas2)

% Para el parámetro (en comp. principal, esto es, una cierta combinación de parámetros [c1 c2]) más difícil de identificar, tenemos que

min(info1)/min(info2)

% Este número es cómo de "más informativo" es el experimento "entradasexperim1" que "entradas2". 150 veces!

%entradas=entradasexperim1; para otras pruebas
entradas=randn(40,2);
svd(entradas)
Nmuestras=size(entradas,1)
Nvariables=size(entradas,2)

% Generemos las salidas con una cierta matriz C, definida arriba.
% Generamos salida "ruidosa": como en el modelo original y=Cu cada dato es columna, pero está en las matrices en "fila", entonces vamos a identificar C' en y'=u'*C'.

%rng(3241) si queremos comparar con la misma secuencia de núms. aleatorios.
salidalimpia=entradas*C';% inaccesible
size(salidalimpia)
salidas= salidalimpia+randn(size(salidalimpia))*.25;
yy=salidas;
X=[entradas];
svd(X)
%X=[entradas entradas.^2 entradas.^3];
%svd(X)

th=pinv(X)*yy;
Cestim=th'

%
% Comprobemos el ajuste:

error=yy-X*th

% Calculamos la varianza y desviación típica del error:

    Varerror=(error'*error)/(Nmuestras-2); %2 params por columna...
    std_error=sqrt(Varerror);
std_error

% Ese std_error será una estimación del ruido de medida en la variable de salida... (si es verdad que ha sido generado con un ruido de medida aleatorio, claro... si no, incluye ruido sistemático --bias--).
% Como el ruido de medida se multiplica por pinv(X), se multiplicará por la inversa de la información para producir parámetros estimados:

pi_svd=svd(pinv(X))

% La matriz de varianzas covarianzas de [c1 c2] es:

Var_th=Varerror(1)*inv(X'*X) %VC param estimados

% Lo cual producirá, en cada uno de los errores una posible corrección si tuviéramos más datos, de varianza, para la medida "k", dada por:
% \nu_k=varianza(X_k\theta)= E(X_k\theta\theta^TX_k^T)=X_kE(\theta\theta^T)X_k^T
% y desv. típica dada por la raíz cuadrada. Por ejemplo, para la fila 7 tendríamos:

correccion_error7=sqrt(X(7,:)*Var_th*X(7,:)');
[error(7) correccion_error7*2]

% Calculando la media de los cuadrados de esas correcciones de error para todos los errores:

desviacion_media_por_vzparametros=sqrt(trace(X*Var_th*X'/Nmuestras)) 
%divido por Nmuestras porque es calculo teórico de vza

% Eso tiene que ser la desviacion típica "en media" que podría añadirse, por la incertidumbre en \theta, a la salida del modelo, al tomar más datos.
%
% La desviación típica de cada parámetro es la raíz cuadrada de la diagonal de la matriz de varianzas-covarianzas:

VMat=[diag(Var_th)];
std_C=sqrt(VMat)'

%% Intervalos de confianza parámetros estimados
% Tenemos un 99% de confianza que si el modelo ha sido originado por lineal+"ruido de medida dist normal", el proceso "real" esté entre:

fprintf('Valores mínimos 99%% [%g %g]',Cestim-3*std_C)
fprintf('Valores máximos 99%% [%g %g]',Cestim+3*std_C)

%% Conclusiones
% El experimento identifica un estimado de la matriz C

Cestim

% Con una desviación típica en cada coeficiente dada por:

std_C

% y una desviación típica del error de predicción en la variable de salida dada por:

std_error

% Que podría reducirse con más datos hasta, aproximadamente, como mucho:

std_no_se_explicaconmasdatos=max(0,std_error-desviacion_media_por_vzparametros)

%