\documentclass{article}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{graphicx}

\title{Processi Gaussiani per la regressione e Kriging}
\author{Maria Veronica Vinattieri}

\begin{document}

\maketitle

\section{Introduzione}
\section{I processi gaussiani}
Teoria sui processi gaussiani e sul modello geostatistico/spaziale. Applicazione tramite il metodo di previsione kriging. \\\\
I processi gaussiani nascono dal problema di imparare la relazione tra input e output da dati empirici. \\\\ Quando l'output assume valori nel continuo allora la relazione tra questo e l'input è la regressione, invece quando l'output è discreto allora il metodo per imparare la relazione è la classificazione. \\\\ Si definisce l'input (o il vettore di input) $x$ e l'output $y$. \\\\ Il problema è di tipo induttivo: dalle osservazioni $D\left\{(x_i,y_i)|i=1,...,n\right\}$ si vuole prevedere il valore delle $y^{*}$ tale per cui il relativo valore di input $x^{*}$ non è stato osservato. Viene definita quindi una funzione $f$ che fa previsioni per ogni input $x$.\\\\ La $f$ può essere scelta restringendosi a un'unica famiglia di funzioni ma incorre il problema dell'overfitting, in quanto la $f$ potrebbe essere troppo fedele ai dati che ha osservato.\\ Oppure può essere assegnata ad ogni $f$ una probabilità a priori, tanto più alta quanto è più probabile che quella funzione sia adatta al caso considerato e scelta una $f$ a posteriori. Il problema con questo approccio è che vi sono infinite funzioni tra le quali scegliere.\\\\ La scelta più giusta sembra essere il processo gaussiano. Infatti tramite esso viene scelta una $f$ che può essere vista come un lungo vettore di valori $f(x)$ per ogni particolare input $x$. Il punto di forza del processo gaussiano è la semplicità unita alla trattabilità computazionale. \\\\ La funzione $(f(x): x \in R^k)$ è un processo gaussiano se la distribuzione congiunta delle $f(x_i)$ è distribuita come una gaussiana multivariata. \\\\
Un processo gaussiano è totalmente descritto da una funzione della media $\mu(x)$ e una funziona della covarianza $K(x',x)$. \\\\
Per facilitare la notazione, il processo gaussiano sarà d'ora in poi scritto $X_t$ con $t = (t_1,t_2,...,t_n)$. Il processo gaussiano è continuo nel tempo e nello spazio. \\\\ Per definizione la funzione $\mu_X:[0, \infty] \rightarrow R$ definita per ogni $t \ge 0$ da $\mu(X_t)=E(X_t)$ è chiamata funzione della media del processo $X_t$. \\\\ Per definizione la funzione $K_X:[0, \infty] \rightarrow R$ definita per ogni $s,t \ge 0$ nel seguente modo $C_X(s,t) = Cov[X_s,X_t]$ è chiamata funzione dell'autocovarianza del processo $X_t$.
Allora un processo gaussiano è governato dalla legge:
\begin{align*}
    p(f(x),f'(x)) \sim N(\mu, \Sigma)
\end{align*}
A seconda della forma di $\mu_X$ e $\Sigma_X$ vengono definiti processi gaussiani diversi. Un esempio tipico è quello con $\mu_X=\mu$ e $\Sigma_X=\sigma^2 I$. \\\\
Il processo gaussiano non è altro che la generalizzazione della distribuzione di densità gaussiana. Il PG governa le proprietà delle funzioni, come una densità di probabilità governa le proprietà della variabili. \\\\
Imparare la relazione input-output tramite i processi gaussiani è esattamente come trovare delle proprietà convenienti per la funzione di covarianza.
\section{Regressione}
\subsection{Regressione lineare Bayesiana}
Si considera il modello di regressione lineare \begin{align*}
f(x) = x'\beta  \ \rightarrow \ y = f(x) + \epsilon
\end{align*}
dove le $x$ sono gli input, $\beta$ i parametri, e $y$ è l'attributo d'interesse, nonchè ouput. $f$ è la funzione che trasforma gli input negli output con un'aggiunta dell'errore di misura $\epsilon$ distribuito normalmente $N(0, \sigma_n^2$. \\\\ La likelihood, densità di probabilità delle osservazioni dati i parametri, del modello è \begin{align*}
p(y|X,\beta) = \dfrac{1}{(2\pi\sigma_n^2)^{n/2}} exp(-\dfrac{1}{2\sigma_n^2}|y-X'\beta|^2) \sim N(X'\beta, \sigma_n^2 \mathbf{I})
\end{align*}
Seguendo un approccio Bayesiano si definisce una a priori sui parametri \begin{align*}
\beta \sim N(0, \Sigma_p)
\end{align*}
E si ricava l'a posteriori \begin{align*}
&p(\beta|y,X) = \dfrac{p(y|X,\beta)p(\beta)}{p(y|X)} \\
&p(\beta|y,X) \sim p(y|X,\beta)p(\beta) \ \mbox{dove} \ p(y|X,\beta)=\int p(y|X,\beta)p(\beta)d\beta \\
&...\\
&p(\beta|y,X) \sim N(\bar{\beta} = \dfrac{1}{\sigma_n^2}A^{-1}Xy, A^{-1}) \\
&\mbox{dove} \ A=\sigma_n^{-2}XX'+\Sigma_p^{-1}
\end{align*}
L'a posteriori sui parametri è ancora una Normale. \\\\
Quindi si può fare previsione pesando le probabilità sui valori previsti con le probabilità a posteriori sui parametri \begin{align*}
&p(f_*|x_*,X,y) = \int p(f_*|x_*,\beta)p(\beta|X,y)d\beta \\
&p(f_*|x_*,X,y) \sim N(\dfrac{1}{\sigma_n^2}x'_*A^{-1}Xy, x'_* A^{-1}x_*)
\end{align*}
La distribuzione predittiva è ancora una Normale. \\\\
Per superare il problema dellea limitata espressività del modello lineare Bayesiano, si può proiettare gli input $x$ su uno spazio di dimensione maggiore $\phi(x)$ e lavorare direttamente con esso. Risulta che \begin{align*}
f(x) = \phi(x)'\beta
\end{align*}
e la distribuzione predittiva è \begin{align*}
&f_*|x_*,X,y \sim N(\dfrac{1}{\sigma_n^2}\phi(x_*)'A^{-1}\phi y, \phi(x_*)A^{-1}\phi(x_*)) \\
&\mbox{dove} \ \phi = \phi(X) \ \mbox{e} \ A=\sigma_n^{-2}\phi \phi'+\Sigma_p^{-1} \\
&f_*|x_*,X,y \sim N(\phi(x_*)\Sigma_p \phi(K+\sigma_n^2I)^{-1} y, \\ &\phi(x_*)\Sigma_p\phi(x_*)-\phi(x_*)\Sigma_p\phi(K+\sigma_n^2I)\phi'\Sigma_p\phi(x_*))
\end{align*}
dove $K=\phi'\Sigma_p\phi$ e $k(x,x')=\phi(x)\Sigma_p\phi(x')$.\\\\ Allora si definisce un processo gaussiano $f(x)\sim GP(m(x),k(x,x'))$ totalmente specificato da una funzione della media \begin{align*}
m(x) = E[f(x)]=\phi(x)E[\beta]=0
\end{align*}
e una funzione di covarianza \begin{align*}
k(x,x')= &E[(f(x)-m(x))(f(x')-m(x'))] \\ =&\phi(x)E[\beta,\beta']\phi(x')=\phi(x)\Sigma_p\phi(x')
\end{align*} La funzione di covarianza fra output è scritta in funzione degli input. Il valore della covarianza è indirettamente proporzionale alla distanza spaziale fra gli input. \\\\
Il processo gaussiano applicato ad un numero di punti $X_*$ è \begin{align*}
f_* \sim N(0,K(X_*,X_*))
\end{align*}
Una volta definito il processos gaussiano $f$ si ricavano 2 diverse distribuzioni predittive: \begin{itemize}
\item \textbf{previsione su dati senza rumore}: la distribuzione congiunta tra i dati $f$ e i dati da prevedere $f_*$ è \begin{align*}
\begin{bmatrix} f \\ f_* \end{bmatrix} \sim N \left( 0,\begin{bmatrix}
K(X,X) \ \ K(X,X_*) \\
K(X_*,X) \ \ K(X_*,X_*) \end{bmatrix} \right)
\end{align*}
Allora la predittiva è la distribuzione congiunta Gaussiana a priori condizionata alle osservazioni, e risulta: \begin{align*}
f_*|X_*,X,f \sim &N(K(X_*,X)K(X,X)^{-1}f, \\
&K(X_*,X_*)-K(X_*,X)K(X,X)^{-1}K(X,X_*))
\end{align*}
\item \textbf{previsione su dati con rumore}: con il modello $y=f(x)+\epsilon$ la covarianza risulta $\mbox{cov}(y)=K(X,X)+\sigma_n^2 I$. \\\\La distribuzione congiunta  dei valori osservati e i valori da predirre è \begin{align*}
\begin{bmatrix} y \\ f_* \end{bmatrix} \sim N \left( 0,\begin{bmatrix}
K(X,X)+\sigma_n^2I \ \ K(X,X_*) \\
K(X_*,X) \ \ K(X_*,X_*) \end{bmatrix} \right)
\end{align*}
Allora la predittiva risulta: \begin{align*}
f_*|X_*,X,f \sim &N(K(X_*,X)[K(X,X)+\sigma_n^2I]^{-1}y, \\
&K(X_*,X_*)-K(X_*,X)[K(X,X)+\sigma_n^2I]^{-1}K(X,X_*))
\end{align*}
Per un solo punto $x_*$ in cui si vuole fare previsione l'equazione sopra è ridotta a\begin{align*}
&\bar{f}_* = k'_*(K+\sigma_n^2I)^{-1}y\\
&V[f_*]=k(x_*,x_*)-k'_*(K+\sigma_n^2I)^{-1}k_*
\end{align*} dove $k_*=k(x_*)$.
\end{itemize}

\section{Classificazione}
Date le osservazioni $D=\left\{(x_i,y_i)|i=1,...,n\right\}$ tali per cui le $y_i$ possono assumere valori discreti $C_1,...,C_c$, si vuole assegnare ogni nuovo input $x^{*}$ a una delle classi. \\\\ La classificazione è definita probabilistica, quindi ogni previsione equivale alla assegnazione ad una classe probabilistica. Vi è quindi una certo grado di incertezza su ogni previsione.\\\\ Nella classificazione, al contrario della regressione, la likelihood è non-gaussiana ma il processo a posteriori è un'approssimazione gaussiana. \\\\ Il problema della previsione si traduce nel fare inferenza su $p(y|x)$. Nell'ambito della classificazione esistono 2 approcci distinti per fare inferenza sulla probabilità predittiva. \begin{itemize}
    \item Generativo: la probabilità è ricavata tramite la formula di Bayes:\begin{align*}
        p(y|x) = \dfrac{p(y)p(x|y)}{\sum_{c=1}^{C}p(C_c)p(x|C_c)}
    \end{align*} \\ La scelta per la densità condizionata può essere una gaussiana $p(x|C_c) \sim N(\mu_c, \Sigma_c)$. Inoltre può essere applicativo un approccio bayesiano ponendo delle a priori appropriate sulla media $\mu_c$ e la covarianza $\Sigma_c$.
    \item Discriminativo: la probabilità predittiva $p(y|x)$ viene modellata direttamente. \\In questo caso si deve scegliere una funzione di risposta (l'inversa della funzione link $g(E(y_i|x_i)) = \eta_i$) che racchiude la risposta $y$ nel range $[0,1]$: \begin{enumerate}
    \item la funzione logit: $ p(C_1|x) = logit(\beta'x) = \dfrac{1}{1+exp(-\beta'x)}$
    \item la funzione probit: $p(C_1|x) = probit(\beta'x) =\phi(\beta'x)$
\end{enumerate}
\end{itemize}
Si segue l'approccio discriminativo per la classificazione tramite i processi Gaussiani.

\subsection{Teoria delle decisioni per la classificazione}
L'approccio probabilistico alla classificazione combina l'implementazione della probabilità a posteriori sulle funzioni con la funzione di perdita della teoria sulle decisioni, proprio per prendere una decisione sulla classificazione. La teoria delle decisione si serve di una funzione di perdita. La funzione di perdita $L(c,c')$ riguarda l'alternativa scelta $c'$ quando la vera classe è $C_c$. \\\\
La perdita attesa nel prendere la decisione $c'$ dato x è $R_{L(c'|x)} - \sum_c L(c,c')p(C_c|x)$ e la scelta ottimale $x^{*}$ è quella che minimizza la $R_{L(c'|x)}$. \\\\ Il classificatore ottimale è conosciuto come il classificatore di Bayes. Una $x^{*}$ non osservata viene assegnata ad una delle possibili $R_1,...,R_c$ regioni di decisione, tali per cui la regione $R_c$ è assegnata alla classe $C_c$. Nel caso $C >2$ la $x^{*}$ viene classificata quando la differenza tra $max(p(C_j|X^{*})$ e la seconda probabilità più alta $\ge \theta$ con $\theta \in [0,1]$. \\\\ Viene quindi approssimata la funzione di classificazione come quella che minimizza il rischio \begin{align*}
    R_c(c)=\int L(y,c(x))p(y,x)dydx
\end{align*}
dove $c(x)$ è la funzione di classificazione che assegna $x$ ad una delle $C$ classi.

\subsection{Modelli lineare per la classificazione}
Viene considerato il caso con 2 classi quindi $y=-1, y=+1$. La funzione di verosimiglianza è: \begin{align*}
    p(y=+1|x,w) = \sigma(x'w)
\end{align*}
dove w è il vettore dei pesi, e $\sigma()$ è una funzione monotona crescente definita da $R$ a $[0,1]$. Allora $\sigma(x)$ può essere pari a $\lambda(x)$: il modello è lineare logit; oppure $\sigma(x)=\phi(x)$ allora il modello è lineare probit. \\\\ Per fare previsioni su $x_*$ da un set osservato $D$ si ha \begin{align*}
    p(y_* =+1|x_*, D) = \int p(y_*=+1|w,x_*)p(w|D)dw
\end{align*}
dove $p(y_*=+1|w,x_*)= \sigma (x'_* w)$ e $p(w|D)=p(w|X,y)$ dove, data l'a priori $w \sim N(0, \Sigma_p)$, il logaritmo della distribuzione a posteriori non normalizzata è definito \begin{align*}
    log(p(w|X,y)) = \dfrac{1}{2} w' \Sigma^{-1}_p w + \sum_{i=1}^n log (\sigma(y_ix'_iw))
\end{align*}
Invece nel caso multi-classe la probabilità di classificazione è \begin{align*}
    p(y=C_c|x,W) = \dfrac{exp(x'w_c)}{\sum_c' exp(x'w_c')}
\end{align*}
dove $W$ è la matrice dei pesi per ogni classe, infatti $w_c$ è il peso della classe generica $C_c$.

\subsection{L'utilizzo dei processi gaussiani per la classificazione} Viene ora considerata la funzione $f(x)$ con un processo gaussiano, quindi assegnata una GP prior adatta. Si ricava quindi una prior per la probabilità predittiva che è l'obiettivo di interesse $\pi(x)=p(y=+1|x)=\sigma(f(x))$. Inoltre viene definita una GP prior anche per il vettore dei pesi $w$. \\\\ Viene prima implementata l'a posteriori del processo gaussiano $f_* = f(x_*)$ sulla nuova osservazione $x_*$ \begin{align*}
    p(f_*|X,y,x_*)=\int p(f_*|X,y,x_*,f)p(f|X,y)df
\end{align*}
dove $p(f|X,y) = \dfrac{p(y|f)p(f|X)}{p(y|X)}$. Dopo di che viene calcolata la probabilità predittiva \begin{align*}
    \pi_* = p(y_*=+1|X,y,x_*) = \int \sigma(f_*)p(f_*|X,y,x_*)df_*
\end{align*}
L'integrale non è risolvibile analiticamente. Vengono quindi utilizzati degli algoritmi di approssimazione: l'approssimazione di Laplace e metodo della propagazione attesa (EP method). Anche algoritmi di approssimazione quali MCMC possono funzionare.

\section{Il modello geostatistico}
Il modello geostatistico è espresso \begin{align*}
Y_i = S(x_i)+U_i
\end{align*} dove $S() = \left\{S(x), \forall x \in R^k \right\}$  con $x_i$ locazioni spaziali dei punti di osservazione, è un processo gaussiano spaziale, e $U_i$ è l'errore di misurazione. \\\\
Allora le componenti sono distribuite nel modo seguente:
\begin{align*}
    &S = (S(x_1),...,S(x_n)) \sim MVN(\mu, K) \\
    &U = (U_1,...,U_n) \sim MVN(0, \tau^2 I)\\
    &Y \sim MVN(\mu, K+\tau^2 I)
\end{align*}
Lo scopo dei modelli geostatistici è quello di fare inferenza su un'area d'interesse $A \in R^k$ a partire da dei dati osservati sull'area stessa. Il processo gaussiano si assume come modello di base utilizzato dalla geostatistica per fare inferenza. \\\\ Innanzitutto si vuole imparare la relazione che esso ha con la variabile di risposta (processo di misura) $Y_i$, con $i=1,...,n$. \\\\
Proprietà che assume per il modello geostatistico: \begin{enumerate}
    \item \textbf{Campionamento non preferenziale}: i punti d'osservazione non dipendono dalla natura del processo di misura. La localizzazione dei punti campionati ($x_i$)è indipendente da cosa viene micurato nei punti ($y_i$).
    \item \textbf{Stazionarietà debole}: i momenti primo e secondo del segnale sono finiti e costanti e la funzione di covarianza $K$ del segnale dipende unicamente dagli input. \begin{align*}
    &E[S(x)] = \mu <  \infty, \ E[S(x)^2] < \infty \ \ \ \forall x \in A\\
    &Cov(S(x_i),S(x_j)) = K(c_i-x_j) \ \ \ \forall x_i,x_j \in A \\
    &Var(S(x))=Cov(S(x),S(x)) = K(0) = \sigma^2 \ \ \ \forall x \in A
    \end{align*} Allora anche il processo stocastico definito dagli incrementi $D_u(x)=S(x)-S(x-h)$ è debolmente stazionario.
    \item \textbf{Isotropia}: La funzione di covarianza del segnale è, definita la distanza fra i punti $h=x_j-x_i$ \begin{align*}
    K(x_i-x_j)=K(h)=K(||h||), \ \ \ \forall x_i,x_j \in A
    \end{align*} ciò significa che la covarianza fra 2 punti dipendende unicamente dalla distanza e non dalla direzione. Soddisfa quindi la proprietà di semidefinita positività tale per cui $\sum_{i,j=1}^n a_ia_jK(x_i,x_j) \ge 0$.
    \item \textbf{Indipendenza}: le variabili $Y_i$ sono indipendenti, condizionatamente al PG spaziale, e hanno media costante \begin{align*}
    E[Y_i|S] = \mu(x_i) = \mu \ \ \ \forall i \in 1,...,n
    \end{align*}
\end{enumerate}
Assunte le proprietà sopra, le $Y_i$ si distribuiscono normalmente con \begin{align*}
E[Y_i|S] = S(x_i), \ \ Var(Y_i|S)= \sigma^2 \ \ \ \forall i \in 1,...,n
\end{align*}
\subsection{Variogramma}
Il variogramma misura la dipendenza spaziale di secondo ordine fra coppie di punti: più la vicinanza è stretta più essi saranno simili se vale una forma di dipendenza spaziale.
Il processo gaussiano S ha varianza costante $\sigma^2$ e funzione di covarianza (o covariogramma) che rispetta le seguenti proprietà: \begin{align*}
&K(0) \ge 0 \\
&K(x) = K(-x) \\
&|K(x)| \le K(0)
\end{align*} Allora viene definita la funzione di correlazione simmetrica e non negativa nell'origine: \begin{align*}
\rho(h) = \dfrac{K(h)}{\sigma^2}
\end{align*}
Da qui viene definito il semivariogramma del PG \begin{align*}
\gamma(x_i,x_j) = \dfrac{1}{2}Var(S(x_i)-S(x_j)).
\end{align*} Il variogramma $2\gamma$ è il doppio del semivariogramma. Il variogramma rispetta le proprietà seguenti: \begin{align*}
&\gamma(0) = 0\\
&\gamma(h) = \gamma(-h)\\
&\gamma(h) \ge K(0)\\
&\sum_{i,j=1}^n k_ik_j\gamma(x_i-x_j) \le 0, \ \sum_{j=1}^n k_j=0
\end{align*}
Valgono le seguenti relazioni tra il variogramma, la covarianza, la varianza e la correlazione del segnale: \begin{align*}
&\gamma(h)=K(0)-K(h) = \sigma^2-K(h) \\
&\gamma(h) = \sigma^2(1-\rho(h))
\end{align*}
Dalle formule sopra si può notare che il variogramma tende alla varianza del processo nel momento in cui la correlazione spaziale è nulla. pag 24. \\\\
La scelta del modello per il variogramma (e di conseguenza per il covariogramma) può essere fatta seguendo un modello per il variogramma empirico: \begin{itemize}
\item \textbf{Modello sferico}: \begin{align*}
\gamma(h) = &\sigma^2(\dfrac{3h}{2r}-\dfrac{h^3}{2r^3}) h \le r \\ &\sigma^2 \ \mbox{altrimenti}
\end{align*}
\item \textbf{Modello esponenziale}: \begin{align*}
&\gamma(h) = \sigma^2(1-exp(-\dfrac{3h}{r}) \\
&\rightarrow K(h)=\sigma^2(exp(-\dfrac{3h}{r}))
\end{align*}
\item \textbf{Modello gaussiano}: \begin{align*}
\gamma(h) =\sigma^2(1-exp(-\dfrac{3h^2}{r^2}))
\end{align*}
\end{itemize}
dove $\sigma^2$ è il valore verso cui tende il variogramma detto sella, e $r$ è il range cioè la distanza alla quale il variogramma comincia ad appiattirsi sulla sella.
\subsection{Kriging}
Il kriging è l'applicazione dei processi gaussiani al modello geospaziale per la previsione della misura d'interesse in punti in cui non è stata osservata, a partire dai dati misurati. Si assume che la misura d'interesse $Y_i$ sia continua nello spazio, quindi i punti che si trovano a distanza minore avranno un valore più simile rispetto ai punti più lontani. In generale la previsione spaziale di $T=\tau(S())$ in funzione del segnale, dove $\tau$ è la funzione identità, è definita come $\hat{T}=t(Y)$ dove Y è in funzione del segnale stesso. L'obiettivo è quindi quello di fare previsione minimizzando l'errore quadratico medio $MSE(\hat{T})=E[(T-\hat{T})^2]$.\\\\ Il valore incognito $\hat{T}$ è combinazione lineare dei punti osservati. I pesi assegnati alle osservazioni note dipendono dalla relazione spaziale tra essi e il punto incognito. I pesi sono ricavati dal variogramma, funzione che espone il grado di dipendenza spaziale tra 2 punti, cioè l'autocorrelazione tra $Y(s_i)$ e $Y(s_j)$, dove $s_i, s_j$ sono le coordinate dei punti osservati. \\\\
Si sceglie di lavorare con la componente del processo gaussiano, allora $Y(s_i)$ coincide con un processo gaussiano spaziale \begin{align*} \left \{ Y(s_i), \forall s_i \in \mathbf{R}^k\right\} \end{align*}
Si consideri quindi il modello geostatistico \begin{align*}
    &Y(s) = S(s) = x'\beta+u(s) \end{align*}
    con $x$ le covariate inserite nel modello e $\beta$ i rispettivi coefficienti.\\\\
Si conoscono 3 tipi di kriging: \begin{itemize}
    \item Semplice: quando $\mu(s)$ è nota.
    \item Ordinario: quando $\mu(s)$ è incognita ma costante.
    \item Universale: quando $\mu(s)$ è incognita e varia al variare delle p covariate inserite nel modello.
\end{itemize}

\subsection{Kriging semplice}
Dato il modello $Y(s) = x'\beta+u(s)$, la media del segnale $\mu(s)$ e il coviarogramma K sono note. \\\\ La stima di una nuova osservazione è pari alla media $\mu(s)$ in addizione all'errore stimato $\hat{u}$ nel nuovo punto, si ipotizza infatti che i residui siano legati alla localizzazione dei punti osservati. Allora per il kriging semplice la quantità da prevedere è $T=u(s)$. \\\\
Il residuo stimato per ogni punto è \begin{align*}
    u(s_i) = y_i-\mu(s)=y_i - x'\hat{\beta}
\end{align*}
Per una nuova osservazione $s$ si vuole stimare l'errore $\hat{u}(s)$ e aggiunto alla media costante del segnale per ottenere la previsione di $\hat{y}(s)$. \\\\ La stima del nuovo residuo è una media pesata dei residui osservati: \begin{align*}
    \hat{u}(s)=\sum_{i=1}^m \lambda_i(s)u(s_i)
\end{align*}
Quindi il previsore del kriging semplice risulta: \begin{align*}
    &\hat{S}(s) = \mu(s) + \sum_{i=1}^m \lambda_i(s)u(s_i) \\
    &=\mu(s) + \sum_{i=1}^m \lambda_i(s)(y_i - \mu(s)) \\
    &=(1-\lambda'(s))\mu(s)+\lambda'(s)y
\end{align*}
Il previsore bilancia le informazioni date dalla media nota e quelle ricavate dai dati osservati. \\\\
I pesi $\lambda_i(s)$ vengono ricavati dalla minimizzazione del MSE sulla componente residua \begin{align*}
    E[(\hat{u}(s)-u(s))^2]
\end{align*}
che porta al risultato in forma vettoriale $\lambda(s)=K^{-1}k(s)$ dove K è il covariogramma fra tutti gli $(s_i, s_j)$ e $k(s)$ è il vettore delle covarianze $K(s,s_i)$ che venogono assunte tra ogni punto $s_i$ e la nuova osservazione $s$.\\\\ Allora la stima del residuo risulta $\hat{u}(s) = \lambda_i(s)'u_i=K^{-1}k(s)u$. \\\\
Il minimo valore per l'MSE = $E[(\hat{u}(s)-u(s))^2]$ risulta essere pari all'errore quadratico medio di previsione (o varianza del kriging) $\sigma^2_e = \sigma^2-k'(s)K^{-1}k(s)$ con $\sigma^2 = VAR[u(s)]=K(s,s)$.

\subsection{Kriging ordinario}
Il kriging ordinario è utilizzato quando si lavora localmente. La scelta del modello per il variogramma in questo caso è fondamentale. PERCHE?. \\\\
La media del segnale non è nota ma costante su tutta la superficie d'interesse. La  previsione viene fatta direttamente sul processo di misura \begin{align*}
    \hat{S}(s)=\sum_{i=1}^n\lambda_i(s)Y(s_i)
\end{align*}
I pesi si scelgono minimizzando l'errore quadratico medio $E[(\hat{S}(s)-S(s))^2]$ sotto il vincolo $\lambda'I=1$ tali per cui la media del segnale è $\mu(s)$. I pesi risultano essere $\lambda(s) = K^{-1}k(s)$ dove $k(s) = K(s, s_i)$ il vettore della covarianza tra il nuovo punto da stimare $s$ e i punti osservati $s_i$. \\\\ Il minimo valore corrisponde alla varianza di kriging pari a: $\sigma_e^2 = \sigma^2-k'(s)K^{-1}k(s)$ che coincide a quello trovato con il kriging semplice.

\subsection{Kriging universale}
La media del segnale è ignota e dipende dal numero p di covariate $x$ inserite nel modello $\mu(s) = x' \beta $.\\\\ La matrice K in questo caso ha dimensioni superiori: $(n+p) \times (n+p)$ perchè contiene i valori di tutte le covariate $x_j \ \mbox{con} \ j \in 1,...,p$  nei punti osservati $s_i$. Così anche il vettore delle covarianze fra i punti osservati $s_i$ e il nuovo punto $s$ in quanto contiene tutti i valori delle covariate nel punto $s$: $x_j(s)$. Il vettore dei pesi $\lambda(s)$ ha dimensioni $(n+p) \times 1$. \\\\ I risultati coincidono con quelli ottenuti nel kriging ordinario: \begin{align*}
    &\hat{S}(s)=\sum_{i=1}^n\lambda_i(s)Y(s_i) \\
    &\mbox{dove} \ \lambda(s) = K^{-1}k(s) \  \mbox{con} \ k(s) = K(s, s_i)
\end{align*}
\\\\
Fenomeni rilevanti nel processo previsionale: \begin{itemize}
    \item \textbf{De-clustering}: i pesi delle osservazioni dipendono dalla distanza tra il punto in cui si vuole prevedere il segnale e i punti già osservati. Da una parte i punti più vicini hanno un peso maggiore perchè hanno una locazione simile al nuovo punto d'interesse. Allo stesso tempo perdono la loro importanza quando sono troppo vicini perchè si ipotizza non possano apportare informazioni aggiuntive rispetto a una singola osservazione.
    \item \textbf{Effetto schermo}: quando 2 osservazioni sono (quasi) spazialmente collineari, a una viene assegnato un peso maggiore rispetto all'altra, sebbene siano posizionate esattamente alla stessa distanza dal nuovo punto $x^*$. Può succedere che a uno venga assegnato un peso positivo e all'altro punto un peso negativo. Questo effetto può essere risolto tramite un vincolo sui pesi stessi.
\end{itemize}
Il minimo valore corrisponde alla varianza di kriging pari a: $\sigma_e^2 = \sigma^2-k'(s)K^{-1}k(s)$. Questa volta il valore differisce da quelli trovati con gli altri tipi di kriging perchè il covariogramma ha un'altra composizione.
\section{Applicazione}
I dati sono stati raccolti in una pianura inondata dal fiume Meuse vicino al villaggio di Stein in Olanda. Le osservazioni in totale sono 155. Sono stati rilevate le coordinate dei punti d'osservazione: coordinate orientali e coordinate settentrionali in RDH, l'unità di misura topografica olandese. Nell'area d'interesse sono state misurate le concentrazioni di metalli pesanti quali cadmio, rame, piombo e zinco. Inoltre è stato misurato per ogni locazione la distanza orizzontale in metri dal fiume e l'elevazione in metri rispetto al letto del fiume.
\begin{figure}[h]
\includegraphics[width=\textwidth]{plot1.jpeg}
\caption{concentrazione di zinco (ppm)}\label{fig:1}
\end{figure}
\\\\ Il modello geostatistico ha come variabile di risposta $Y(s_i)$ che misura la concentrazione di zolfo in parti per milione - \emph{figura 1} - e $s_i$ le coppie di coordinate dei punti di osservazione \begin{align*}
Y(s_i) = x'(s_i)\beta + u(s_i)
\end{align*}
Inizialmente la concentrazione di zinco viene trasformata in scala logaritmica per semplicità e maggiore chiarezza nella visualizzazione dei dati stessi: la distribuzione in scala log è infatti più concentrata rispetto alla lineare \begin{align*}
&log(Y(s_i)) = x'(s_i)\beta + u(s_i)\\
&Y(s_i) = exp(x'(s_i)\beta+u(s_i))
\end{align*}
Infatti lo zinco ha i seguenti quantili e un range nell'intervallo $[113.0, \ 1839.0]$.
\begin{align*}
&Min. &1st Qu. \ \ \ \ \ \ \ &Median &Mean \ \ \ \ \ \ \ &3rd Qu. &Max.\\
&113.0 &198.0 \ \ \ \ \ \ \ &326.0 &469.7 \ \ \ \ \ \ \ &674.5 &1839.0
\end{align*}
Mentre lo zinco in scala log varia nell'intervallo $[4.727, \ 7.517]$.
\begin{align*}
&Min. &1st Qu. \ \ \ \ \ \ \ &Median &Mean \ \ \ \ \ \ \ &3rd Qu. &Max.\\
&4.727   &5.288 \ \ \ \ \ \ \ &5.787   &5.886 \ \ \ \ \ \ \ &6.514  &7.517
\end{align*}
L'obiettivo è quello di prevedere le concentrazioni di zinco nei punti non osservati. Per la previsione non vengono considerate ulteriori covariate, ma unicamente la dipendenza spaziale del processo di misura: le concentrazioni di zolfo in punti non osservati verranno totalmente previste tramite le concentrazione di zolfo in punti vicini. \\\\ La scelta di un modello senza trend sulla media porta ad applicare un kriging ordinario per la previsione.
\subsection{Variogramma}
Il modello sul variogramma può essere sferico, esponenziale o gaussiano. Viene analizzato l'andamento del variogramma in funzione della distanza tra coppie di punti (\emph{figura 2}). \begin{figure}[h]
\includegraphics[width=\textwidth]{plot2.jpeg}
\caption{modello per il variogramma}\label{fig:3}
\end{figure}
\\\\Si nota come per il modello sferico il variogramma si appiattisce sulla sella in modo più deciso rispetto agli altri modelli. \\\\
La sella e il range per i modelli risultano:
\begin{center}
\begin{tabular}{|l|l|l|}
\hline
Modello &$\sigma^2$ &range\\
\hline
sferico &0.5906         &897.0209 \\
esponenziale &0.7186    &449.7668 \\
gaussiano   &0.4975     &386.5347 \\
\hline
\end{tabular}
\end{center}
I modelli rendono risultati simili per quanto riguardo il sill ma hanno una differenza importante nel range. \\\\
Per valutare la qualità della previsione utilizzando i diversi modelli viene applicato il metodo della cross validation. Il dataset viene diviso in training set e test set. Tramite le informazioni fornite dal training set si fanno previsioni sulle coordinate dei punti inclusi nel test set. In seguito vengono confrontate le vere concentrazioni di zinco con quelle previste tramite MAE e MSE \begin{align*}
&MAE = \dfrac{1}{n}\sum_{i=1}^n|\hat{Y}(s_i)-Y(s_i)| \\
&MSE = \dfrac{1}{n}\sum_{i=1}^n(\hat{Y}(s_i)-Y(s_i))^2
\end{align*}
Risultano: \begin{center}
\begin{tabular}{|l|l|l|}
\hline
Modello &MAE    &MSE\\
\hline
sferico &4.440892e-15   &3.944305e-30 \\
esponenziale &2.398082e-14  &3.23433e-29 \\
gaussiano   &1.154632e-14   &1.025519e-29\\
\hline
\end{tabular}\end{center}
I valori in entrambi i modelli non sono significativi, questo può essere dovuto alla dimensione campionaria ridotta. Sebbene fra loro i modelli per il variogramma non differiscono in modo significativo, viene scelto un modello sferico.
\subsection{Kriging ordinario}
Viene scelto il metodo del kriging ordinario perchè la media del processo gaussiano non è nota ma si ipotizza costante. E non si inseriscono altre covariate nel modello, in caso contrario sarebbe stato applicato un metodo di kriging universale
\begin{center}
\begin{figure}[h]
\includegraphics[width=\textwidth]{plot3.jpeg}
\caption{localizzazioni delle previsioni}\label{fig:2}
\end{figure}
\end{center}
I punti in cui sono state rilevate le concentrazioni di zinco sono molto sparsi rispetto alle previsioni che si vogliono fare per coprire l'intera area d'interesse (\emph{figura 3}).
\\\\
Secondo l'utilizzo del kriging ordinario la previsione per un punto $s$ non osservato è combinazione lineare dei dati osservati nei punti $s_i$: \begin{align*}
\hat{Y}(s) = \sum_{i=1}^n\lambda_i(s)Y(s_i)
\end{align*}
Viene di seguito rappresentata la concentrazione osservata e quella totale prevista in tutta l'area d'interesse \begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{plot4.jpeg}
    \caption{Previsione della concentrazione di zinco}
    \label{fig:4}
\end{figure}
\secntion{Conclusioni}
La previsione tramite il kriging ordinario sembra essere buona. Uno sviluppo ulteriori potrebbe essere quello di applicare il metodo del kriging universale inserendo delle covariate del modello. Infatti misure come la distanza dal fiume e l'elevazione rispetto al letto del fiume possono avere un'influenza significativa nel determinare la concentrazione di zinco in diverse locazioni.
\end{document}