RIA

laboria

ARCHIVES

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 5 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

<figure><img src="image_1.png" /></figure>

LISSAGE TYPOLOGIQUE

E. DIDAY
Y. OK
M. JACOTTET
R. TOMASSONE

Rapport de Recherche N° 164

Juin 1976
# LISSAGE TYPOLOGIQUE

Y. Ok, E. Diday

IRIA

M. Jacottet*, R. Tomassone**

## Résumé :

Etant donné un ensemble d'objets caractérisés par un paramètre continu (le temps par exemple) avec des données manquantes, le problème est de trouver simultanément une partition des objets et des fonctions polynominales telles que l'adéquation entre chaque classe de la partition et la fonction polynominale correspondante, soit la meilleure possible au sens des moindres carrés. Un algorithme efficace est proposé, ses propriétés de convergence sont étudiées, l'aspect statistique est considéré. Le programme est illustré par une application en agro-biologie.

## Abstract :

Given a set of objects characterized by a continuous parameter (time for instance) with missing data, the problem is to find simultaneously a partition of the objects and a set of polynomial functions such that the fit between each class of the partition and the corresponding polynomial function is optimal in the sense of the mean square criterion. An efficient algorithm is given, its convergence properties are studied, the statistical aspect is considered. The program is illustrated by an application in agro-biology.

<footnote>* Ecole Polytechnique Fédérale de Lausanne - Département de Mathématiques, 61, avenue de Cour - 1007 Lausanne.
I.N.R.A. - Laboratoire de Biométrie du C.N.R. 2 - Jouy en Josas.</footnote>

<footnote>** I.N.R.A. - Laboratoire de Biométrie du C.N.R. 2 - Jouy en Josas.</footnote>
LISSAGE TYPOLOGIQUE

Une méthode de traitement des séries chronologiques incomplètes
# SOMMAIRE

## INTRODUCTION

1. - NOTATIONS UTILISEES

2. - LE FORMALISME DE LA METHODE

3. - L'ALGORITHMEE ET SA CONVERGENCE

4. - RECHERCHE DU POLYNOME DE REGRESSION D'UNE CLASSE k : interprétations statistique et algébrique de la fonction g

4.1. - Hypothèse sur la variance et mesures de proximité

4.2. - Centre de gravité d'une classe

4.3. - Recherche des estimateurs au sens des moindres carrés des coefficients du polynôme de régression et aspect algébrique

4.3.1. - Estimateurs au sens des moindres carrés, pour la norme $V^{-1}$ des coefficients du polynôme de régression du centre de gravité $G_k^J$ de la classe k

4.3.2. - Aspect algébrique de la recherche du polynôme du centre de gravité d'une classe

4.3.3. - Proposition 2

4.4. - Aspect statistique

4.4.1. - Les modèles

4.4.2. - Tests sous l'hypothèse de normalité

4.5. - Lien entre l'aspect statistique et l'aspect algébrique

5. - LA BASE D'INTERPOLATION

5.1. - Choix de la base

5.2. - Construction de la base de polynômes orthogonaux sur le support $J_k$

6. - QUALITE DE LA REGRESSION. TEST D'AJUSTEMENT DU DEGRE DU POLYNOME

6.1. - Le problème

6.2. - Procédure utilisée

7. - LE PROGRAMME

8. - APPLICATIONS

9. - CONCLUSION

<figure><img src="image_2.png" /></figure>
# INTRODUCTION

A l'origine de ce travail se trouvent des problèmes concrets, d'une part le traitement de données décrivant l'évolution d'un phénomène susceptible d'être représenté par une courbe, c'est-à-dire des données pour lesquelles un modèle polynômial convient, d'autre part le problème des données incomplètes car il est fréquent que les mesures n'aient pas été faites dans les mêmes conditions pour tous les individus.

De telles données sont, par exemple, des courbes de croissance pondérale d'animaux en fonction du temps ou des courbes de teneur en minerai en fonction de la profondeur ...etc.

Il s'agit donc d'un problème très précis : la recherche dans une population donnée de groupes d'individus ayant le même type de courbe d'évolution. On veut détecter dans un échantillon multidimensionnel donné, l'existence possible de sous-ensembles susceptibles d'être caractérisés par une fonction polynômiale.

La solution présentée ici constitue une approche orientée vers la classification automatique du traitement des chroniques, plus spécialement des chroniques incomplètes.

L'algorithme proposé "Lissage Typologique" se propose d'atteindre ce but en réalisant en quelque sorte une synthèse ou tout au moins en utilisant simultanément la méthode de classification automatique des nuées dynamiques et une méthode de lissage polynômial.

C'est une technique de la famille des nuées dynamiques (M.N.D., on emploie l'abréviation plus loin sans l'expliciter), i.e. une technique qui, au moyen d'un processus itératif convergent, améliore au sens d'un critère donné une partition accessible à partir d'une partition quelconque aléatoire ou estimée.

Dans le cas particulier de données incomplètes, au lieu de faire précéder la classification proprement dite de procédures d'approximation visant à compléter les données en estimant au mieux les mesures manquantes, l'algorithme LISTYP effectue simultanément classification et approximation.

Ceci permet d'éviter plusieurs inconvénients. Tout d'abord, le procédé qui consisterait à effectuer autant d'approximation que d'individus serait très coûteux, d'autre part, pour estimer les mesures manquantes d'un individu, il est
souhaitable de tenir compte des autres individus, en tout cas de ceux qui lui
ressemblent le plus, c'est-à-dire de substituer à la notion classique d'interpo-
lation, une notion d'_interpolation statistique_ (i.e. par groupes d'individus)
puisqu'en effet on s'intéresse davantage aux traits caractéristiques et distinc-
tifs du groupe qu'aux individus.

D'un point de vue pratique, on a considéré que les lignes du tableau des données étaient les valeurs connues sur un certain support de fonctions d'une variable définies par un procédé expérimental.

On s'est limité au cas où le support J est discret, naturellement ordonné
et à valeurs régulièrement espacées.

Il ne s'agit donc nullement, du moins dans l'état actuel des travaux, d'une
méthode générale de traitement des données incomplètes, mais uniquement du traï-
tement des séries chronologiques incomplètes pour lesquelles un modèle polynômial
peut être retenu.

Comme, généralement, on ne connaît pas les fonctions qui ajusteraient au
mieux les données, il a semblé raisonnable de chercher des polynômes de bas
degré pour caractériser les classes.

Les formules choisies pour mesurer les proximités entre individus et poly-
nômes (D) et entre classes et polynômes (R) conduisent à appliquer la méthode d'appro-
ximation au sens des moindres carrés à chaque classe et à la recherche de son
polynôme.

comme
$$
\sum_{i \in P_k} \sum_{j \in J} \mu_j (x_i^j - \mathcal{L}^j)^2 = R(\mathcal{L}, k, P), \text{ on est amené à minimiser}
$$

$$
\sum_{j \in J} (\mu_j / |J|) (g_k^j - \mathcal{L}^j)^2 \text{ où } g_k^J \text{ est le centre de gravité de la}
$$

classe P_k i.e. à rechercher le polynôme de meilleure approximation au sens des
moindres carrés du centre de gravité. Les formules ci-dessus correspondent à des
données complètes.

Dans la situation où les données sont incomplètes, on montre qu'on peut
également, à condiction de postuler convenablement la forme des erreurs, se
ramener à ce cas. L'algorithme LISTYP a pour résultat une partition de la popu-
lation en classes caractérisées, chacune par une fonction polynômiale dont le
degré peut être ajusté grâce à des tests de validité. La contribution de la classe à la valeur du critère mesure sa dispersion autour de la courbe qui la caractérise et permet d'apprécier la qualité de cette caractérisation.

On peut donc aboutir à la reconnaissance de sous-ensembles qui seraient des échantillons correspondant à différents types de courbes d'évolution du phénomène étudié.

# 1. - LES NOTATIONS UTILISEES

La présence de "trous", i.e. de données manquantes, dans le tableau à traiter permet difficilement de concilier les exigences de clarté et de concision auxquelles doivent satisfaire de "bonnes" notations. Pour être suffisamment explicites, les notations adoptées, sont relativement lourdes.

Soient $I$ et $J$ deux ensembles finis d'indices.

$i \in I$ désigne un individu, $j \in J$ un paramètre.

$x_i^j$ est le résultat de la mesure du paramètre $j$ sur l'individu $i$.

Vi $\in I$ on note $J_i$ l'ensemble des paramètres effectivement mesurés sur l'individu $i$.

$J_i = \{j \in J | x_i^j \text{ non manquant}\}; |J_i| = \text{card } J_i$

$V_j \in J$ on note $I_j$ l'ensemble des individus pour lesquels on a effectué la mesure du paramètre $j$.

$I_j = \{i \in I | x_i^j \text{ non manquant}\}; |I_j| = \text{card } I_j$

Par suite tout individu $i$ peut être représenté par le vecteur de $\mathbb{R}$ noté $x_i^J = \{x_i^j | j \in J_i\}$ ; tout paramètre par un vecteur de $\mathbb{R}$ noté $x_i^j = \{x_i^j | i \in I_j\}$

L'ensemble des individus peut être représenté par :

$$E = \{x_i^J | i \in I\}$$

