\documentclass{article}
\begin{document}

\'A partir d'\'echantillons relev\'es sur le terrain, on cherche
\`a d\'eterminer la concentration par exemple d'un m\'etal dans
toute la zone. On utilise pour cela des mod\`eles de distribution
de la quantit\'e int\'eressante.

Pour fixer les id\'ees, notons $z(s)$ la quantit\'e de m\'etal (par
exemple en parties par million) en fonction de la position $s$
(en dimension 2 ou 3). Les mod\`eles \'etudi\'es ici approchent
$z(s)$ par $Z(s)$ avec~:
\[ Z(s) = m(s) + e(s) \]
o\`u $m(s)$ est la tendance ({\em trend} en anglais)
qui est d\'eterministe (par exemple $m(s)$ peut
\^etre constant ou \'egal \`a la distance \`a un fleuve, etc.)
et $e(s)$ repr\'esente la partie al\'eatoire de la distribution.

On dispose de $n$ \'echantillons, donc de $n$ mesures $z(s_i)$ 
($1\leq i \leq n$)
et il s'agit de trouver $m(s)$ et $e(s)$ qui permettent de pr\'evoir
$z(s)$ le plus pr\'ecis\'ement possible (krigeage)


\section{Interpolation} \label{sec:inv_distance}
Une m\'ethode d'estimation simple consiste \`a interpoler $z(s)$
par les $z(s_i)$ en mettant un poids proportionnel \`a l'inverse
de la distance $|s-s_i|$ au carr\'e~:
\[ z(s)= \frac{1}{C} \sum _{i=1}^n \frac{1}{|s-s_i|^2} z(s_i),
\quad C=\sum _{i=1}^n \frac{1}{|s-s_i|^2} \]

\section{Approcher la tendance par des mod\`eles lin\'eaires au sens des moindres carr\'es.} \label{sec:moindres_carres}
On se donne $p$ fonctions $f_j(s)$ et on cherche des coefficients 
$\beta _j$ tels que~:
\[ m(s)= \sum _{j=1}^p \beta _j f_j(s) \]
approche $z(s)$ le mieux possible, ``mieux'' \'etant d\'efini au sens
des moindres carr\'es, on veut minimiser $(z(s)-m(s))^2$ aux points
d'\'echantillonage. Le plus \'equilibr\'e consiste donc \`a minimiser~:
\[ g(\beta _1,...,\beta _p)=\sum _{i=1}^n (z(s_i)-m(s_i))^2 \]
Pour trouver les $\beta _j$, on d\'erive $g$ par rapport \`a 
$\beta _1,...,\beta _p$, on doit trouver 0 puisque la fonction poss\`ede
un extr\^emum en ce point.

Exemple si $p=1$~:
\[ g(\beta _1)=\sum _{i=1}^n (z(s_i)-\beta _1 f_1(s_i))^2 \]
On d\'erive par rapport \`a $\beta _1 $~:
\[ g'(\beta _1)=\sum _{i=1}^n -2 f_1(s_i) (z(s_i)-\beta _1 f_1(s_i)) \]
et on cherche $\beta _1$ tel que $g'(\beta _1)=0$~:
\[ \beta _1 = \frac{\sum _{i=1}^n f_1(s_i) z(s_i) }{\sum _{i=1}^n f_1(s_i)^2} \] 
Si $p=2$, on a~:
\[ g(\beta _1,\beta _2)=
\sum _{i=1}^n (z(s_i)-\beta _1 f_1(s_i)-\beta _2 f_2(s_i))^2  \]
et en d\'erivant par rapport \`a $\beta _1$ et $\beta _2$,
on obtient un syst\`eme lin\'eaire de deux \'equations 
\`a deux inconnues~:
\[ \left\{ \begin{array}{rcl}
 \sum _{i=1}^n -2 f_1(s_i) (z(s_i)-\beta _1 f_1(s_i)-\beta _2 f_2(s_i)) 
  & = & 0  \\
 \sum _{i=1}^n -2 f_2(s_i) (z(s_i)-\beta _1 f_1(s_i)-\beta _2 f_2(s_i)) 
 & = & 0
\end{array} \right. \]
donc~:
\[ \left(\begin{array}{cc}
 \sum_i f_1(s_i)^2 & \sum_i f_1(s_i) f_2(s_i) \\
 \sum _i f_1(s_i) f_2(s_i) & \sum _i f_2(s_i)^2
\end{array}\right) \left(\begin{array}{c}
 \beta _1 \\
  \beta _2
\end{array}\right)=\left(\begin{array}{c}
 \sum _i f_1(s_i) z(s_i) \\
 \sum _i f_2(s_i) z(s_i)
\end{array}\right) \]

Dans le cas g\'en\'eral, on a int\'er\^et \`a utiliser
les notations de l'alg\'ebre lin\'eaire (vecteurs et matrices). On note
$F$ la matrice \`a $n$ lignes et $p$ colonnes dont la $j$-i\`eme colonne
est form\'ee des valeurs de $f_j$ aux points $s_1$, ..., $s_n$ et on
note $F^*$ sa transpos\'ee (c'est-\`a-dire la matrice
obtenue par sym\'etrie par rapport \`a la diagonale). On note
$z$ le vecteur colonne form\'e par la valeur de $z(s)$ en $s_1$,..., $s_n$
et $\beta $ le vecteur colonne form\'e par $\beta _1$,..., $\beta _p$.
On montre alors que $\beta $ v\'erifie l'\'equation~:
\[ F^* F \beta = F^* z \]
Pour la r\'esoudre il faut donc inverser la matrice carr\'ee $F^* F$.

Remarques~:
\begin{itemize}
\item 
On peut bien sur se limiter \`a faire des calculs avec les \'echantillons
$s_i$ de la strate dans laquelle on veut faire des estimations si
on connait l'existence de strates.
\item 
Dans certaines situations, on met des poids diff\'erents $w_i$ aux
diff\'erents sites d'observations et on minimise plutot~:
\[ \sum _{i=1}^n w_i (z(s_i)-m(s_i))^2 \]
\end{itemize}

\section{Mod\'eliser les erreurs} \label{sec:variogrammes}
Dans une approche naive, on peut supposer que la partie al\'eatoire
$e(s)$ en un lieu est ind\'ependante de $e(s')$ en un autre lieu,
et que la distribution de ces variables al\'eatoires est identique,
on parle alors de mod\`ele IID ({\em independent identically distrbuted}).

On peut aussi consid\'erer que les parties al\'eatoires sont corr\'el\'ees
entre elles et ce d'autant plus que la distance entre $s$ et $s'$ est
petite. Dans les mod\`eles g\'eostatistiques, on consid\`ere
en fait que la corr\'elation ne d\'epend que de la distance entre
$s$ et $s'$ et on estime cette corr\'elation \`a l'aide des
\'echantillons $s_i$ en calculant ou en mod\'elisant le
{\em variogramme\/} (ou le covariogramme) des $s_i$. Une fois
ce calcul fait, on s'en sert pour pond\'erer le calcul de tendances au sens
des moindres carr\'es.

Covariogramme:\\
On veut estimer la covariance des variables al\'eatoires $e(s)$ et $e(s')$.
On suppose donc qu'elle ne d\'epend que de $h=|s-s'|$, la distance
de $s$ \`a $s'$, on la note Cov$(h)$.
On fait alors la moyenne de $e(s_i) e(s_j)$ pour tous les couples
$(i,j)$ de l'\'echantillonage tels que $|s_i-s_j| \in [h,h+\delta [$
o\`u $\delta $ est une constante, choisie de mani\`ere judicieuse~:
$\delta $ doit \^etre petit pour que la fonction covariogramme varie peu
sur l'intervalle $[h,h+\delta [$, mais il faut aussi que
cet intervalle contienne suffisamment de couples
de l'\'echantillon pour que la moyenne soit repr\'esentative.

Une autre fonction, appel\'ee variogramme, est souvent utilis\'ee
pour repr\'esenter la corr\'elation de $e(s)$ et $e(s')$, il s'agit
de la moyenne de $\frac{1}{2} (e(s_i)-e(s_j))^2$ pour 
$|s_i-s_j| \in [h,h+\delta [$.

Il existe diff\'erents mod\`eles de fonctions variogramme et covariogramme:
le mod\`ele de la p\'epite par exemple qui d\'ecorr\`ele un \'echantillon
du reste, des mod\`eles continus (sph\'eriques par exemple).
On peut \`a partir du variogramme de l'\'echantillon d\'eterminer
les coefficients d'un mod\`ele de variogramme approchant au mieux
ce vraiogramme, par exemple au sens des moindres
carr\'es, cf. l'appendice A.1 de la documentation du logiciel \verb|gstat|
disponible sur~:\\
\verb|http://www.gstat.org|

Lorsqu'on prend comme mod\`ele des erreurs variables al\'eatoires
$e(s)$ telles que~:
\[ \forall s, E(e(s))=0 , \quad
\mbox{Cov}(e(s),e(s'))= \mbox{Cov}(h) \]
($E(e(s))$ d\'esigne l'esp\'erance, c'est-\`a-dire la moyenne de $e(s)$),
on adapte la m\'ethode d'estimation par les moindres carr\'es \`a poids.
On note $V$ la matrice de coefficients~:
\[ v_{ij}=\mbox{Cov}(e(s_i),e(s_j))=\mbox{Cov}(|s_i-s_j|) \]
et $z(s)$, $m(s)$ les vecteurs~:
\[ z(s)=(z(s_1),...,z(s_n)), \quad m(s)=(m(s_1),...,m(s_n)) \]
On minimise alors le produit scalaire~:
\[ <z(s)-m(s)|V^{-1}(z(s)-m(s))> \]
par rapport aux coefficients $\beta _j$ pour d\'eterminer ces coefficients
$\beta _j$ (et donc $m$).
La partie erreur en un lieu $s_0$ est estim\'ee par~:
\[ \mbox{Cov}(e(s_0),e(s)) V^{-1} (Z(s)-m(s)) \]
Les calculs sont un peu plus
compliqu\'es, la valeur optimale de $\beta $ est alors~:
\[ \beta =(F^* V^{-1} F)^{-1} F^* V^{-1} z(s) \]
et on peut aussi calculer la variance de l'erreur $e(s_0)$, 
cf. la documentation logiciel {\tt gstat} (appendice A.2).

\section{Programmation.} \label{sec:data}
Il existe des formats standard pour repr\'esenter les donn\'ees des 
\'echantillons. Par exemple le format \verb|.eas| utilis\'e par \verb|gstat|,
cf. sur \verb|pcm2| le fichier \verb|/usr/local/gstat/cmd/zinc.eas|~:
\begin{verbatim}
Zinc measurements on River Meuse flood plains
3
xcoord, m
ycoord, m
zinc, ppm
181072 333611 1022
181025 333558 1141
181165 333537 640
...
\end{verbatim}
La premi\`ere ligne est un commentaire, la seconde indique le nombre
de donn\'ee de chaque \'echantillon (3 ici: 2 coordonn\'ees $x$, $y$ en
m\`etres, et la valeur de $z$ en ppm, ce qu'expliquent les 3 lignes suivantes),
puis chaque ligne correspond aux donn\'ees d'un \'echantillon.

On lit les donn\'ees dans le fichier et on les stocke dans un vecteur.
L'\'el\'ement de base du vecteur sera compos\'e de $s$ et de $z(s)$,
par exemple ici $s=(x,y)$~:\\
\verb|struct echantillon { double x; double y; double z; } ;|\\
(\verb|x| et \verb|y| sont les coordonn\'ees, \verb|z| repr\'esente $z(s)$).
Pour lire les donn\'ees \`a partir du fichier ci-dessus, par exemple~:
\begin{verbatim}
  ...
  ifstream infile("/usr/local/gstat/cmd/zinc.eas");
  // on saute les 5 premieres lignes
  char s[256];
  for (int i=0;i<5;i++)
    infile.getline(s,255);
  // on lit les donnees ligne par ligne
  vector<echantillon> v;
  echantillon e;
  for (;;){
    infile >> e.x >> e.y >> e.z;
    if (!infile) // si le fichier est termine on arrete la boucle
      break;
    v.push_back(e);
  }
  ... 
\end{verbatim}
Les calculs pr\'esent\'es ne posent pas de difficult\'es
algorithmiques sauf le calcul de l'inverse d'une matrice pour les
mod\`eles les plus sophistiqu\'es.

Calcul d'inverse de matrice~: soit on dispose d'une librairie externe,
soit on programme l'algorithme du pivot de Gau\ss(rappel de math\'ematiques: 
on colle la matrice identit\'e \`a la matrice \`a inverser, on
r\'eduit, lorsqu'on a la matrice identit\'e \`a gauche, alors la
matrice \`a droite est l'inverse de la matrice de d\'epart). Pour
le choix du pivot \`a l'int\'erieur d'une colonne donn\'ee, on
choisit le nombre de valeur absolue maximale afin de diminuer
les erreurs d'arrondi.

Pour le calcul du covariogramme, on pourra cr\'eer un type~:\\
\verb|struct couple { echantillon un; echantillon deux; double longueur;};|\\
On cr\'ee un vecteur \verb|v| de $n^2$ couples (\verb|vector<couple> v;|) 
contenant tous les couples d'\verb|echantillon| (en calculant la longueur)
et on ordonne ce vecteur par ordre croissant en utilisant la fonction
de tri \verb|sort| (ajoutez \verb|#include <algorithm>| pour pouvoir
l'utiliser). On d\'efinit la fonction de tri \verb|est_plus_pres|~:\\
\verb|bool est_plus_pres(couple a,couple b){|\\
\verb|  return a.longueur<b.longueur;|\\
\verb|}|\\
On trie le vecteur par~:\\
\verb|sort(v.begin(),v.end(),est_plus_pres);|\\
Il est ensuite ais\'e de d\'eterminer pour des valeurs de $h$ et $\delta $
donn\'ees les couples dont la \verb|longueur| appartient \`a $[h,h+\delta [$.

Remarque: la m\'ethode la plus efficace pour trouver le premier couple 
dont la longueur est sup\'erieure \`a $h$ consiste \`a couper en deux 
moiti\'es \'egales le vecteur de couples. Si le couple du milieu a
une longueur sup\'erieure \`a $h$, on recommence avec la premi\`ere
moiti\'e du vecteur, sinon avec la deuxi\`eme moiti\'e, on s'arr\^ete
lorsqu'il n'y a plus qu'un couple dans le vecteur.

\end{document}