\documentclass[a4paper,12pt,twoside]{article}
%\documentclass[10pt,journal,a4paper]{IEEEtran}
\usepackage[left=2cm,right=2cm,top=2.5cm]{geometry}%bottom=1.7cm

%\usepackage[utf8]{inputenc}
\setlength{\parindent}{0pt}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc} % codifica windows
%\usepackage[utf8]{inputenc}    % codifica linux
%\usepackage{textcomp}
%\usepackage[applemac]{inputenc} 
%\usepackage[italian]{babel}
\usepackage[dvips]{graphicx}
\usepackage{float}
\usepackage{epsf}
\usepackage{epsfig}
\usepackage{verbatim}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsbsy}
\usepackage{fancyhdr}
\pagestyle{fancy}
\usepackage{hyperref}

\usepackage[final]{listings}
\usepackage[dvips]{graphicx}
\usepackage[usenames,dvipsnames]{xcolor}
%\usepackage[numbered,framed]{pcks/mcode}

\hypersetup{colorlinks,%
citecolor=blue,%
filecolor=black,%
linkcolor=black,%
urlcolor=blue,%
pdftex}

\numberwithin{equation}{section}
\numberwithin{figure}{section}

\renewcommand{\sectionmark}[1]{\markright{\thesection.\ #1}}
\fancyhf{}
\fancyhead[LE,RO]{\thepage}
\fancyhead[LO,RE]{\scshape\rightmark}
\renewcommand{\headrulewidth}{0.4pt}
\renewcommand{\footrulewidth}{0pt}
%\oddsidemargin = 1.3 cm
%\evensidemargin = 0.5 cm
\title{Parte relativa ai Tmote}
\author{Riccardo}
\date{14/02/2010}
\begin{document}

\section{Raccolta dati}

\subsection{I dispositivi utilizzati}
Per eseguire le misurazioni delle grandezze di interesse sono stati utilizzati i dispositivi  \textit{MoteIV Tmote Sky} equipaggiati con sensori per il rilevamento di temperatura, umidità e luminosità. Sebbene questi sensori siano predisposti per una comunicazione dei dati mediante rete wireless si è scelto di effettuare l'acquisizione off-line sfruttando la memoria flash per immagazzinare i dati raccolti più delle informazioni aggiuntive che sono stati letti alla fine mediante usb. La scelta di utilizzare questa modalità di lettura è dipesa principalmente dall'esigenza di un attento utilizzo del consumo energetico, infatti sebbene questi dispositivi siano progettati per garantire un basso consumo ed abbiano un processore in grado di funzionare con tensioni che arrivano fino a 1.9 V i componenti necessari all'esperimento (i sensori e la memoria flash) necessitano di una tensione costante di circa 3 V fornita da due batterie di tipo AA alcaline. Se la tensione in ingresso diminuisse eccessivamente potrebbero accadere dei malfunzionamenti al convertitore ADC con una perdita di bit significativi (da 14 bit a 12 bit upper) e alla memoria flash che per una scrittura sicura richiede una almeno 2.6 V. Si fa notare che la scelta di una alimentazione mediante porta usb per questa esperienza non è possibile poichè in questo caso l'alimentazione a 5 V provoca un surriscaldamento del processore posto in prossimità del sensore di temperatura, sfalsando così la misura di temperatura. In questo modo si è riusciti ad avere una alimentazione corretta durante tutta la raccolta dei dati, in Figura ??? è riportato l'andamento della tensione di ogni dispositivo raccolta proprio per monitorare il buon funzionamento di ogni dispositivo.\\
Viene riportata qui di seguito una breve descrizione dei sensori utilizzati dai \textit{Tmote} e successivamente del software necessario:
\begin{enumerate}
\item \textit{Umidità e temperatura}:  Il sensore di temperatura e umidità utilizzato sul mote è il sensore Sensiron SHT11. Esso consta di un sensore di temperatura il cui principio di funzionamento si basa sulla dipendenza dalla temperatura della banda di non conduzione di un semiconduttore ed ha range di funzionamento: -40/+123.8 C. La sua risoluzione è a 14 bit (può essere ridotta a 12 bit per applicazioni che richiedono velocità elevate o un utilizzo con potenza in ingresso molto bassa). La calibrazione è gestita internamente in base a delle tabelle, la formula di conversione dalla parola a 14 bit $T_{raw}$ presente nella memoria flash a gradi centigradi è
\begin{equation}
T = b_1+b_2T_{raw}-b_3(T_{raw}-b_4)^2 \quad [\text{\textdegree}C]
\end{equation}
con costanti $b_1=-39.6$, $b_2=0.01$, $b_3=-2 \cdot 10^{-8}$ e $b_4=7000$.
 \\ Il sensore per rilevare l'umidità è basato sulla variazione di capacità di un polimero in funzione dell'umidità presente nell'ambiente circostante e ha range di funzionamento 0/100 \% RH. La misura viene registrata in una parola a 12 bit $H_{raw}$ che può scendere ad 8 bit se necessario. La formula per la conversione è:
\begin{equation}
RH_{true} = (T-25)(c_1+c_2H_{raw})+c_3+c_4H_{raw}+c_5H_{raw}^2 [\%]
\label{rh_true}
\end{equation}
con constanti $c_1=0.01$, $c_2=8\cdot 10^{-6}$, $c_3=-4$,$c_4=0.0405$ e $c_5=-2.8\cdot 10^{-6}$. Per maggiori informazioni si rimanda al datasheet ???.
\item \textit{Luminosità} I sensori presenti nel mote per misurare la quantità di luce sono due fotodiodi prodotti dalla Hamamatsu. In particolare è presente il modello S1087 sensibile alla sola luce visibile (320 - 730nm) e il modello S1087-01 che rivela anche lo spettro infrarosso (320 - 1100nm). Le misure sono memorizzate in parole da 16 bit unsigned e la conversione in lux é data da:
\begin{equation}
AR = 625AR_{raw} \quad [lux]
\label{ar}
\end{equation}
con $AR_{raw}$ il valore in bit fornito dal sensore S1087. Analogamente a partire dal dato in bit $TR_{raw}$ del sensore S1087-01 si implementa
\begin{equation}
TR = 76.9TR_{raw} \quad{[lux]}\,.
\label{tr}
\end{equation}
Per la fase di identificazione la misura che viene utilizzata è quella fornita dal secondo sensore, ovvero quella comprensiva della radiazione visibile e dell'infrarosso poichè entrambe influiscono sul comportamento dinamico dell'edificio.
\item \textit{Memoria flash} Il mote è equipaggiato con una unità di memoria flash esterna di tipo seriale di capienza 1024kb. Essendo tuttavia configurata in lettura e scrittura solo nei primi quindici settori lo spazio effettivamente disponibile è di 983040 byte. La struttura dati presente consta di 6 campi da 16 bit ciascuno, il che significa una capacità totale di 81920 campioni. Questo valore è in realtà ridotto a 65535 campioni poichè parte dello spazio è utilizzato per controllare se dei dati non vengono scritti in memoria (di fatto l'ultimo dei sei campi è occupato da un numero che aumenta di un passo ad ogni istante di campionamento). Contando che il tempo di campionamento adatto per le dinamiche termiche di un edificio è dell'ordine dei 5-10 minuti si capisce che la memoria disponibile nel mote è più che sufficiente per lo svolgimento dell'esperienza.\\
\end{enumerate}
Si fornisce qui di seguito una breve descrizione del codice che viene utilizzato per la gestione dei mote. La sessione di acquisizione viene suddivisa in tre fasi dal firmaware caricato su ogni mote.  La prima parte è dedicata al setup che viene effettuato comunicando con il mote tramite usb e pc. In questa fase vengono gestite le operazioni principali di lettura della memoria flash e di cancellazione della stessa. La seconda fase è dedicata alla sincronizzazione e il passaggio avviene allo scadere di un timer opportunamente impostato. In questa fase le comunicazioni usb sono interrotte e gli unici eventi che il mote è può rilevare sono ricezione del pacchetto radio di sincronizzazione o la pressione del tasto user (che causa la partenza di un nuovo timer e l'invio di un nuovo pacchetto di sincronizzazione). Con la terza fase si inizia l'acquisizione vera e propria dei dati. Il modulo radio viene disattivato ed un timer posto uguale al tempo di campionamento comanda la lettura dei dati rilevati dai sensori e la scrittura su flash. Il principio di funzionamento è riportato in Figura \ref{strutturaStati}.\\
L'interfaccia con il PC è sviluppata in Java su piattaforma Linux ed utilizza le librerie di comunicazione seriale e le classi per la gestione dei pacchetti fornite con TinyOS 2.x. La lettura del mote fornisce in outuput un file di testo contenente i valori registati dai sensori che vengono convertiti in unità ingegneristiche da uno script in Matlab secondo (\ref{rh_true})-(\ref{tr}). La programmazione dei mote è avvenuta eseguendo lo script $EggMote.h$ in cui si precedentemente impostato il tempo di campionamento.

\begin{figure}[H]
\begin{center}
\includegraphics[width=7cm]{figure/pp.eps}
\caption{Struttura degli stati e delle funzionalità del mote}
\label{strutturaStati}
\end{center}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Metodi di identificazione ai sottospazi}
I metodi di identificazione a sottospazi sono particolarmente adatti per l'identificazione di sistemi MISO (o MIMO) di grandi dimensioni per la loro efficienza computazionale. Vengono implementati utilizzando strumenti derivati dall'algebra lineare (per effettuare opportune proiezioni ortogonali e oblique) come la decomposizione QR e SVD di cui esistono implementazioni molto efficienti. Si noti che in questo modo si evita l'utilizzo di minimizzazioni ricorsive e quindi si evitano i problemi di convergenza che esse comportano (e che possono capitare con i metodi PEM).\\
Il problema è dunque, date le osservazioni $\{y_1,\ldots,y_T,u_1,\ldots,u_T\}$, $y_t\in \mathbb{R}^p$, $u_t\in \mathbb{R}^m$ costruire il modello:
\begin{equation}
\left\{ \begin{array}{rcl}
x(t+1) & =& Ax(t)+Bu(t)+Ke(t)\\
y(t) &=& Cx(t)+e(t)
\end{array}\right.
\end{equation}
che descriva i dati (con $e(t)$ rumore bianco).
%Si richiede che il sistema abbia tutti i modi eccitabili dall' ingresso $u$ (stocastico o deterministico) ovvero si richiede la controllabilità della coppia $(A,[BQ^{1/2}])$. Inoltre tutti i modi devono essere osservabili dall'uscita $y$, cioè la coppia $(A,C)$ deve essere osservabile. Il problema è dunque quello di trovare un ordine $n$ per il modello e le matrici $A,B,C,D,Q,R,S$ dati gli ingressi $u_0,\ldots,u_{s-1}$ e le uscite $y_0,\ldots,y_{s-1}$. A partire da una proiezione attraverso una decomposizione ai valori singolari si ottiene una stima dell'ordine e la matrice di osservabilità (e se necessario la sequenza degli stati del sistema). Nel caso di un ingresso deterministico le matrici $A,B,C,D$ del sistema vengono trovate risolvendo in maniera esatta di un opportuno sistema lineare. Nel caso di ingresso stocastico il sistema lineare viene risolto tramite i minimi quadrati e successivamente si ricavano dai residui le matrici di covarianza $Q,S,R$ e di conseguenza gli ingressi stocastici del sistema. In questa trattazione verranno usati alcuni strumenti la cui trattazione completa esula lo scopo prefissato. Per approfondimenti si rimanda a ??. 

A partire dal processo stostico $\mathbf{y}(t)$ si scrivono i dati nella forma $\{y_t,\ldots,y_{t+N-1}\}=\mathbf{Y}_t $ per creare una corrispondenza biunivoca (al crescere del numero $N$ di dati) tra lo spazio astratto di variabili aleatorie $span\{y(t),\ldots,y(s) \}$ e lo spazio generato dalle righe $\mathbf{Y}_t$ della matrice di Henkel diagonale a blocchi:
\begin{equation}
Y_{[t,s]}=
\begin{bmatrix}
y_t & y_{t+1} & \ldots & y_{t+N-1}\\
y_{t+1} & y_{t+2} & \ldots & y_{t+N} \\
\vdots & \vdots & \ddots & \vdots \\
y_s & y_{s+1} & \ldots & y_{s+N-1}
\end{bmatrix}
\end{equation}
di dimensione $Y_{[t,s]}\in \mathbb{R}^{p(s-t)xN}$ e da questa si definiscono le matrici delle uscite passate e future
\begin{equation}
Y_{[0,2t-1]}=
\begin{bmatrix}
Y_p \\
Y_f
\end{bmatrix}
=
\begin{bmatrix}
Y_{[0,t-1]} \\
Y_{[t,2t-1]}
\end{bmatrix}
\end{equation}
e si svolgono delle operazioni analoghe con i dati d'ingresso $\mathbf{u}(t)$. A questo punto si può costruire una base per lo spazio di stato $X_t$ attraverso una proiezione adatta delle sequenze dei dati e successivamente da questa base si risolve nel senso dei minimi quadrati il sistema:
\begin{equation}
\left\{ \begin{array}{rcl}
\hat{X}_{t+1} & \simeq & A\hat{X}_t+B\hat{U}_t+KE_t \\
Y_t & \simeq & C\hat{X}_t+E_t
\end{array}\right.
\end{equation}
dove $\hat{E}_t:=Y_t-\hat{C}X_t$. 
\\Si definisce poi la matrice di osservabilità del sistema come:
\begin{equation}
\Gamma:=
\begin{bmatrix}
C\\
CA\\
CA^2\\
\vdots \\
CA^f
\end{bmatrix}
\end{equation}
e matrici triangolari inferiori a blocchi di Toeplitz con i parametri di Markov del sistema:
\begin{equation}
H_d=
\begin{bmatrix}
0 & 0 & \ldots & 0 \\
CB & 0 & \ldots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
CA^{f-1}B & CA^{f-2}B & \ldots & 0
\end{bmatrix} \quad,
\end{equation}
\begin{equation}
H_s=
\begin{bmatrix}
I & 0 & \ldots & 0 \\
CK & I & \ldots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
CA^{f-1} & CA^{f-2}K & \ldots & I
\end{bmatrix} \quad .
\end{equation}
 \\ Utilizzando ora la relazione
\begin{equation}
\mathbf{Y}_{[p,p+f]}=\Gamma \mathbf{X}_p+H_d\mathbf{U}_{[p,p+f]}+H_s\mathbf{E}_{[p,p+f]}
\end{equation}
si può effettuare la proiezione obliqua 
\footnote{
La proiezione obliqua \'e la scomposizione in due sottospazi tra loro non ortogonali. Si denota come $\hat{E}_{||B}[\cdot|C_t]$, ossia la proiezione lungo lo spazio generato delle righe di $B$ sullo spazio generato dalle righe di $C_t$. La formula estesa \'e:
\begin{equation}
 \hat{E}_{\|B_t}[A_t|C_t] = \hat{\Sigma}_{ac|b}\hat{\Sigma}_{cc|b}^{-1}C_t
\label{obl}
\end{equation}
Dove sotto le ipotesi di ergodicit\'a le matrici di covarianza campionaria (condizionata e non) sono: $ \hat{\Sigma}_{ab} := \frac{1}{N}\sum_{i=0}^{N-1}a_{t+i}b_{t+i}^{\top} $ e $ \hat{\Sigma}_{ab|c} := \hat{\Sigma}_{ab} - \hat{\Sigma}_{ac}\hat{\Sigma}_{cc}^{-1}\hat{\Sigma}_{cb} $
}
 (mediante decomposizione QR)
\begin{equation}
\Gamma \mathbf{X}_p \simeq \hat{\mathbf{Y}}_{[p,p+f]}:=E_{||\mathbf{U}_{[p,p+f]}}\left[\mathbf{Y}_{[p,p+f]}|\mathbf{Z}_{[0,p)}\right].
\label{proiez_obl}
\end{equation}
%inserire figura

\begin{figure}[H]
\begin{center}
\includegraphics[width=7cm]{figure/proiezione_obliqua.eps}
\caption{Rappresentazione geometrica della proiezione obliqua effettutata in (\ref{proiez_obl})}
\label{proiezione_obliqua}
\end{center}
\end{figure}

Per stimare $\Gamma$  ( e se necessario $\mathbf{X}_p$ ) si utilizza una decomposizione "`SVD pesata"'
\begin{equation}
U_nS_nV_n^T \simeq W_r \hat{\mathbf{Y}}_{[p,p+f]}W_c
\end{equation}
e si ottiene
\begin{equation}
\hat{\Gamma}:=W_r^{-1}U_nS_n^{1/2} \qquad \hat{\mathbf{X}}_p:=\hat{\Gamma}_r^{-L}\hat{\mathbf{Y}}_{[p,p+f]}
\end{equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Stato dell'arte}
L'identificazione black-box è ampiamente utilizzata per ricostruire la dinamica di modelli da un insieme finito di dati ingresso-uscita. In questo scenario una situazione particolare che si verifica in svariati ambiti (ad esempio in ingegneria chimica, in econometria, nella visione computazionale, e così via) è quella in cui si richiede l'identificazione di sistemi di grandi dimensioni in cui entrano in gioco un gran numero di variabili. In un contesto simile un punto cardine delle procedure di identificazione dovrebbe essere quello di favorire la sparsità dei risultati, ovvero di estrarre dal gran numero di sottosistemi presenti nel sistema principale solo il sottoinsieme che influenza l'uscita del sistema in maniera più significativa. Il fatto che ogni variabile possa dipendere solo da un sottoinsieme di altre variabili si esprime nel caso Gaussiano statico con delle condizioni di indipendenza condizionata fra sottoinsiemi di variabili. Nel caso dinamico, cioè quando i dati osservati sono traiettorie di processi stocastici (possibilmente stazionari) le condizioni di indipendenza condizionale si esprimono nel fatto che la predizione dell'andamento di una variabile (che verrà chiamata uscita) richieda solo il passato di qualche altra variabile (che verranno chiamate ingressi) più eventualmente il proprio. Si può rappresentare ciò con un grafo in cui i nodi sono le variabili, gli archi diretti sono le funzioni di trasferimento non nulle e i cappi (archi che si chiudono sullo stesso nodo da cui partono) codificano la dipendenza della variabile dal proprio passato, come si può vedere in Figura ??. In generale sia il sistema dinamico sia la struttura che rappresenta il legame fra le variabili dovranno essere estratti dai dati. In questo contesto di identificazione gli usuali metodi di minimizzazione dell'errore di predizione (PEM) possono diventare computazionalmente problematici poichè la teoria sottostante ai metodi di stima dell'ordine del modello (ad esempio AICC e BIC) inizia a vacillare. Questo problema insieme al problema legato alla presenza di non-convessità è stato risolto con l'introduzione degli algoritmi ai sottospazi che invece di fornire la stima utilizzando metodi iterativi di ottimizzazione non lineare (come in PEM), si basano su degli strumenti efficienti e robusti derivati dall'algebra lineare come la decomposizione QR e SVD. Tuttavia anche per questi metodi la presenza di un gran numero di parametri da stimare rispetto al numero di dati disponibili resta un problema. In questa situazione dunque è possibile che sia gli algoritmi parametrici tradizionali che quelli ai sottospazi possano fornire risultati non soddisfacenti.\\
Oggigiorno sono stati introdotti dei nuovi approcci parametrici (come LAR e LASSO) che inducono sparsità nel risultato sfruttando la norma $l_1$. Un siffatto termine peso porta con sè la cosiddetta proprietà di bi-separazione, cioè la preferenza di soluzioni con molti termini nulli a spesa di pochi termini ma di modulo maggiore. Sono inoltre stati ampliati con delle versioni "a gruppi" (GLars e GLasso) in cui interi gruppi di termini sono posti a zeri ottenendo così una sparsità in termini di funzioni di trasferimento piuttosto che di singoli termini. Tuttavia la maggior parte del lavoro in questo ambito è stato svolto in uno scenario statico a discapito del caso dinamico. \\
Gli algoritmi che verranno testati maggiormente in questo lavoro adottano un punto di vista Bayesiano e prendono come punto di partenza un nuovo paradigma sviluppato in [??] che si basa sulla stima non parametrica delle risposte impulsive. In questo modo invece di dover postulare una struttura finito dimensionale per le funzioni di trasferimento del sistema (con ordine quindi da stimare) si cerca la risposta impulsiva all'interno di uno spazio infinito dimensionale. Il mal condizionamento intrinseco a un problema così formulato viene risolto utilizzando dei metodi di regolarizzazione Bayesiana. Più nel dettaglio si può dire che le risposte impulsive sono modellate come dei processi Gaussiani la cui autocovarianza viene chiamata $stable\, spline\, kernel$ ed include il vincolo di BIBO stabilità del risultato.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Preliminari: nuclei per l'identificazione}


\subsection{Definizione del problema}
\label{sec:DefinizioneDelProblema}

Sia $\{z_t\}_{t\in \mathbb{Z}},\, z_t\in \mathbb{R}^m $ un processo stocastico, stazionario che descrive l'evoluzione temporale di qualche variabile di interesse. Con un piccolo abuso di notazione in questo contesto $z_t$ indica sia la variabile aleatoria che la sua realizzazione. Come mostrato in figure \ref{fig:rete} ogni componente del processo $\{z_t\}_{t\in \mathbb{Z}}$ può essere pensato come un nodo di una grafo. L'obiettivo di questi algoritmi di identificazione è di trovare modelli lineari dinamici che descrivono tutti i componenti di $\{z_t\}$ in funzione degli altri componenti. A questo scopo si definisce $y_t:=z_t^{[i]}$ (la componente i-esima di $z_t$) come l'uscita e tutte le altre componenti $u_t:=z_t^{[-i]}\in \mathbb{R}^{m-1}$ come gli ingressi. Ovviamente il procedimento può essere ripetuto per ogni $i=1,\ldots,m$ ottenendo in questo modo una descrizione completa di tutte le variabili di $z_t$ in funzione delle rimanenti. In questa sezione però per questioni di chiarezza si sceglie un valore specifico per $i$, per esempio $i=1$. In questo caso $z_t$ può essere rappresentato come
\begin{equation*}
z_t = \begin{bmatrix} y_t \\ u_t \end{bmatrix}
\end{equation*}

Si definiscono gli insiemi delle variabili passate all'istante $t$
$$ Y^t=[y_{t-1},y_{t-2},\ldots], \qquad U^t=[u_{t-1},u_{t-2},\ldots] $$

Dalla stazionarietà del processo $\{z_t\}_{t\in \mathbb{Z}}$ segue che anche i processi $\{y_t\}_{t\in \mathbb{Z}}$ e $\{u_t\}_{t\in \mathbb{Z}}$ sono processi stocastici, congiuntamente stazionari che possono quindi essere pensati rispettivamente come ingressi e uscita di un sistema dinamico tempo-invariante. Possono essere quindi descritti dalla seguente relazione:
\begin{equation}
y_t = \sum_{k=1}^{\infty}f_k u_{t-k}+\sum_{k=0}^{\infty}g_k e_{t-k}
\label{eq:SisDin}
\end{equation}
dove $f_k\in\mathbb{R}^{1\times m}$ (avendo assunto assenza di retroazione diretta, cioè $f_0=0$) e $g_k\in\mathbb{R}$ sono coefficienti incogniti della risposta impulsiva e $e_t$ rappresenta un processo di innovazione. Esso è equivalente all'errore di predizione di un passo
\begin{align}
	e_t &:= y_t - \hat{y}_{t|t-1} = y_t - E[y_t|Y^t,U^t] \\
	E[y_t|Y^t,U^t] &:= \sum_{j=1}^{m-1}\left[ \sum_{k=1}^{\infty}h_k^{j}u_{t-k}^{[j]} \right] +\sum_{k=1}^{\infty}h_k^{[m]}y_{t-k}.
\label{predittore_gay}
\end{align}
Le sequenze $h_k:=[h_k^{[1]},h_k^{[2]},\ldots,h_k^{[m]}]\in\mathbb{R}^{1\times m}, k\in \mathbb{Z}^+$ rappresentano i coefficienti della risposta impulsiva del predittore e servono a descrivere sistemi (BIBO) stabili e si ha quindi $h^{[m]}\in l_1(\mathbb{Z}^+)$.\\
Nell'ambito della PEM la predizione di sistemi dinamici può essere vista come predizione dei coefficienti della risposta impulsiva $h_k$ del predittore. Si vuole prestare particolare attenzione a situazioni in cui non tutte le variabili sono necessarie per predire $y_t$, il che in termini matematici equivale ad avere $h_k^{[i]}=0$ per ogni $k\in\mathbb{Z}^+$. In una rappresentazione grafica c'è una arco diretto $u_k^{[i]}$ al nodo $y_k$ se e solo se esiste $k\in\mathbb{Z}^+$ tale che $h_k^{[i]}\neq0$ per $i=1,\ldots,m-1$; inoltre non c'è un capo se e solo se esiste $k\in\mathbb{Z}^+$ tale che $h_k^{[m]}\neq0$. Per esempio il grafo rappresentato in Figura \ref{fig:rete} ha $h_k^{[5]}=h_k^{[1]}=0, \, \forall k\in\mathbb{Z}+$, mentre $h_k^{[2]}=h_k^{[3]}=h_k^{[4]}$ e $h_k^{[6]}$ non sono identicamente nulle. Questo significa che per predire $y_t$ serve il passato di $u^{[2]},\ u^{[3]},\ u^{[4]}$ e di $y$ stesso.\\
\begin{figure}[t]
    \begin{center}
      \includegraphics[width=0.7\textwidth]{figure/rete}
    \caption[rete]
    {\label{fig:rete}
    Grafo che rappresenta l'interconnessione di $m=6$ variabili. Gli archi solidi rappresentano le connessioni nel modello dinamico con il nodo $y_t:=z_t^{[1]}$ da tutti gli altri nodi. L'assenza degli archi da $u_t^{[i]},\ i=1,5$ a $y_t$ indica che $h_k^{[i]}=0,\ i=1,5,\ \forall k \in \mathbb{Z}^+$. Il nodo $y_t$ ha l'arco entrante indicando la dipendenza dall'errore di predizione a un passo $e_t$. Gli archi punteggiati indicano dipendenze tra gli ingressi $u_t$.}
    \end{center}
\end{figure}
In pratica non si sa a priori quali variabili saranno significative per predire l'uscita $y_t$. I metodi PEM e N4SID non effettuano alcuna selezione di ingressi e stimano un modello ``completo'' che tiene conto di tutte le misure. Questo è un motivo per cui questi due algoritmi a volte producono stime scarse.\\
Qui di seguito si presentano degli algoritmi che favoriscono la sparsità e che sono quindi capaci di individuare una struttura del sistema dinamico simile a quella riportata in Figura \ref{fig:rete} ed allo stesso momento riescono a stimare le risposte impulsive $h_k^{[i]}$ nulle.





\subsection{La regolarizzazione}
Un approccio ampiamente utilizzato per ricostruire una funzione da un insieme di misure $\{y_t\}$ consiste nel minimizzare un funzionale di regolarizzazione appartenente ad un RKHS (Reproducing Kernel Hilbert Space) $\mathcal{H}$ associato ad un nucleo $K$ simmetrico e definito positivo.\\
Date $N$ misure, la regolarizzazione ai minimi quadrati in $\mathcal{H}$ stima la funzione sconosciuta come
\begin{equation}
\hat{h}=arg \min_{h}\sum_{t=1}^{N}(y_t-\Gamma_t[h])^2+\eta \|h\|^2_{\mathcal{H}} 
\label{kernel based regularization}
\end{equation}
dove $\{\Gamma_t\}$ sono dei funzionali lineari e limitati in $\mathcal{H}$ legati al modello di misura mentre il parametro $\eta$ è uno scalare positivo che permette di settare il trade-off fra l'errore di predizione e la regolarità della soluzione. La funzione che minimizza (\ref{kernel based regularization}) può essere interpretata da un punto di vista Bayesiano come la stima a minima varianza di $h$ quando $h$ è visto come un processo Gaussiano a media nulla con autocovarianza $K$ e $\{y_t-\Gamma_t[h] \}$ un rumore bianco indipendente da $h$. Spesso la conoscenza a priori è limitata al fatto che i segnali in gioco siano continui e a energia limitata. In questo caso la $h$ viene modellata con un integrale $p$-dimensionale di un rumore bianco di intensità unitaria. L'autocorrelazione $W_p$ di $h$ si può così scrivere come:
\begin{equation}
W_p(s,t)=\int_0^1G_p(s,u)G_p(t,u)du\,
\label{stable spline kernel}
\end{equation}
\begin{equation}
G_p(r,u)=\frac{(r-u)_+^{p-1}}{(p-1)!},\quad (u)_+=
\left\{ \begin{array}{l}
u\quad \text{if} \, u\geq 0 \\
0\quad \text{if} \, u < 0
\end{array}\right.
\label{id spline kernel}
\end{equation}
dove $G_p$ è l'autocovarianza associata all'interpretazione Bayesiana di una spline di regolarizzazione di ordine $p$. Se $p=2$ si ottiene un nucleo a spline cubica.\\
Il problema principale nell'utilizzo di un nucleo come (\ref{id spline kernel}) è che non tiene conto della stabilità della risposta impulsiva. Infatti se l'autocovarianza di $h$ è proporzionale a $W_p$ (come nel caso cubico) la varianza di $h(t)$ è nulla nell'istante iniziale e poi diverge all'infinito. Per questo motivo è stato definito un nucleo apposta per l'identificazione di sistemi lineari in cui $h$ rappresenti la risposta impulsiva di un sistema stabile (con covarianza non nulla nell'instante iniziale che converge esponenzialmente a zero). Questo viene garantito prendendo la covarianza proporzionale alla classe di nuclei
\begin{equation}
K_p(s,t) = W_p(e^{-\beta s},e^{-\beta t})\, \quad s,t\in \mathbb{R}^+
\label{classe di nuclei stabili}
\end{equation}
dove $\beta$ è uno scalare positivo che regola il decadimento della varianza; essendo in pratica un parametro sconosciuto sarà trattato come un iperparametro da stimare a partire dai dati. Tenendo presente le formule precedenti si vede che con $p=2$ l'autocovarianza diventa ciò che in letteratura viene chiamato $Stable\,Spline\,kernel$:
\begin{equation}
K_2(t,\tau)=\frac{e^{-\beta(t+\tau)}e^{-\beta \max(t,\tau)} }{2}-\frac{e^{-3\beta \max(t,\tau)}}{6}
\label{stable spline kernel}
\end{equation}
E' stato dimostrato (vedi[??]) che se $h$ è un processo Gaussiano a media nulla con autocovarianza $K_2$ allora con probabilità uno le realizzazioni di $h$ sono risposte impulsive di sistemi BIBO stabili.\\
E' chiaro che lavorando in pratica con segnali a tempo discreto si parlerà di versioni campionate $h_k$ di $h(t)$ così che $\{h_k\}_{k\in \mathbb{Z}^+}$ verranno viste come realizzazioni campionate del nucleo $K(i,j)$.\\

Si considera ora il problema di stimare il parametro $\theta \in \mathbb{R}^m$ dal modello lineare
\begin{equation}
Y=X\theta +W
\label{arx modello}
\end{equation}
dove si ha $Y\in \mathbb{R}^N$ è il vettore dei dati di uscita, $X\in \mathbb{R}^{N\times m}$ è il vettore di regressione e $W\in \mathbb{R}^N$ è un rumore bianco vettoriale a media nulla con $\mathbb{E}[WW^T]=\sigma^2I$. Ciò che verrà fatto è considerare ogni componente di $\theta$ appartenente ad uno spazio infinito dimensionale RKHS con nucleo $K(s,t)$ come in (\ref{classe di nuclei stabili}). In altre parole $\theta_i$ viene visto come un processo gaussiano con covarianza $K(s,t)$. La matrice di regressione $X$ sarà un operatore lineare contente la storia passata di $u$ e $y$.\\

Il problema di selezionare le variabili più significative può essere affrontato come un problema di ricerca sparsità nella soluzione con alcune importanti differenze. In particolare si tratta ora di un problema a gruppi, poichè per escludere dalla procedura di identificazione delle variabili è necessario che vengano poste a zero delle funzioni di trasferimento. Ogni gruppo di parametri descrive dunque una risposta impulsiva $h^{[i]}$. Per evitare il problema relativo alla selezione del modello e della parametrizzazione non lineare si preferisce non utilizzare algoritmi a gruppi parametrici scegliendo di lavorare in un contesto non parametrico. I parametri vivono dunque in un RKHS che come già detto possono essere trattati mediante probabilità a priori e regolarizzazione.\\
I due approcci che verranno utilizzati sono:
\begin{enumerate}
\item SS-GLAR: Versione a gruppi del problema posto dal Lasso
\footnote{Nel Lasso la scelta del regressore avviene risolvendo il problema nella forma:
\begin{equation}
\hat{\theta}:= arg \min_{\theta}\|Y-X\theta\|_2^2+\gamma_1\|\theta\|_1
\end{equation}
con $\gamma_1$ parametro di regolarizzazione compreso fra 0 e 1.
}
portato in un contesto non-parametrico dove i gruppi $h^{[i]}$ appartengono ad uno spazio infinito dimensionale. Per includere il vincolo del problema è necessario risolvere un problema di regolarizzazione misto $l_1-l_2$.
\item SSEH: un modello gerarchico dove $h^{[i]}$ è un processo Gaussiano con covarianza $\lambda_iK(s,t)$ e il fattore di scala $\lambda_i$ ha una distribuzione esponenziale per favorire la sparsità sullo spazio dei fattori di scala.

\end{enumerate}

A questo punto si rende necessario introdurre della notazione aggiuntiva per poter spiegare gli algoritmi che verranno utilizzati. Si definisce dunque
\begin{equation}
y_t^+:=
\begin{bmatrix}
y_t \\ 
\vdots \\
y_{t+N-1}
\end{bmatrix}
\quad
e_t^+:=
\begin{bmatrix}
e_t\\
\vdots  \\
e_{t+N-1}
\end{bmatrix}
\quad
h:=
\begin{bmatrix}
h^{[1]} \\
\vdots \\
h^{[m]}
\end{bmatrix}
,
\label{def yeh}
\end{equation}
dove $h^{[i]}\in l_1(\mathbb{Z}^+)$, $i=1,\ldots,m$. Il predittore (\ref{predittore_gay}) può essere così riscritto
\footnote{Il prodotto tra matrici semi-infinite va inteso come il limite della sequenza finita, che sotto l'ipotesi $h^{[i]}\in l_1{\mathbb{Z}^+}$, è ben definito.}
come:
\begin{equation}
y_t^+=\underbrace{[A_{t1}\,\ldots \, A_{tm}]}_{:=A_t}h+e_t^+
\label{predittore riscritto}
\end{equation}
con $A_{ti}\in \mathbb{R}^{N \times \infty}$, $i=1,\ldots,m$ definiti da:
\begin{equation}
\begin{array}{rcl}
A_{ti}^{[jk]} & := & u_{t-j-k}^{[i]}\,, \quad i=1,\ldots,m-1\ ,\quad j,k \in \mathbb{Z}^+ \\
A_{tm}^{[jk]} &:=& y_{t-j-k},\quad j\,,k \in \mathbb{Z}^+
\end{array}
\label{gli A}
\end{equation}
e il problema di identificazione può essere dunque visto come la stima di $h$ in (\ref{def yeh}),(\ref{predittore riscritto}) soggetto al vincolo che $h^{[i]}\in l_1(\mathbb{Z}^+)$, $i=1,\ldots,m$. Come già detto lo stimatore dovrà automaticamente selezionare tra $\{u^{[1]},\ldots,u^{[m]},y\}$ le variabili necessarie alla predizione di $y$ e quali no e questo farà in modo che alcune risposte impulsive $\hat{h}^{[i]}$ siano poste esattamente a zero.






\section{Predizione in RKHS e approssimazione finito-dimensionale}
\label{sec:PredizioneInRKHSEApprossimazioneFinitoDimensionale}


Si assume che le risposte impulsive $h^{[i]}$ siano (le versioni campionate) di funzioni in $\mathcal{H}_K$, lo spazio di Hilbert riproduttivo associato al nucleo campionato $K$. Questo corrisponde al fatto che le risposte impulsive $h^{[i]}$ sono processi Gaussiani con funzione di covarianza $K$.\\
Il problema di stimare le risposte impulsive $h^{[i]}$ dai dati misurati può essere formulato come il seguente problema di tipo Tikhonov di regolarizzazione:
\begin{equation}
\{\hat{h}^{[i]}\}=\arg \min_{h^{[i]}\in\mathcal{H}_K}\sum_{t=t_0}^N(y_t-\hat{y}_{t|t-1})^2+ \gamma_2^2\left(\sum_{i=1}^{m}\|h^{[i]}\|^2_{\mathcal{H}_K}\right)
\label{eq:Tikhonov}
\end{equation}
soggetto a
\begin{equation}
\hat{y}_{t|t-1}=\sum_{i=1}^{m-1}\left[\sum_{k=1}^{\infty}h^{[i]}_k u^{[i]}_{t-k}\right] +\sum_{k=1}^{\infty}h^{[m]}_k y_{t-k}.
\label{eq:PredictorTikhonov}
\end{equation}
Il parametro $\gamma_2$ è il cosiddetto parametro di regolarizzazione che deve regolare il trade-off tra fit e regolarità di $h^{[i]}$.\\

Per portare questo problema in dimensione finita si modifica il nucleo $K$ in modo che la realizzazione del nucleo modificato soddisfi la risposta impulsiva alla condizione $h^{[i]}_k=0, \ \forall k>J$. È facile vedere che il nucleo (la covarianza, indicata con $K_J$) che soddisfa la condizione è $K_J(h,k)=0, \  \forall(h,k): \ k>J \ \text{oppure}  \ h>J$.\\
Raccogliendo i coefficienti della risposta impulsiva in un vettore colonna si ha
$$ \underline{h}^{[i]}:=[h^{[i]}_1,h^{[i]}_2,\ldots,h^{[i]}_{t_0}]^T, $$
dove $t_0=J$. Il numero $J$ non ha influenza sul trade-off tra bias e varianza, ma serve solo per ragioni di complessità computazionale. Con un piccolo abuso di notazione si indica con $K_J \in \mathbb{R}^{J\times J}$ la matrice simmetrica, definita positiva formato dai nuclei campionati fino al ritardo $J$. Sotto queste condizione la norma $\|h^{[i]}\|^2_{\mathcal{H}_{K_J}}$ può essere scritta come
$$ \|h^{[i]}\|^2_{\mathcal{H}_{K_J}}=(\underline{h}^{[i]})^T K_J^{-1}\underline{h}^{[i]}. $$
Il predittore a un passo $\hat{y}_{t|t-1}$ diventa
\begin{equation}
\hat{y}_{t|t-1}=\sum_{i=1}^{m-1}\left[\sum_{k=1}^{J}h_k^{[i]}u_{t-k}^{[i]}\right] +\sum_{k=1}^{J}h_k^{[m]}y_{t-k}.
\label{eq:PredictorTikhonovApprox}
\end{equation}
Risulta che sotto la restrizione $h^{[i]}\in\mathcal{H}_{K_J}$ la soluzione del problema \ref{eq:Tikhonov} è
\begin{equation}
\arg \min_{\underline{h}^{[i]}\in\mathbb{R}^{t_0}}\sum_{t=t_0}^N(y_t-\hat{y}_{t|t-1})^2+ \gamma_2^2\sum_{i=1}^{m}\left[(\underline{h}^{[i]})^T K_J^{-1}\underline{h}^{[i]}\right]
\label{eq:TikhonovApprox}
\end{equation}
soggetto a \ref{eq:PredictorTikhonovApprox}. Si vede quindi che il problema può essere riformulato come quello di stimare un modello ARX con un peso di regolarizzazione incluso nel problema della stima.


\section{Stabil Spline Group LAR (SSGLAR) algoritm}
\label{sec:StabilSplineGroupLARSSGLARAlgoritm}

Viene ora discusso come si può risolvere il problema \ref{eq:TikhonovApprox} invocando sparsità nei gruppi $h_{t_0}^{[i]}$ usando l'algoritmo GLAR. Per fare questo si assume che i parametri $\theta$ siano stati fissati; senza alcuna informazione a priori questo parametro viene posto uguale a zero.\\
Il problema \ref{eq:TikhonovApprox} è un problema di regressione lineare con una funzione peso in $l_2$ che può essere riscritto nella forma
\begin{equation}
\overline{y}^+_{t_0}=\sum_{i=1}^{m}\overline{A}_{Ji}h_{t_0}^{[i]}+W 
\label{SchwuleForml1}
\end{equation}
avendo definito
\begin{align*}
\overline{y}^+_{t_0} &:= \begin{bmatrix} y^+_{t_0} \\ 0_{1\times (t_0\cdot m)} \end{bmatrix} \\
\overline{A}_{Ji} &:= \left[A_{Ji}\quad v_i\otimes \Lambda\right] \quad \Lambda:=\gamma_2 K_J^{-1/2} \\
v_i &:= [\underbrace{0\ldots 0}_{i-1} 1 \underbrace{0\ldots 0}_{m-i}]^T	
\end{align*}

La selezione degli ingressi come descritto prima può essere eseguito dall'algoritmo GLAR applicato al problema di regressione. Nella prossima sezione si definisce nel dettaglio questo algoritmo che viene chiamato SS-GLAR.

\subsection{Algoritmo SS-GLAR}
\label{sec:AlgoritmoSSGlar}
I passi che compie l'algoritmo SS-GLAR sono:
\begin{enumerate}
	\item stabilire il parametro $\beta$ del nucleo;
	\item stabilire il parametro $\gamma_2$ di \ref{eq:TikhonovApprox} e formare $\overline{A}_{Ji}$ in \ref{SchwuleForml1};
	\item stima di $h^{[i]}$ applicando l'algoritmo GLAR al problema \ref{SchwuleForml1};
\end{enumerate}


\subsection{Stima degli iperparametri}
\label{sec:StimaDeiHyperParametri}
Prima di poter eseguire l'algoritmo appena visto bisogna fissare i seguenti parametri:
\begin{enumerate}
	\item la funzione peso $\gamma_2$ in $l_2$ nell'equazione \ref{SchwuleForml1} (regolarità di $[h^{[i]}_1,h^{[i]}_2,\ldots,h^{[i]}_{t_0}]^T$ nello spazio $\mathcal{H}_{K_{t_0}}$);
	\item il parametro $\beta$ del nucleo (velocità di decrescita del nucleo);
	\item Il numero di blocchi diversi da zero stimati dall'algoritmo GLAR;
\end{enumerate}
Questi parametri possono essere stimati usando il seguente approccio di validazione: disponendo delle misure $\{y_t,u_t\}_{t=1,\ldots,N}$ si suddivide l'insieme in due parti. Il primo $\{y_t,u_t\}_{t=1,\ldots,2N/3}$ viene chiamato l'insieme di identificazione, il secondo con i dati rimanenti viene chiamato l'insieme di validazione. Si esegue l'algoritmo per iperparametri fissati usando l'insieme di identificazione e successivamente si valuta il modello ottenuto sull'insieme di validazione. Dall'insieme degli iperparametri $(\beta\in \mathbb{R}^+ ,\gamma \in \mathbb{R}^+ )$ si sceglie un insieme con un numero finito (possibilmente piccolo) di alternative, sul quale poi si deve testare l'algoritmo. Vengono inoltre testati diversi livelli di sparsità. Alla fine viene scelto il modello identificato che ha dato le miglior prestazioni sul set di validazione in termini di $RMS_1$ (Root Mean Square) dell'errore di predizione a un passo, dove il generico $RMS_k$ ($k$ numero di passi) è definito da:
$$ RMS_k:=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_t-\hat{y}_{t|t-k})^2}. $$
Vengono quindi fissati gli iperparametri e il livello di sparsità con i quali viene rifatta la stima del modello usando l'intero insieme di dati $\{y_t,u_t\}_{t=1,\ldots,N}$.


\subsection{Stable Spline with Exponential Hyperpprior (SSEH) Algoritm}
\label{sec:StableSplineWithExponentialHyperparametriSSEHAlgoritm}

In questa sezione il predittore di equazione \ref{eq:PredictorTikhonov} viene formulato in un contesto Bayesiano assumendo che le risposte impulsive $h^{[i]}$ siano processi Gaussiani con covarianza $\lambda_i K(t,s)$. Serve specificare gli iperparametri $\beta$, $\theta$ e $\lambda_i$ e la varianza di rumore $\sigma^2$ che descrivono interamente le densità di probabilità delle misure di $y$ e le risposte impulsive dei predittori $h^{[i]}$.\\


\subsubsection{Iperprobabilità a priori per gli iperparametri}
\label{sec:HyperprobabilitaAPrioriPerIHyperparametri}

La varianza di rumore $\sigma^2$ viene stimata in un passo preliminare usando un modello ARX per poi supporla nota nella stima Bayesiana. Gli iperparametri $\beta$, $\theta$ e $\lambda_i$ vengono invece modellati come variabili aleatorie indipendenti. $\beta$ ha una densità di probabilità non informativa (cioè inizialmente non si hanno informazioni a priori sul parametro) su $\mathbb{R}^+$, mentre $\theta$ è distribuita uniformemente. Ogni $\lambda_i$ è una variabile aleatoria esponenziale con media $1/\xi \in \mathbb{R}^+$, quindi vale
$$ p(\lambda_i)=\xi exp (-\xi\lambda_i)\chi(\lambda_i\geq0), \qquad i=1,\ldots,m $$
con $\chi$ la funzione indicatrice. $\xi$ può essere interpretata come una variabile aleatoria non informativa su $\mathbb{R}^+$. Tutti gli iperparametri possono essere raccolti in un vettore aleatorio indicato con $\zeta:=[\lambda_1,\ldots,\lambda_m,\theta_1,\theta_2,\beta,\xi]$.\\
Si noti che i parametri $\lambda_i$, che hanno lo stesso ruolo di $\sigma^2 / \gamma^2$ con $\gamma^2$ definito in \ref{eq:Tikhonov}, possono ora essere tutti differenti aumentando così la flessibilità del modello.\\


\subsubsection{Il modello Bayesiano intero}
\label{sec:IlModelloBayesianoIntero}

Sia $A_{ti}\in\mathbb{R}^{N\times\infty}$ come definito in \ref{gli A} in modo che per le uscite future $y^+:=y_1^+$ risulta
$$ y^+=\left[\sum_{i=1}^{m}A_{1i}h^{[i]}\right]+e^+, $$
dove $e^+ := e^+_1$.\\
In pratica $y^-:=[y_0,y_{-1},y_{-2},\ldots]$ non è mai interamente noto e una soluzione consiste nel porre a zero le componenti non note. Viene quindi imposta la seguente approssimazione:
\begin{align*}
	p(y^+,\{h^{[i]}\},y^-,u|\zeta) &\propto \left[\prod_{t=1}^{N}p(y_t|\{h^{[i]}\},y_t^-,\zeta,u_t^-)\right]p(y^-,\{h^{[i]}\},u^-|\zeta)\\
	&\approx \left[\prod_{t=1}^{N}p(y_t|\{h^{[i]}\},y_t^-,\zeta,u_t^-)\right]p(\{h^{[i]}\}|\zeta)p(y^-,u^-)
\end{align*}
La proporzionalità deriva dall'assunzione che il predittore di $u_t$ dato il passato $u_t^-$ e $y_{t+1}^-$ non dipenda da $\zeta$ mentre l'approssimazione deriva dall'assunzione che il passato $u_t^-$ e $y_{t}^-$ porti informazione sulla risposta impulsiva del predittore e sugli iperparametri. Questo modello statistico è descritto dalla rete Bayesiana riportata in Figura \ref{fig:ReteBayesiana} (sinistra).\\
Come visto sopra, per questioni di complessità computazionale bisogna troncare la risposta impulsiva $h^{[i]}$ a una certa lunghezza $t_0$.\\

\subsubsection{Stima degli iperparametri}
\label{sec:StimaDegliHyperparametri}

Il vettore degli iperparametri $\zeta$ viene stimato minimizzando la probabilità marginale, cioè la densità di $y^+$, $\zeta$ e $\{h^{[i]}\}$. Può essere dimostrato (si guardi \cite{Chiuso10}) che la stima si ottiene da
\begin{equation}
 \hat{\zeta}=arg \min_{\zeta}J(y^+;\zeta)
 \label{eq:MinimizZeta}
\end{equation}
con $\xi,\beta>0$, $\lambda_i>0, \text{ per } i=1,\ldots,m$ e $J$ dato da
\begin{align}
	J(y^+;\zeta) &=\frac{1}{2} log (det[2\pi V[y^+]])+\frac{1}{2}(y^+)^T(V[y^+])^{-1}y^++\xi\sum_{i=1}^{m}\lambda_i-log(\xi)\\
 V[y^+] &=\sigma^2I_N+\sum_{k=1}^{m}\lambda_i A_{1i} KA_{1i}^T\,.
 \label{eq:SchwuleForml2}
\end{align}

\begin{figure}[t]
    \begin{center}
      \includegraphics[width=0.7\textwidth]{figure/ReteBayesiana}
    \caption[ReteBayesiana]
    {\label{fig:ReteBayesiana}
    Rete Bayesiana che descrive il nuovo modello non parametrico per l'identificazione di sistemi dinamici sparsi, dove $y^-:=[y_{l-1},y_{l-2},\ldots]$ e nel modello ridotto $\lambda:=\lambda_1,\ldots,\lambda_m$.}
    \end{center}
\end{figure}

Una questione molto importante per l'utilizzo pratico di questo procedimento matematico è l'utilizzo di un buon valore iniziale. Nel seguito viene descritto un schema con il quale è possibile ottenere una soluzione sub-ottimale risolvendo il problema di ottimizzazione in $\mathbb{R}^4$ in relazione al modello Bayesiano di Figura \ref{fig:ReteBayesiana} (destra).

\begin{description} 
  \item[(i)] Si ottengono $\{\hat{\lambda_i}\}$, $\hat{\theta}$ e $\hat{\beta}$ risolvendo il problema \ref{eq:MinimizZeta} modificato nella seguente versione:
  $$ \arg \min_{\zeta}[J(y^+;\zeta)-\xi\sum_{i=1}^{m}\lambda_i+log(\xi)] $$
  con $\beta > 0, \ \lambda_1=\ldots=\lambda_m\geq0$
  \item[(ii)] Si fissano $\hat{\xi}=1/\hat{\lambda_1}$ e $\hat{\zeta}:=[\hat{\lambda}_1,\ldots,\hat{\lambda}_m,\hat{\theta}_1,\hat{\theta}_2,\hat{\beta},\hat{\xi}]$. Poi, per $i=1,\ldots,m$ si fissa $\bar{\zeta}=\hat{\zeta}$ tranne la i-esima componente di $\bar{\zeta}$ che viene messa a zero. Se invece $J(y^+;\bar{\zeta}) \leq J(y^+;\hat{\zeta})$ si fissa $\hat{\zeta}=\bar{\zeta}$.
\end{description} 



La procedura descritta al passo (ii) potrebbe non essere adatta nel caso in cui gli ingressi sono altamente correlati, perché potrebbe risultare sensibile all'ordine con cui le componenti di $\hat{\zeta}$ convergono a zero. Per aggirare questo problema si propone di usare l'algoritmo SS-GLAR per effettuare questa selezione. In particolare si può sostituire il passo (ii) con il passo seguente:
\begin{description}
	\item[(ii')] Si fissa $\gamma_2^2=\sigma^2/\hat{\lambda}_1$ e si aggiusta il nucleo visto in sezione \ref{sec:StabilSplineGroupLARSSGLARAlgoritm} con $K(:,:,\hat{\theta},\hat{\beta})$. Si esegue poi l'algoritmo SS-GLAR selezionando il numero di componenti non nulli usando l'approccio di validazione. A questo punto i componenti $\{j_1,\ldots,j_k\}\subseteq\{1,\ldots,m\}$ di $\hat{\zeta}$ vengono annullati, il che equivale a annullare $\hat{h}^{[j]}=0,\ j\in\{j_1,\ldots,j_k\}$.
\end{description}
Il valore $\hat{\zeta}$ finale può essere usato come punto iniziale in problema \ref{eq:SchwuleForml2}.


\subsubsection{Stima delle risposte impulsive del predittore per $\zeta$ noto}
\label{sec:StimaDiRisposteImpulsivePerZetaNoto}

Sotto certi condizioni (vedi \cite{Chiuso10}) si può dimostrare:
$$	\{\hat{h^{[i]}}\}^{m}_{i=1}=arg \min_{\{f^i\in\mathcal{H}\}^{m}_{i=1}} \|y^+ -\sum_{i=1}^{m}A_{1i}f^i\|^2 + \sigma^2\sum_{i=1}^{m}\frac{\|f^i\|^2_{\mathcal{H}_K}}{\lambda^2_i} $$
Inoltre per $k=1,\ldots,m+1$ si ha
$$ \hat{h^{[i]}} = \lambda^2_i K A_{1i}^T c, \quad c=\left(\sigma ^2I_n+\sum_{i=1}^{m+1}\lambda_iA_{1i}KA_{1i}^T\right)^{-1}y^+ $$
Dopo aver ottenute le stime di $\{h^{[i]}\}$ semplice formule possono essere usate per ottenere le risposte impulsive $f$ e $g$ in
\ref{eq:SisDin} e quindi anche il predittore a $k$ passi (per dettagli vedere \cite{Ljung99}).










\end{document}