$E$ correspond au tableau noté $x_I^J$ comportant des "trous" et $\sum_{i \in I} |J_i|$ données non manquantes.
Soit $P = (P_1, \dots, P_k, \dots, P_K)$ une partition de $I$ en $K$ classes.

Considérant la classe $P_k$ on définit :

- l'ensemble $J_k$ des paramètres mesurés au moins une fois sur la classe k

$$J_k = \{j \in J \mid \exists i \in P_k \text{ où } x_i^j \text{ est non manquant}\} = \bigcup_{i \in P_k} J_i$$

- l'ensemble $I_{j,k}$ des individus de la classe $P_k$ dont on a mesuré le paramètre j

$$I_{j,k} = \{i \in P_k \mid x_i^j \text{ non manquant}\}$$

## 2. - LE FORMALISME DE LA METHODE

Pour reprendre les notations habituelles aux M.N.D. l'espace des noyaux $\mathcal{L}$ est ici l'ensemble des polynômes de degré inférieur ou égal à n.

Pour $\mathcal{L} \in \mathcal{L}$ on note $\mathcal{L}^J$ l'ensemble $\{\mathcal{L}(j) \mid j \in J\}$ des valeurs prises par le polynôme sur le support fini d'approximation $J$; $\forall j, \mathcal{L}(j)$ sera noté $\mathcal{L}^j$.

On sait que $\mathcal{L}$ est un espace vectoriel de dimension $n+1$. Tout polynôme de $\mathcal{L}$ peut donc être décomposé suivant une base $\Lambda = \{\Lambda_0, \Lambda_1, \dots, \Lambda_n\}$ en $\mathcal{L} = \sum_{\alpha=0}^{n} C_{\alpha} \Lambda_{\alpha}$ et caractérisé par le vecteur des coefficients relatifs à cette base. $\mathcal{L}$ peut donc être identifié au vecteur $C = (C_0, C_1, \dots, C_n)$.

On définit :

a) une mesure de proximité $D$ entre individus et polynômes par :

$$D : E \times \mathcal{L} \to \mathbb{R}^+$$

$$D(x_i^J, \mathcal{L}) = \frac{1}{\sum_{j \in J} \mu_j} \left( \sum_{j \in J_i} \mu_j (x_i^j - \mathcal{L}^j)^2 \right) \quad (2.1)$$

où $\mu_j$ est une pondération de l'ensemble des paramètres

(par exemple : $\mu_j = 1$; $\mu_j = \frac{1}{|J|}$; $\mu_j = \frac{1}{|J|} \cdot \frac{|I_j|}{|I|}$ etc ...)
b) une mesure de proximité R entre classes et polynômes par :

$$
R : \mathcal{L} \times \mathcal{K} \times \mathcal{P}_{\mathcal{K}} \to \mathbb{R}^{+} \qquad (2.2)
$$

$$
R (\mathcal{L}, k, P) = \sum_{i \in P_k} D(x_i^J, \mathcal{L})
$$

c) la fonction f de la même manière que dans toutes les variantes M.N.D. par :

$$
f : \mathcal{L}_K \to \mathcal{P}_K \quad \forall \mathcal{L} = (\mathcal{L}_1, \dots, \mathcal{L}_K), f(\mathcal{L}) = P \text{ où}
$$

$$
P_k = \{i \mid i \in I; D(x_i^J, \mathcal{L}_k) = \inf \{D(x_i^J, \mathcal{L}_{k'}) : k' \in K\}\}^*
$$

d) la "fonction" g par :

$$
g : \mathbb{P}_K \to \mathbb{L}_K \text{ est}
$$

g (P) = L = (L_1, ..., L_K, ...) où L_K est le polynôme de degré inférieur
où égal à n qui minimise R(., k, P) on peut montrer que pour une base fixée, il
existe et est unique (cf. 4.4 et 4.5) et que c'est un polynôme de régression
(cf. 4.2) nous l'appellerons n-polynôme de régression de la classe P_k.

e) le critère W mesure de proximité entre K-noyaux et K-partitions par :

$$
W : \mathcal{L}_K \times \mathcal{P}_K \to \mathbb{R}^+
$$

$$
\begin{align*}
W (L, P) &= \sum_{k \in \mathcal{K}]} \sum_{i \in P_k} D(x_i^J, \mathcal{L}_k) \\
&= \sum_{k \in \mathcal{K}]} \sum_{j \in J} \sum_{i \in I_{j,k}} \frac{(x_i^j - \mathcal{L}_k^j)^2}{\sum_{j' \in J} (j - j_{i,k})^2}
\end{align*}
$$

3. - L'ALGORITHME ET SA CONVERGENCE

Le processus d'itérations utilisé est décrit dans l'organigramme suivant.

<footnote>* le cas suivant : $\exists k, \exists k'$, $D(x_i^J, \mathcal{L}_k) = D(x_i^J, \mathcal{L}_{k'})$ est traité comme en 2.1.2.</footnote>
3.1. - Organigramme

<figure><img src="image_3.png" /></figure>
L'étape 1 est le départ aléatoire de l'algorithme, i.e. une partition quel-
conque en K classes; K est fixé arbitrairement ou d'après la connaissance que
l'on peut avoir des données.

L'ensemble des étapes 2, 3, 4, 5 forme une itération; l'étape 2 (en parti-
culier 2.1 et 2.2) constitue l'originalité de cet algorithme par rapport aux
autres méthodes de la famille des Nuées Dynamiques, son contenu sera précisé au
paragraphe 4.

Le test de convergence et l'arrêt de l'algorithme sont justifiés par la pro-
priété de convergence démontrée ci-dessous (cf. 3.2).

Les résultats obtenus pour chaque tirage différent (étape 6) consistent en
une partition de l'ensemble I en K classes (avec éventuellement une ou plusieurs
classes vides) que l'on caractérise chacune par un polynôme de degré inférieur
ou égal à un degré n fixé a priori et estimé grâce à un test d'ajustement
(voir § 6).

3.2. - <u>Convergence</u>

Au calcul du critère effectué à l'étape 3, pour chaque itération,
correspond la suite $u_n = W (L^{(n)}, P^{(n)})$.

On peut montrer que sous réserve de l'_existence_ et de l'_unicité_ du
n-polynôme de régression de la classe P<sub>k</sub>, ∀<sub>k</sub> (démontrées en 4.3), les défini-
tions choisies pour f, g et R entraînent la décroissance de la suite u<sub>n</sub> et par
suite la convergence de l'algorithme.

La démonstration est classique dans la M.N.D.

4. - <ins>RECHERCHE DU POLYNOME DE REGRESSION D'UNE CLASSE k : INTERPRETATIONS</ins>
<ins>STATISTIQUE ET ALGEBRIQUE DE LA FONCTION g</ins>

L'objet de ce paragraphe est de montrer que la formule choisie pour l'écart
quadratique, correspond à la méthode usuelle des moindres carrés appliquée à
une classe P<sub>k</sub> et à la recherche de son polynôme à condition de postuler conve-
nablement la forme des erreurs.

<ins>RAPPEL</ins> : On appelle n-polynôme de régression de la classe P<sub>k</sub>, le polynôme de degré inférieur ou égal à n qui minimise R(., k, P) (cf. 2 d).
### 4.1. - Hypothèse sur la variance et mesures de proximité

Considérons un individu i représenté par le vecteur

$x_i^J$ des $|J_i|$ observations faites sur lui.

Considérons un polynôme $\mathcal{L} = \Lambda C$ et le vecteur correspondant $\mathcal{L}_i^J$ des $|J_i|$ valeurs de ce polynôme sur le support $J_i$.

Munir $\mathbb{R}^{|J_i|}$ d'une métrique $(V_i)^{-1}$ définissant une distance entre éléments de $\mathbb{R}^{|J_i|}$ par :

$$ 
\begin{aligned}
d^2(x_i^J, \mathcal{L}_i^J) &= \left[ x_i^J - (\Lambda C)^{J_i} \right] \left[ V_i^{-1} \right] \left[ x_i^J - (\Lambda C)^{J_i} \right] \\
&= \left\| x_i^J - (\Lambda C)^{J_i} \right\|_{V_i^{-1}}^2
\end{aligned}
$$

et chercher les coefficients $C$ minimisant cette distance équivaut à poser pour l'individu i le modèle statistique :

$$ 
\begin{aligned}
\left| E(x_i^J) - (\Lambda C) \right. \\
\left. \mathrm{var}(x_i^J) = \sigma^2 V_i \right.
\end{aligned}
$$

et réciproquement.

On remarque alors que les mesures de proximité définies en 2 s'interprètent de la façon suivante :

$$ (1) \qquad D(x_i^J, \mathcal{L}_i^J) = \left\| x_i^J - \mathcal{L}_i^J \right\|_{V_i^{-1}}^2 $$

$$ (2) \qquad R(\mathcal{L}, k, P) = \sum_{i \in P_k} \left\| x_i^J - \mathcal{L}_k^J \right\|_{V_i^{-1}}^2 $$

$$ (3) \qquad W(L, P) = \sum_{k \in \mathbb{K}} \sum_{i \in P_k} \left\| x_i^J - \mathcal{L}_k^J \right\|_{V_i^{-1}}^2 $$

$$ \text{avec } V_i^{-1} = \frac{1}{\sum_{j \in J_i} \mu_j} \left[ \cdot \mu_j \cdot \right]_{j \in J_i} $$
Remarques :

1 - La définition choisie pour ces fonctions nous impose l'hypothèse faite sur la variance dans le modèle ci-dessus.

2 - Un polynôme $\mathcal{L}$ minimisant (2) est bien un polynôme de régression. Le problème est de trouver $\mathcal{L}$ - i.e. Les coefficients $C_0, \dots, C_n = C$.

3 - Si $V_j \mu_j = 1$, $\Sigma\{\mu_j | j \in J_i\} = |J_i|$ on prend ce cas en 4.3.

4.2. - Centre de gravité d'une classe

Définition

On appelle centre de gravité de la classe k noté $G_k^{J(*)}$ le vecteur à
$|J_k|$ composantes définies par :

$$
\forall j \in J_k, G_k^j = \sum_{i \in I_{j,k}} \frac{1}{|J_i|} \frac{x_i^j}{\sum_{i \in I_{j,k}} \frac{1}{|J_i|}}
$$

ou encore en posant $v_{j,k} = \sum_{i \in I_{j,k}} \frac{1}{|J_i|}$ où $v_{j,k}$ est la masse totale

des éléments de P_k dont on a mesuré le paramètre j

$$
\forall j \in J_k, G_k^j = \frac{1}{v_{j,k}} \sum_{i \in I_{j,k}} \frac{x_i^j}{|J_i|}
$$

Notons $V(k) = \left[ \begin{array}{c|c} 1/\mu_{j_i} & v_{j_1,k} \\ \hline \vdots & \vdots \\ 1/\mu_{j_m} & v_{j_m,k} \end{array} \right]$

<footnote>(*) $G_k^J$ est le centre de gravité au sens usuel si l'on considère la classe k comme le nuage $\{x_i^{J,k} , \frac{1}{|J_i|} | i \in P_k\}$ où $x_i^{J,k}$ a pour poids $\frac{1}{|J_i|}$ comme $J_i \subset J_k$, si $j \in J_k$ et $j \notin J_i$, i.e. si $x_i^j$ est manquant, la définition adoptée revient à supposer $x_i^j = 0$.</footnote>
### 4.3. - Recherche des estimateurs au sens des moindres carrés des coefficients du polynôme de régression et aspect algébrique

#### 4.3.1. - Estimateurs au sens des moindres carrés, pour la norme $v^{-1}$ des coefficients du polynôme de régression du centre de gravité $G_k^J$ de la classe k.

Soit $\hat{C}$ un ensemble des estimateurs au sens des moindres carrés de $C$, c'est-à-dire un ensemble de fonctions de $G_k^J$ dont les valeurs minimisent

$$ ||G_k^J - \Lambda C||^2_{V^{-1}} \quad (\text{cf. Scheffé [2] définition page 9}). $$

$\hat{C}$ est donc la solution du système d'équations normales écrit matriciellement sous la forme $\Lambda' V^{-1} G - \Lambda' V^{-1} \Lambda = 0$ et obtenu en dérivant $||G - \Lambda C||^2_{V^{-1}} = \mathcal{G}(G, C)$

Si $\Lambda$ est de rang $n+1$, $S^{-1} = (\Lambda' V^{-1} \Lambda)^{-1}$ existe et la solution unique est donnée par : $\hat{C} = S^{-1} \Lambda' V^{-1} G$.

#### 4.3.2. - Aspect algébrique de la recherche du polynôme du centre de gravité d'une classe

Soit $n$ le degré maximum du polynôme cherché, soit $m = \text{card } J_k$ (m très grand relativement à n)

on considère les espaces :

$$ \mathbb{R}^m = F; \mathbb{R}^{n+1} = E $$

et leurs duaux algébriques :

$F^*$ et $E^*$

On suppose $\mathbb{R}^m$ muni d'une métrique $M$ soit $\Lambda = (\Lambda_0, \dots, \Lambda_n)$ une base de l'espace vectoriel des polynômes de degré inférieur ou égal à n (où $\Lambda_\alpha$ est de degré $\alpha$).

Soit $\Lambda(k) = (\Lambda_0^J, \dots, \Lambda_n^J)$ la matrice constituée par $(n+1)$ vecteurs de $\mathbb{R}^m$, où $\Lambda_\alpha^J$ est la valeur prise au point du support par le polynôme $\Lambda_\alpha$ de la base.
Soit $H(\Lambda)$ le sous-espace vectoriel de $F = \mathbb{R}^m$ engendré par
$(n+1)$ vecteurs ci-dessus, soit à estimer $G$.

Le problème posé, à savoir trouver $(C_0, \dots, C_n)$ tel que
$$ \sum_{\alpha=0}^{n} \Lambda_{\alpha} C_{\alpha} $$ soit le point le plus proche de $G$ au sens de la métrique $M$, équivaut à:

Trouver $\hat{G} \in H(\Lambda)$ tel que $d_M^2(G, \hat{G}) = \inf_{v \in H(\Lambda)} d_M^2(G, v)$

On sait que le point $\hat{G}$ est la projection M-orthogonale de $G$ sur $H(\Lambda)$ et que cette projection est unique.

On a donc le résultat suivant :

$$ (1) \quad \forall G \in \mathbb{R}^m \exists! \hat{G} \in H(\Lambda) \mid |G - \hat{G}|_M^2 = \inf_{v \in H(\Lambda)} \mid |G - v|_M^2 $$

**Remarque** : Bien que cet élément soit unique, sa décomposition selon $\Lambda(k)$ n'est unique que dans le cas où $\Lambda(k)$ est une base de $H(\Lambda)$ sinon $\hat{G}$ admet une infinité de décomposition suivant $\Lambda(k)$.

On peut montrer le lemme et la proposition suivants qui prouvent l'existence et l'unicité sous certaines hypothèses du polynôme de régression d'une classe.

**Lemme 1** : Si $\Lambda$ est une base de polynômes M-orthogonaux sur le support $J_k$, alors $\Lambda(k)$ est une base M-orthogonale de $H(\Lambda)$.

**Proposition 1** : Soit $G \in \mathbb{R}$, une base $\Lambda$ de polynômes M-orthogonaux sur le support $J_k$ étant choisie, il existe $\hat{C} = (\hat{C}_0, \dots, \hat{C}_n)$ unique et par suite il existe un polynôme $\hat{\mathcal{L}} = \Lambda \hat{C}$ unique tel que :

$$ ||G - \hat{\mathcal{L}}^J_k|^2_M = \inf_{\mathcal{L} \in L} ||G - \mathcal{L}^J_k||_M^2 $$

$(\hat{\mathcal{L}} \in L)$ est défini $\forall x \in \mathbb{R}$ par $\hat{\mathcal{L}}(x) = \sum_{\alpha=0}^{n} \hat{C}_{\alpha} \Lambda_{\alpha}(x)$ et tel que $\hat{G} = \hat{\mathcal{L}}^J_k$ (notations § 2).

## Aspect algébrique de l'estimation de G

$\Lambda(k)$ peut être considérée comme une application linéaire de rang $(n+1)$ de $E^*$ dans $F$, et $\Lambda' M \Lambda$ définit un isomorphisme de $E^*$ dans $E$.
On a donc le diagramme suivant :

<figure><img src="image_4.png" /></figure>

On note :

$$ (3) \qquad A_n = \Lambda(\Lambda' M \Lambda)^{-1} \Lambda' M $$

$A_n$ est une application de $\mathbb{R}^m$ dans $H(\Lambda)$ dont on peut montrer facilement qu'elle est idempotente et M-symétrique. C'est donc bien un opérateur de projection M-orthogonale.

Estimer $G$ par un polynôme de degré inférieur ou égal à $n$ équivaut à déterminer par (3) l'opérateur $A_n$, quel que soit $n$, pour toute classe $P_k$.

**Remarque** : On montre dans le paragraphe suivant une proposition qui nous permettra d'affirmer que pour un degré $n$ fixé, l'application $g$ de l'algorithme décrit au § 3 est une fonction qui s'écrit :

$$ g(P) = [A_n^1(G_1), \dots, A_n^i(G_i), \dots, A_n^k(G_k)] $$

$G_i, V_i^{-1}, \Lambda_i$ sont respectivement le centre de gravité, la métrique, la matrice $\Lambda$, associés à la classe $P_i$.

$A_n^i$ est le projecteur défini par (3) avec $M = V_i^{-1}$
$\Lambda = \Lambda_i$
4.3.3. - **Proposition 2**

<figure><img src="image_5.png" /></figure>

Le $n$-polynôme de régression de la classe $P_k$ (i.e. le polynôme minimisant $R(., k, P)$) est le polynôme de meilleure approximation au sens des moindres carrés, du centre de gravité de la classe $P_k$.

Démonstration :

$\hat{\mathcal{L}}$ étant estimé par $\hat{\mathcal{C}} = (\Lambda' M \Lambda)^{-1} \Lambda' M G$, soient
$\hat{\mathcal{L}} = \Lambda \hat{\mathcal{C}}$ le polynôme correspondant, $\hat{\mathcal{L}}_k = \Lambda(k) \hat{\mathcal{C}}(k) = (\hat{\mathcal{L}}^j_1, \dots, \hat{\mathcal{L}}^j_m)$,

les $|J_k|$ valeurs estimées du polynôme du centre de gravité sur le support $J_k$ et
$\mathcal{L}_k^* = \Lambda_{(k)}^* \cdot C_k$ les valeurs estimées sur le support $J$.

Par construction $\hat{\mathcal{L}}$ est tel que (1)

$$
(G - \hat{\mathcal{L}})^{-1} (G - \hat{\mathcal{L}}) = \sum_{j \in J_k} (G_k^j - \hat{\mathcal{L}}^j)^2 \mu_j v_{j,k} =
$$

$$
\begin{align*}
(E) &= \sum_{j \in J_k} \mu_j \left( \sum_{i \in I} \frac{x_i^j}{|J_i|} \right)^2 / v_{jk} + v_{jk} (\hat{\mathcal{L}}^j)^2 - 2 \sum_{i \in I} x_i^j \frac{\hat{\mathcal{L}}^j}{|J_i|} \\
&= \inf_{\mathcal{L} \in \mathcal{L}} \left| |G - \mathcal{L}| \right|^2 v^{-1}
\end{align*}
$$

L'expression (E) est une fonction de $\hat{\mathcal{L}}$ de la forme :

$$
F(\hat{\mathcal{L}}) = \sum_{j \in J_k} \left[ \alpha^j + f(\mathcal{L}^j) \right] \quad \text{où } \alpha^j \text{ est une constante correspondant au premier} \\
\text{terme entre crochets et où } f \text{ est une fonction définie par les deux derniers} \\
\text{termes. } F(\hat{\mathcal{L}}) \text{ peut encore s'écrire } F(\hat{\mathcal{L}}) = \alpha + \sum_{j \in J_k} f(\hat{\mathcal{L}}^j) = \alpha + h(\mathcal{L}).
$$

Le résultat démontré précédemment (Proposition 1) s'écrit
$$ \exists! \hat{\mathcal{L}} \forall \mathcal{L} \in \mathcal{L} F(\hat{\mathcal{L}}) = \inf_{\mathcal{L} \in \mathcal{L}} F(\mathcal{L}). $$

<footnote>(1) Pour ne pas alourdir les notations, on omet, dans la suite, l'indice k, mais on raisonne évidemment pour une classe donnée.</footnote>
Puisque $\alpha$ est une constante, on a donc :

$$h(\hat{\mathcal{L}}) = \sum_{j \in J_k} \mu_j \left[ \sum_{i \in I} \frac{(\hat{L}^j)^2 - 2 \sum_i x_i^j \hat{L}^j}{|J_i|} \right] = \inf_{\mathcal{L} \in \mathcal{L}} h(\mathcal{L})$$

et par suite comme : $\forall \mathcal{L}, \forall P$,

$$R(\mathcal{L}, k, P) = \sum_{j \in J_k} \mu_j \sum_{i \in I} \frac{1}{|J_i|} (x_i^j - \mathcal{L}^j)^2 \\ = \sum_{j \in J_k} \mu_j \sum_{i \in I} \frac{(x_i^j)^2}{|J_i|} + h(\mathcal{L})$$

Le polynôme $\hat{\mathcal{L}}$ minimisant $h(\mathcal{L})$ minimise $R(\mathcal{L}, k, P)$.

c.q.f.d.

### 4.4. - Aspect statistique

**Notations :**

La population E est un ensemble de T individus connus par des mesures $y_t(x_j)$, $t = 1, ..., T$ sur un paramètre x dont $x_j$ est le $j$-ième niveau.

$\forall \mathcal{P} \in \mathcal{L}$, ensemble des polynômes de degré d

$$D(y, \mathcal{P}) = \frac{1}{|J|} \sum_{j \in J(y)} q_j \left[ y(x_j) - \mathcal{P}(x_j) \right]^2$$

où $J(y) = \{j | y(x_j) \text{ existe}\}$ ; $|J| = \text{card } J(y)$ ; $q_j$ est un poids quelconque.

#### 4.4.1. - Les modèles

La fonction R à minimiser est donc :

$$R(\mathcal{P}, k, P) = \sum_{y \in C_k} \frac{1}{|J_k|} \cdot \sum_{j \in J(y)} q_j \left[ y(x_j) - \mathcal{P}(x_j) \right]^2 = \sum_{y \in C_k} \left\| y - \mathcal{P} \right\|^2 = \sum_{y \in C_k} D(y, \mathcal{P}).$$

Cette norme doit correspondre, pour se rattacher à l'aspect statistique de la recherche d'un polynôme de régression, à la norme de la covariance des variables y; c'est-à-dire qu'il faut faire l'hypothèse que :
$$ \mathrm{Var}(y(x_j)) = \frac{J}{q_j} \sigma_k^2 \quad y \in E, \ j \in J(y) \quad (4-4-1) $$

Les différentes mesures sont supposées indépendantes. L'hy-
pothèse sur la variance est raisonnable: plus $J_y$ est grand, c'est-à-dire que
plus il y a de mesures $y(x_j)$ pour la suite $y$, moins le polynôme de régression
est précis (à degré fixé).

$q_j$ peut être un facteur de proportionnalité avec le nombre de mesures effectuées
au point $x_j$.

$\sigma_k^2$ est un paramètre inconnu dépendant de la classe.

Nous pouvons donc maintenant poser le modèle statistique
suivant :

$$ \begin{gathered} \Omega : y_t(x_j) = \sum_{i=0}^{d} c_i P_i(x_j) + e_{tj}; (t : y_t \in C_k) \\ e_{tj} \text{ indépendantes et telles que } E(e_{tj}) = 0, \mathrm{Var} \ e_{tj} = \frac{J}{q_j} \sigma_k^2 \end{gathered} $$

(4-4-2)

Les $P_i$ sont des polynômes orthogonaux sur le support
$\{x_j, j \in J(C)\}$ de degré $i$ exactement (cf. § 4.4.1)

$$ J(C) = \bigcup_{y \in C} J(y). $$

On a donc :

$$ E(y_t(x_j)) = \sum_{i=0}^{d} c_i P_i(x_j) = \mathcal{P}(x_j). $$

Cette espérance ne dépend pas de $y_t$; on peut donc simpli-
fier le modèle en "sommant" sur $t$.

Réduction du modèle

Soit $J^{-1}(j)$ l'ensemble des $y \in C$ définis en $x_j$. (Il se peut
que $\exists j$ et $\exists C$, tels que $J^{-1}(j) = \emptyset$).
Affectons aux valeurs $y(x_j)$ la masse $1/J_y$ et considérons

la suite :

$$ \{G(x_j)\}, \quad j \in J(C) $$

des centres de gravité, correspondant au système de masse

$$ \{1/J_y\}, \quad y \in C $$

$$ G(x_j) = \frac{\sum_{y \in J} \frac{y(x_j)}{y_j}}{\sum_{j \in J} \frac{1}{y_j}}, \quad j \in J(C) \quad (4-4-3) $$

**Proposition** : Le polynôme de régression de degré d de la classe C est le même que le polynôme d'interpolation de la suite $\{G(x_j)\}$, $j \in J(C)$, pour la norme de la covariance, respectivement des $y_t(x_j)$ et des $G(x_j)$.

**Démonstration** :

$$ E(G(x_j)) = \mathcal{P}(x_j) $$

$$ \text{Posons } N_j = \sum_{y \in J} \frac{1}{y_j} $$

$$ \text{Var}(G(x_j)) = \frac{\sigma_k^2}{N_j q_j} $$

Comme les $y_t(x_j)$ sont indépendantes, les $G(x_j)$ le sont aussi.

Le modèle $\Omega$ se réduit donc au modèle :

$$ \begin{gathered} \hat{\Omega} : G(x_j) = \mathcal{P}(x_j) + \hat{e}_j ; \quad j \in J(C_k) \quad (4-4-4) \\ \hat{e}_j \text{ indépendantes, } E(\hat{e}_j) = 0, \text{ Var } \hat{e}_j = \frac{\sigma_k^2}{N_j q_j} \end{gathered} $$

$$ \text{Le polynôme } \mathcal{P}(x_j) = \sum_{i=0}^{d} c_i P_i(x_j) \text{ est le polynôme qui} $$

minimise :

$$ \begin{align} \hat{D}(G, \mathcal{P}) &= \sum_{j \in J(C)} (\mathcal{P}(x_j) - G(x_j)) N_j q_j (\mathcal{P}(x_j) - G(x_j)) \nonumber \\ &= \langle G - \mathcal{P}, G - \mathcal{P} \rangle = \hat{R}(\mathcal{P}, k, P) \tag{4-4-5} \end{align} $$

$$ = R(\mathcal{P}, k, P) \text{ comme on vient de le démontrer.} $$
Estimateur de la variance $\sigma^2$

Avec les hypothèses du modèle $\Omega$ :

$$ \operatorname{Var} e_{tj} = \frac{y_t}{q_j} \sigma^2. $$

Avec celles du modèle $\tilde{\Omega}$ :

$$ \operatorname{Var} \hat{e}_j = \frac{1}{N_j q_j} \sigma^2 $$

$\sigma^2$ est théoriquement le même dans les deux cas.

$R(\mathcal{P}, k, P)$, qui exprime la distance du polynôme $\mathcal{P}$ à la classe $C_k$ dans le premier cas et à la suite $G$ dans le second, est un estimateur biaisé de $\sigma^2$.

On obtient un estimateur sans biais (c'est-à-dire tel que son espérance mathématique soit $\sigma^2$) en divisant $R(\mathcal{P}, i, P)$ par le nombre de degrés de liberté.

1) Sous $\Omega$ :

$$ R(\mathcal{P}, k, P) = \sum_{y_t \in C_k} \frac{1}{y_t} \sum_{j \in J(y_t)} q_j \left[ y_t(x_j) - \mathcal{P}(x_j) \right]^2 $$

Le nombre de degrés de liberté est la dimension du sous-ESPACE ORTHOGONAL À CELUI ENGENDRÉ PAR $\{ \sum_{i=0}^{d} c_i P_i(x_j) \}$ dans l'espace engendré par les $y_t(x_j)$.

Dimension de l'espace entier : $\sum_{y_t \in C_k} J(y_t) = N$ ($J(y_t) = \text{card } J(y_t)$).

Dimension de l'espace engendré par $\{ \sum_{i=0}^{d} c_i P_i(x_j) \}$, $j \in J(C)$, $d=1$.

d'où le nombre de degrés de liberté de $R(\mathcal{P}, k, P)$ sous $\Omega$ :

$$ N - (d+1) $$

(4-4-6)

2) Sous $\tilde{\Omega}$ :

$$ R(\mathcal{P}, k, P) = \sum_{j \in J(C)} N_j q_j (\mathcal{P}(x_j) - G(x_j))^2. $$
L'espace entier est engendré par {$G(x_j)$}$_{J \in J(C)}$ , il est donc de dimension $card(J(C)) = n$ (par hypothèse d'indépendance sur les $\hat{e}_j$)

d'où le nombre de degrés de liberté de $R(\mathcal{P},k,P)$ sous $\hat{\Omega}$ :

$$n - (d+1) \qquad (4-4-7)$$

On a : $n < N$.

Les estimateurs sans biais de $\hat{\sigma}^2$, notés $s^2$ sous $\hat{\Omega}$ et $\hat{s}^2$ sous $\hat{\Omega}$ valent donc :

$$s^2 = \frac{R(\mathcal{P},k,P)}{N - (d+1)} ; \quad \hat{s}^2 = \frac{R(\hat{\mathcal{P}},k,P)}{N - (d+1)} \qquad (4-4-8)$$

On en conclut que le modèle $\hat{\Omega}$ n'est équivalent à $\hat{\Omega}$ que pour l'estimation des $c_i$.

### 4.4.2. - Tests sous l'hypothèse de normalité

Soit le modèle $\hat{\Omega}$ (4-4-2) et supposons de plus que

$$e_{tj} : N(0, \frac{y_t}{q_j} \sigma^2)$$

$$e_{tj} = \left(\frac{q_j}{J_y}\right)^{1/2} \cdot e_{tj} : N(0,\sigma^2) \text{ et sont indépendantes} \Rightarrow$$

$$\sum_{t,j} (e_t)^2 = R(\hat{\mathcal{P}},i,P) : \sigma^2 \chi^2_{N-(d+1)} \qquad (4-4-9)$$

#### Test de validité du modèle $\hat{\Omega}$

Ce test permet de savoir si l'approximation d'une classe de mesures par un polynôme de degré d est réaliste; il donne une idée de la "ressemblance" des éléments de la classe avec le polynôme de régression.

Intervalle de confiance pour chaque $y_t(x_j)$ :

$$y_t(x_j) - \mathcal{P}(x_j) \text{ suit la loi } N(0, \frac{y_t}{q_j} \sigma^2).$$
$$ \text{Posons } \gamma_{tj} = \frac{y_t}{q_j} \sigma^2 . $$

Soit $\Gamma_{tj}$ un estimateur sans biais de $\gamma_{tj}$.

$$ \Gamma_{tj} = \frac{y_t}{q_j} \cdot \frac{R(\hat{\mathcal{P}}_{k,P})}{N - (d+1)} \quad 4-4-10 $$

Alors $\frac{(N-d-1)\Gamma_{tj}}{\gamma_{tj}}$ suit la loi du $\chi^2_{N-d-1}$ et est indépendant de:

$$ \frac{y_t(x_j) - \mathcal{P}(x_j)}{(\gamma_{tj})^{1/2}} $$ qui suit la loi normale réduite $N(0,1)$ d'où le quotient :

$$ \frac{y_t(x_j) - \mathcal{P}(x_j)}{(\Gamma_{tj})^{1/2}} : \frac{N(0,1)}{(\chi^2_{N-d-1} / (N-d-1))^{1/2}} = T_{N-d-1} \quad 4-4-11 $$

$$ \text{Remarque : } \sum_{\substack{y_t \in \mathcal{C} \\ j \in \mathcal{J}(y_t)}} \frac{(y_t(x_j) - \mathcal{P}(x_j))^2}{\Gamma_{tj}} = \frac{R(\hat{\mathcal{P}}_{k,P})}{R(\hat{\mathcal{P}}_{k,P})} . $$

4-4-11 implique que $\mathcal{P}(x_j)$ doit se trouver dans l'intervalle de confiance suivant, avec la probabilité $1-\alpha$ :

$$ P \left[ \frac{|y_t(x_j) - \mathcal{P}(x_j)|}{((N-d-1)\Gamma_{tj})^{1/2}} < t_{\alpha/2, (N-d-1)} \right] = 1-\alpha \quad 4-4-12 $$

ou :

$$ P \left[ y_t(x_j) - ((N-d-1)\Gamma_{tj})^{1/2} \quad \begin{array}{l} t_{\alpha/2, (N-d-1)} < \mathcal{P}(x_j) \\ < y_t(x_j) + ((N-d-1)\Gamma_{ij})^{1/2} \quad t_{\alpha/2, (N-d-1)} \end{array} \right] = 1-\alpha $$

L'intervalle est donc de longueur

$$ 2 \left\{ (N-d-1) \Gamma_{tj} \right\}^{1/2} t_{\alpha/2, (N-d-1)}. $$

### Comparaison multiple par la méthode de Bonferroni
(cf. Miller [11])

Le modèle $\Omega$ est jugé valable – et la formule suivante donne une idée de sa précision – lorsque :
$$
P \left\{ \bigcup_{\substack{j \in J(y) \\ y \in C}} \left[ \left| y_t(x_j) - \mathcal{F}(x_j) \right| < \sqrt{\Gamma_{tj}} \cdot t_{\alpha/2}, (N-d-1) \right] \right\} = (1-\alpha)^N > 1 - Na
$$

Interprétation statistique du critère W

$$
W(v) = \sum_{i=1}^{k} R(\hat{\mathcal{F}}, i, p) \quad (\text{cf. §})
$$

De plus (formule 4-4-9)

$$
\frac{1}{\sigma_i^2} R(\hat{\mathcal{F}}, i, p) \text{ suit la loi } \chi_{N_i - (d+1)}^2
$$

$$
où N_i = \sum_{y \in C_i} J_y.
$$

$$
W(v) : \sum_{i=1}^{k} \sigma_i^2 \chi_{N_i - (d+1)}^2
$$

(4-4-13)

4.5. - Lien entre l'aspect statistique et l'aspect algébrique

Le problème peut donc se résoudre en ne faisant intervenir que deux éléments arbitraires.

1) Le choix d'une métrique sur $\mathbb{R}^m$

2) La construction d'une base (dépendant de la métrique choisie et du support) qui n'est qu'une étape du calcul de régression et n'influe pas sur la solution.

Si l'on choisit M = I, alors on utilise le critère classique des
moindres carrés.

Cela revient, si on veut rattacher cet aspect à l'aspect statistique
à faire l'hypothèse var (e<sup>J</sup>) = σ<sup>2</sup> I.

Si l'on adopte "l'idée de la régression pondérée" (cf. bibliographie
[5] chap. 8) on peut munir ℝ<sup>m</sup> d'une métrique différente de I et minimiser :

$$
\sum_{j \in J_k} M(j) \left( G^j - \sum_{\alpha=0}^{n} \Lambda^{j}_{\alpha} C_{\alpha} \right)^2
$$

C'est ce que nous avons fait en munissant ℝm de la métrique V⁻¹
(correspondant à l'hypothèse statistique var eJ = σ² V) Vj·V⁻¹(j) peut s'inter-
préter, comme l'inverse de var G^j, ce qui est choisi couramment comme poids associé à l'observation j.

Que l'on aborde le problème sous l'aspect statistique ou sous l'aspect algébrique on est amené à introduire la même hypothèse arbitraire : hypothèse sur la variance dans le 1er cas, choix d'une métrique dans le second.

On a posé :

$$ \forall j \in J_k \quad M(j) = v^{-1}(j) = \mu_j v_{j,k} $$

$$ \text{var } G^j = \sigma^2 / \mu_j v_{j,k} $$

en tenant compte des définitions des "distances" D et R et de la définition du centre de gravité sous la contrainte des données manquantes.

La pondération $\mu_j$ de l'ensemble des paramètres permet de modifier la métrique en respectant ces contraintes (le soin de choisir $\mu_j$ est laissé à l'utilisateur).

# 5. - LA BASE D'INTERPOLATION

## 5.1. - Choix de la base

Cette base peut être constituée des monômes $l, x, x^2, \dots, x^n$ de polynômes de degré inférieur ou égal à n, ..., etc, ...

Notre choix motivé par le désir de simplifier les calculs notamment l'inversion d'une matrice, sera une base de polynômes orthogonaux.

Remarquons que ce problème d'inversion se pose pour chaque classe, que la métrique est associée à la classe $(v_k^{-1})$ et que par conséquent il faudra construire une base pour chaque classe ce qui est relativement peu gênant dans la mesure où le calcul des estimateurs $\hat{c}_k$ se fait en même temps.

A la base correspond une matrice $\Lambda_k$ de taille $|J_k| \times (n+1)$ des valeurs sur le support $J_k$ des différents polynômes $\Lambda_o, \dots, \Lambda_n$ de cette base.

$$ \text{Soit } v_k^{-1} \text{ la matrice } |J_k| \times |J_k| \text{ induisant une métrique sur } \mathbb{R}^{|J_k|} $$
$$ v_k^{-1} = \begin{bmatrix} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & v_{j,k} \mu_j & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix} $$

$$ \text{l'estimateur } \hat{c}_k = (\Lambda'_k v_k^{-1} \Lambda_k)^{-1} \Lambda_k $$

(Dans la suite nous allègeons les notations en supprimant l'indice k, mais il reste entendu que nous raisonnons pour une classe fixée).

La matrice $S = (\Lambda' V^{-1} \Lambda)$ est une matrice $(n+1) \times (n+1)$ de terme général $S_{\alpha\beta} = \sum_{j \in J_k} \Lambda_{\alpha}^j \Lambda_{\beta}^j v_{j,k} \mu_j = \sum_{j \in J} \Lambda_{\alpha}^j \Lambda_{\beta}^j v_{j,k} \mu_j$ (car $v_j \notin J_k$ $v_{j,k} = 0$) que l'on note encore :

$S_{\alpha\beta} = < \Lambda_{\alpha}, \Lambda_{\beta} > V^{-1}$ ou encore $< \Lambda_{\alpha}, \Lambda_{\beta} >$ on voit que si les polynômes

$\{\Lambda_{\alpha}, \alpha = 0, \dots, n\}$ sont orthogonaux sur le support $J$ pour la métrique associée à $V^{-1}$, la matrice $S$ sera diagonale.

De plus comme $cov(\hat{c}) = E[(c-\hat{c})(c-\hat{c})']$, si $S$ est diagonale on a $\forall_{\alpha}$ var $\hat{c}_{\alpha} = \sigma^2/<\Lambda_{\alpha}, \Lambda_{\alpha}>$ et $\forall_{\alpha}, \forall_{\beta} \neq \alpha$ cov $(\hat{c}_{\alpha}, \hat{c}_{\beta}) = 0$ i.e. les coefficients sont non corrélés.

### 5.2. - Construction de la base de polynômes orthogonaux sur le support $J_k$ (cf. [3])

**Rappel** : On définit sur l'ensemble des polynômes considérés comme éléments de $\mathbb{R}^J$ un produit scalaire sur $\mathbb{R}^J$

$< P^J, Q^J > = \sum_{j \in J} p^j Q^j v_{j,k} \mu_j$ associé à la métrique induite par la classe $k.*$

<footnote>* Comme $v_j \notin J_k$, $I_{j,k} = \emptyset$ et $v_{j,k} = 0$</footnote>

<footnote>$< P^J, Q^J > = \sum_{j \in J_k} p^j Q^j v_{j,k} \mu_j$</footnote>
On encore si P, connu par ses valeurs sur le support $J_k$ est consi-
déré comme un élément de $\mathbb{R}^{|J_k|}$ ainsi que Q.

On note :

$$
\langle P, Q \rangle = P' V^{-1} Q .
$$

On cherche une base $\Lambda_0, \Lambda_1, \dots, \Lambda_n$ de l'espace vectoriel, de
dimension $n+1$, de l'ensemble des polynômes de degré inférieur ou égal à $n$,
telle que :

$\forall o < \alpha < n \Lambda_{\alpha}$ est un polynôme de degré $\alpha$

$$
\text{et } \forall \alpha, \beta \ 0 \le \alpha < \beta < n \Rightarrow \langle \Lambda_{\alpha}, \Lambda_{\beta} \rangle = 0
$$

$\Lambda_o$ est évidemment une constante que l'on posera égale à 1.

On peut démontrer la proposition suivante :

(cf. Forsythel [12] ; Peck [13])

Proposition 3 :

1) $\forall \alpha \ge 0$ le polynôme $\Lambda_{\alpha+1}$, qui s'exprime en fonction des polynômes $\Lambda_\alpha$ et $\Lambda_{\alpha-1}$ par la relation

(R) $\Lambda_{\alpha+1}(x) = (x-a_{\alpha}) \Lambda_{\alpha}(x) - b_{\alpha-1} \Lambda_{\alpha-1}(x)$ est orthogonal aux précédents à condition que $a_{\alpha}$ et $b_{\alpha-1}$ aient pour valeurs respectives :

$$
a_{\alpha} = \frac{\langle \Lambda_{\alpha}, \Lambda_{\alpha} \rangle}{\langle \Lambda_{\alpha}, \Lambda_{\alpha} \rangle}; \quad b_{\alpha-1} = \frac{\langle \Lambda_{\alpha}, \Lambda_{\alpha} \rangle}{\langle \Lambda_{\alpha-1}, \Lambda_{\alpha-1} \rangle}
$$

2) $\{\Lambda_o, \{\Lambda_{\alpha+1}, \alpha=0,\dots,n-1\}\}$ est une base de l'espace vectoriel de l'ensemble des polynômes de degré inférieur ou égal à n.

3) Cette base étant choisie, l'estimateur $\hat{c}$ des coefficients du polynôme de régression s'exprime de manière simple par :

$$
\hat{c}_{\alpha} = \frac{\langle \Lambda_{\alpha}, G \rangle}{\langle \Lambda_{\alpha}, \Lambda_{\alpha} \rangle} , \forall \alpha
$$
Le polynôme de degré n minimisant la variance sera donc

$$
\hat{L}_{n+1} = \sum_{\alpha=0}^{n} \hat{C}_{\alpha} \Lambda_{\alpha} = \sum_{\alpha=0}^{n} \frac{\langle \Lambda_{\alpha}, G \rangle}{\langle \Lambda_{\alpha}, \Lambda_{\alpha} \rangle} \Lambda_{\alpha}
$$

Remarque : Si l'on désire un polynôme de degré supérieur $n+1$, on fera la même construction jusqu'à $\Lambda_{n+1}$ et on a :

$$
\hat{L}_{n+1} = \hat{L}_n + \frac{\langle \Lambda_{n+1}, G \rangle}{\langle \Lambda_{n+1}, \Lambda_{n+1} \rangle} \Lambda_{n+1} \quad (5-1)
$$

L'égalité ci-dessus est utilisée lorsque les tests sur le degré du polynôme montrent que le degré atteint est insuffisant. Elle permet d'élever le degré sans refaire tous les calculs (cf. 6.2).

6. - QUALITE DE LA REGRESSION. TEST D'AJUSTEMENT DU DEGRE DU POLYNOME

6.1. - Le problème

Le problème que nous avons résolu ici est celui de trouver le
polynôme (C₀, ..., Cₙ) tel que ||G - ΛC||² V⁻¹ soit minimum.

A degré n fixé, ce problème admet une solution Ĉ unique. Cependant il est intéressant de se demander quels sont les coefficients Cᵢ qui ont le "meilleur pouvoir explicatif", quels sont ceux dont le "pouvoir explicatif" est négligeable (en un sens qui reste à préciser).

En d'autres termes on peut se demander si l'estimation de G par :

$$
(1) \quad \hat{G} = \sum_{\alpha=0}^{n} \Lambda_{\alpha} \hat{C}_{\alpha} \text{ notée } G_n
$$

est "significativement" meilleure que son estimation par :

$$
(2) \quad \hat{G} = \sum_{\alpha \in A} \Lambda_{\alpha} \hat{C}_{\alpha} \quad \text{avec } A \subset \{0, 1, \dots, n\}
$$

ou par

$$
(3) \quad \hat{G} = \sum_{\alpha=0}^{r} \Lambda_{\alpha} \hat{C}_{\alpha} \quad \text{avec } r < n \text{ notée } \hat{G}_{r}
$$
Notons que les coefficients $C_\alpha$ des trois estimations précédentes sont identiques. Celà est dû à la propriété d'indépendance statistique des coefficients signalée en 5.2.

Pour juger de la qualité d'une estimation, on peut penser à utiliser un indice normé tel $R = \frac{\|G - \hat{G}\|}{\|G\|}$

En raison de la propriété rappelée ci-dessus on a :

$$ \forall r, \forall n, r < n \Rightarrow \hat{G}_n = \hat{G}_r + \sum_{\alpha=r+1}^{n} \Lambda_{\alpha} \hat{C}_{\alpha} $$

et par suite :

$$ (R_n)^2 = \frac{\|G - \hat{G}_n\|^2}{\|G\|^2} < (R_r)^2 = \frac{\|G - \hat{G}_r\|^2}{\|G\|^2} $$

La meilleure estimation au sens de R est donc toujours celle de degré supérieur (1).

On voit que ce critère n'est pas très satisfaisant. Il tient compte de la précision de l'estimation, mais la plus ou moins grande complexité du modèle n'est pas considérée; or on préfère un modèle plus simple (i.e. ici une estimation par un polynôme de plus bas degré) à condition de ne pas modifier "significativement" le "pouvoir explicatif" de l'estimation.

Si l'on veut savoir s'il y a intérêt à inclure les n-r polynômes de degré supérieur à r et inférieur ou égal à n dans la régression le calcul de

$$ F = \frac{\| \hat{G}_n - \hat{G}_r \|^2}{\| G - \hat{G}_n \|^2} \times \frac{m - n - 1}{n - r} $$

nous permet d'apprécier le "pouvoir explicatif" de ces n-r polynômes lorsque les r+1 premiers sont utilisés: sa comparaison à la valeur correspondante d'une table du F de Fisher-Snedecor, pour un seuil de confiance fixé, nous permet de juger si l'estimation de G par un polynôme de degré n est significativement meilleure que son estimation par un polynôme de degré r.

<footnote>(1) Le contraire, c'est-à-dire une estimation par r polynômes meilleure qu'une estimation par n, ne pourrait se produire que si les r polynômes n'étaient pas tous contenus dans les n. Ce n'est pas le cas ici.</footnote>
**Remarque** : Le cas le plus général (celui où $\gamma_1$ serait un ensemble quelconque de coefficients $C_i$ et non les r+1 premiers) nous amènerait à considérer les $(2^{n+1} - 1)$ estimations possibles par des polynômes de degré inférieur ou égal à n. Nous nous sommes bornés à tester de manière systématique l'estimation par le meilleur polynôme de degré n contre l'estimation par le meilleur polynôme de degré r ($r < n$).

## 6.2. - Procédure utilisée

On a choisi une technique d'introduction "ascendante" des polynômes qui peut s'effectuer au fur et à mesure de la recherche des coefficients.

On pose $r = n - 1$, un seul polynôme de la base, celui de degré $n$ n'est pas utilisé dans l'estimation de $G$.

$$F = (m - n - 1) \frac{\left| \sum_n \hat{c}_n \right|^2}{\left| \left| G - \hat{G}_n \right| \right|^2}$$

mesure le "pouvoir explicatif" du n<sup>ième</sup> degré conditionnellement lorsque $G$ est estimé par un polynôme de degré $n-1$. On arrête le processus lorsque la dernière variable introduite n'est pas "significative" au sens de $F$.

Cependant, lorsque l'algorithme est modifié pour ajuster le degré du polynôme de cette manière, on risque de se trouver devant une situation de non-décroissance de la suite $u_n$. En effet, si lors du passage d'une itération à la suivante le degré diminue pour certaines classes, les "distances" classes-polynômes correspondantes peuvent augmenter même si l'on est en droit de juger qu'en un certain sens on a une meilleure estimation.

Précisons :

Soient $(i)$ et $(i+1)$ deux itérations successives, $G^{(i+1)}$ le centre de gravité de la classe $P^{(i+1)}$

$$\epsilon(i) = \left| G^{(i+1)} - \hat{G}_n^{(i)} \right|^2 \quad \text{et} \quad \epsilon(i+1) = \left| G^{(i+1)} - \hat{G}_n^{(i+1)} \right|^2$$

Les erreurs correspondant aux estimations par le degré $n$.

En raison des propriétés de l'algorithme, il est certain que :

$$\epsilon(i) \ge \epsilon(i+1) \quad (1)$$
Mais si à l'itération (i+1) le processus d'itération progressive
s'est arrêté au degré r < n i.e. si $\epsilon(i+1)$

$$ = ||G^{(i+1)} - \hat{G}_r(i+1)||^2 \quad \text{comme on a toujours} $$

$$ ||G - \hat{G}_r||^2 > ||G - \hat{G}_r||^2 \quad \text{si } r < n \text{ l'inégalité (1) n'est plus certaine, mais} $$

elle est probable, il est intuitivement évident que les polynômes éliminés
contribuent relativement peu (au sens de R) à l'explication de G. Le test a
été programmé de cette manière (notamment en raison du moindre coût).

Très peu de cas de non-décroissance ont été observés généralement
lors des premières itérations seulement. L'algorithme a toujours convergé. En
effet, même si une augmentation du critère est observée au moment où le degré
diminue, quand les classes et leur degré sont stabilisés (i.e. dans les itéra-
tions proches de la convergence) le critère recommence à décroître et l'algo-
rithme converge.

Il est néanmoins possible de procéder différemment et de poursuivre
le processus d'estimation jusqu'au degré maximum fixé.

## 7. LE PROGRAMME

Le programme LISTYP mis au point s'inspire du programme ANATYP (cf. [10] [11])
en ce qui concerne l'organigramme général, la construction de la partition aléa-
toire de départ, l'affection des individus aux classes, le calcul du critère,
le test de convergence.

Les deux programmes diffèrent essentiellement dans la recherche des noyaux.

Les sous-programmes MATDIA et DIAG effectuant les analyses factorielles
dans ANATYP sont ici remplacés par LPORT, VPORT (cf. [3]) qui construisent une
base de polynômes orthogonaux, calculent les coefficients des polynômes de
meilleure approximation au sens des moindres carrés et estiment les données man-
quantes.

Le sous-programme TEST a pour objet d'ajuster le degré polynôme, son utili-
sation est facultative.
## 8. - APPLICATIONS

### 8.1. - Evolution du nombre de cochenilles du pin sylvestre au cours d'une période de cinquante semaines dans 10 stations agronomiques **

#### Données *

<table><tbody><tr><td rowspan="2">1</td><td>30</td><td>24</td><td>24</td><td>19</td><td>34</td><td>16</td><td>13</td><td>9</td><td>5</td></tr><tr><td>32</td><td>32</td><td>41</td><td>37</td><td>41</td><td>40</td><td></td><td>22</td><td>36</td></tr><tr><td rowspan="2">2</td><td>32</td><td>11</td><td></td><td>22</td><td></td><td>18</td><td>20</td><td>18</td><td>21</td><td>20</td><td>8</td></tr><tr><td>23</td><td>37</td><td>83</td><td>84</td><td>46</td><td></td><td>30</td><td></td><td>30</td><td></td><td>16</td></tr><tr><td rowspan="2">3</td><td>7</td><td>14</td><td></td><td>9</td><td></td><td>13</td><td>12</td><td>13</td><td>8</td><td>10</td><td>13</td></tr><tr><td>12</td><td>20</td><td>19</td><td>15</td><td>18</td><td>22</td><td></td><td></td><td>20</td><td></td><td>12</td></tr><tr><td rowspan="2">4</td><td>19</td><td>12</td><td></td><td>17</td><td></td><td>12</td><td>13</td><td>15</td><td>10</td><td>14</td><td>6</td></tr><tr><td>15</td><td>12</td><td>20</td><td>13</td><td>16</td><td></td><td>12</td><td></td><td>10</td><td></td><td>15</td></tr><tr><td rowspan="2">5</td><td>34</td><td>29</td><td></td><td>35</td><td></td><td>20</td><td>26</td><td>13</td><td>13</td><td>16</td><td>38</td><td>49</td></tr><tr><td>35</td><td>36</td><td>38</td><td>24</td><td>20</td><td>28</td><td></td><td>46</td><td>31</td><td>31</td><td>21</td><td>23</td><td>23</td></tr><tr><td rowspan="2">6</td><td>13</td><td>26</td><td>18</td><td>14</td><td>17</td><td>15</td><td>19</td><td>17</td><td>11</td><td>14</td><td>21</td><td>16</td><td>9</td></tr><tr><td>30</td><td>9</td><td>18</td><td>18</td><td>21</td><td>46</td><td>25</td><td>26</td><td>21</td><td>20</td><td>21</td><td>20</td><td>18</td></tr><tr><td rowspan="2">7</td><td></td><td></td><td></td><td></td><td>9</td><td>4</td><td>6</td><td>3</td><td>3</td><td>5</td><td>17</td><td>18</td><td>7</td></tr><tr><td>7</td><td>10</td><td></td><td></td><td>8</td><td>9</td><td></td><td>35</td><td>28</td><td>21</td><td>82</td><td></td><td></td></tr><tr><td rowspan="2">8</td><td></td><td></td><td></td><td></td><td></td><td></td><td>18</td><td>24</td><td></td><td></td><td>9</td><td>11</td></tr><tr><td>7</td><td>23</td><td>26</td><td>16</td><td>38</td><td>38</td><td>19</td><td>23</td><td></td><td>36</td><td></td><td></td></tr><tr><td rowspan="2">9</td><td></td><td></td><td></td><td>56</td><td>34</td><td>77</td><td>68</td><td>66</td><td>39</td><td>36</td><td>26</td><td>19</td><td>40</td></tr><tr><td>34</td><td>37</td><td>69</td><td>185</td><td>153</td><td>116</td><td>158</td><td>169</td><td>117</td><td>87</td><td>106</td><td></td><td>9</td></tr><tr><td rowspan="2">10</td><td>30</td><td>20</td><td>16</td><td>14</td><td>23</td><td>21</td><td>13</td><td>13</td><td>11</td><td>22</td><td>13</td><td></td><td></td></tr><tr><td></td><td></td><td>41</td><td>34</td><td>81</td><td>61</td><td>31</td><td>37</td><td>18</td><td>10</td><td>37</td><td>22</td><td></td></tr></tbody></table>

Les dix stations correspondent à des conditions agronomiques différentes ; les conditions expérimentales n'ont pas permis d'avoir des séries complètes pour les 50 semaines.

<footnote>* Données fournies par J.C. ONILLON, INRA, station de zoologie et lutte biologique, 06602 - ANTIBES.</footnote>

<footnote>* Ce travail a été réalisé avec Melles Cloutier et Lestrade [10].</footnote>
a) le nombre maximum de classes (NKMAX) est fixé à 3 et le degré maximum des polynômes à 2.

Pour les 5 tirages, la convergence a été atteinte en 4 itérations maximum et on a obtenu 5 fois la même partition. Les valeurs des polynômes qui ont été calculées pour les différentes classes nous donnent la représentation

b) Le nombre maximum de classes est NKMAX = 3 et le degré maximum des polynômes est 5.

On obtient pour les 5 tirages, convergeant en six itérations maximum les mêmes classes que précédemment.

Les polynômes associés aux différentes classes sont représentés sur la figure II.5. Afin de comparer les ajustements obtenus pour les degrés 2 et 5 nous avons représenté les valeurs réelles de tous les individus appartenant à une même classe et ceci pour chaque classe (voir. fig. II.1, II.2, II.3). L'ap-proximation par les polynômes de degré 5 semble plus raisonnable que celle par les polynômes de degré 2 car les courbes (II.5) ajustent mieux que les courbes (II.4) le profil des individus des différentes classes. Le critère correspon-dant à n = 5 est inférieur au critère obtenu pour n = 2, ce qui est normal (cf. propriété signalée au § 6); cependant, comme on l'a remarqué en 6.1, le critère ne permet pas de comparer des estimations par des polynômes de degré différent. De plus, il peut être illusoire de vouloir fixer un degré trop élevé si les mesures n'ont pas une précision suffisante. L'utilisation du test décrit en 6. devrait aider au choix du degré.

c) L'individu 9 formant toujours une classe à part, on l'a éliminé et refait une analyse avec KMAX = 2; d = 5. On a obtenu pour 4 tirages sur 5 le résultat des analyses précédentes qui, de plus correspond au meilleur critère (voir II.6).

Nous voyons que dans les données étudiées ici, les individus 1,2,5,10 d'une part et les individus 3,6,4,7,8 d'autre part, se trouvent toujours dans les mêmes classes. Le choix, parmi les optima locaux, n'a donc posé aucun problème.

**Remarque** : Influence du support J

Ici J = {1,2,3,...,49,50} est un codage discret du temps,
(numérotation de cinquante semaines consécutives) dans la deuxième application
c'est un codage discret de la profondeur, (numérotation de vingt-cinq niveaux
distants de un mètre). La méthode a été conçue à l'origine pour de telles don-
nées i.e. à support ordonné et régulièrement espacé. Elle est sans doute sen-
sible à l'influence du support, notamment à l'espacement des valeurs de la
variable; cette influence reste à étudier.

8.2. - *Données géologiques : Teneur en minerai (Fer, nickel, cobalt) de 120 parcelles à 25 niveaux de profondeur.*

Le support J = (1,2, ..., 25) représente vingt cinq niveaux équi-
distants les uns des autres d'un mètre. Les données étant complètes, nous avons
voulu essayer de déterminer un "seuil de validité"de la méthode i.e. et de
trouver le pourcentage de données manquantes pour lequel les résultats sont très
différents de ceux obtenus pour les données complètes. Pour cela après traite-
ment des données complètes, on a supprimé aléatoirement 10%, puis 20% de mesures
et appliqué la méthode de lissage typologique.

On a effectué 5 tirages, demandé 5 classes et fixé le degré des
polynômes à 5*. On donne les résultats concernant le nickel dans les tableaux
1 à 3 et les graphiques correspondants.

On remarque, que les individus, que nous avions à classer, se re-
groupent très souvent dans les mêmes classes, ceci étant vérifié, nous l'avons
vu précédemment, pour les cochenilles mais aussi pour les données géologiques.
Si l'on compare les résultats obtenus pour ces données complètes, avec 10% puis
avec 20% de trous, on constate que les classes varient assez peu d'une analyse
à l'autre, ce qui semble prouver une certaine efficacité de la méthode.

<footnote>* Il nous a semblé, a priori, souhaitable de choisir un degré assez élevé (cf. 6.1). En effet, les estimateurs des coefficients des monômes de plus fort degré sont de plus en plus petits ce qui modifie peu les valeurs manquantes estimées. A l'opposé, si l'on demande un degré trop faible, il peut arriver que soit négligée une certaine quantité importante pour la détermination des polynômes.</footnote>
## 9. - CONCLUSION

Le travail présenté ici ne constitue qu'une première approche du traitement des données manquantes au moyen de l'algorithme des "Nuées Dynamiques".

On a considéré que les lignes du tableau des données étaient les valeurs connues sur un certain support de fonctions polynômiales d'une variable.

On s'est limité au cas où le support J, de définition expérimentale de ces fonctions, est discret, naturellement ordonné et à valeurs régulièrement espacées.

On a choisi le procédé de "meilleure approximation au sens des moindres carrés", parce qu'il est : "... peu sensible au choix du support dès que le nombre de points est suffisant et qu'ils sont raisonnablement répartis..." et parce que c'est "... une technique ... préférée chaque fois qu'il faudra approcher une fonction définie par un procédé expérimental". (cf. LEGRAS [7]).

Dans le cadre de cette approche beaucoup de problèmes intéressants restent à étudier : influence du pourcentage de données manquantes et détermination d'un "seuil" de validité de la méthode : influence d'un support "irrégulier"; influence de la pondération des paramètres (i.e. des métriques choisies (cf. 4.4.); utilisation d'autres bases de polynômes (Legendre etc...) ou même d'autres fonctions (gaussiennes, sinusoïdales, etc...); utilisation d'autres procédés d'approximation (au sens de Tchebitchev, en moyenne quadratique, par sommes partielles de séries entières, etc...).

Enfin, le problème du choix des paramètres d'entrée de la méthode, c'est-à-dire ici essentiellement le choix du nombre de classes, puisque le degré peut-être ajusté, reste entier.

Ce choix, arbitraire en l'absence d'indications sur les données, influe sans doute sur les résultats.
<figure><img src="image_6.png" /></figure>

Classe 4
Classe 2
Classe 3

Polyèmeres de dégâts 2
{ 3, 4, 6, 7, 8
[ ]

<figure><img src="image_7.png" /></figure>

Polyèmeres de dégâts 5 a
{ 4, 2, 5, 40
[9]
{ 3, 4, 6, 7, 8

Hormon de croissance

<figure><img src="image_8.png" /></figure>

Classe 1 : Staphens { 1, 2, 5, 10 }

Classe 2 : Staphens { 1, 2, 5, 10 }

Temps

(1) D. 6 en semaines
Teneur en Nickel N

NICKEL
(données complètes)
Polynômes de après 5
(sans test d'ajustement du
Test d'ajustement du
défini)

Teneur en Nickel

Teneur en Nickel

<figure><img src="image_9.png" /></figure>

Sondages géologiques : teneur en Nickel

0 % de données manquantes

Nombre maximum de classes demandées : 5

Dégré maximum des polynômes : 5 (sans test d'ajustement)

Partition correspondant au "meilleur" tirage sur 5

effectués, (valeur du critère : W = 5.8885)

Niveau

15

Fig. II-8

Fig. II-14

25

<figure><img src="image_10.png" /></figure>

Sondages géologiques : teneur en Nickel
10 % de données manquantes

Nombre maximum de classes demandées : 5

Dégré maximum des polynômes : 5 (sans test d'ajustement)

Partition correspondant au "meilleur" tirage

effectués (valeur du critère : W = 5.5897

<figure><img src="image_11.png" /></figure>

Sondages géologiques : teneur en Nickel

20 % de données manquantes

Nombre maximum de classes demandées : 5

Dégré maximum des polynômes : 5 (sans test d'ajustement)

Partition correspondant au "meilleur" tirage sur 5

effectués (valeur du critère : W = 5.56287)

Niveaux de
Profondeur

Niveaux de
Profondeur
BIBLIOGRAPHIE

[1] GRAYBILL - An introduction to linear statistical models
Mc GRAW HILL 1961

[2] SCHEFFE - The Analysis of variance
WILEY 1959

[3] Bibliothèque mathématique de programmes FORTRAN CII

[4] Y.V. LINNICK - Méthode des moindres carrés
DUNOD 1963

[5] CAILLIEZ, MAILLE, NAKACHE, PAGES - Analyse des données multidimensionnelles -
C.3E 1971

[6] DIDAY E. - Thèse de Doctorat d'état - 1972

[7] DIDAY E. - Introduction à l'analyse factorielle typologique -
Rapport de Recherche n° 27 - IRIA-Laboria

[8] OK-SAKUN Y. - Thèse de 3ème cycle 1975

[9] LEGRAS - Méthodes et techniques de l'analyse numérique
DUNOD 1971

[10] CLOUTIER et LESTRADE -
Rapport de stage D.E.A. Laboratoire de statistique mathématique

[11] MILLER R.G. Jr - Simultaneous statistical inference
Mc Graw Hill, New-York 1966

12 FORSYTHE G.E. - Generation and use of orthogonal polynomials for data fitting
with a digital computer. 1957
J. SIAM 5, 74-88

13 PECK J.E.L. - Polynomial curvefitting with constraint
Siam Review 4, 2, 135-141 - 1962