Dw-03477

RR257

RO 2961

Côle: 115-R 257

IRIA

laboria

laboratoire de recherche
en informatique
et automatique

RÉGRESSION TYPOLOGIQUE

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

Christian CHARLES

Rapport de Recherche N° 257

Octobre 1977

Institut de Recherche
d'Informatique
et d'Automatique

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

L44p
RÉGRESSION TYPOLOGIQUE

Christian Charles

Résumé :

Dans les années récentes l'attention s'est portée sur les problèmes de discontinuités en régression, discontinuités intervenant en des points inconnus. Cette étude approche ces problèmes par l'uti- lisation de la classification non hiérarchique des observations. Ainsi nous abordons l'hypotèse la plus complexe dans laquelle les régimes de régression qu'il s'agit de découvrir peuvent se rompre et reprendre à des instants quelconques. Les aspects inverses généralisés et ridge régression sont étudiés. Le programme de régression typologique appelé REGTYPO est présenté. On donne enfin une méthode de prévision associée aux modèles linéaires locaux qui devrait être plus précise que les techniques habituellement utilisées pour le modèle linéaire global.

Abstract :

In recent years much attention has been focussed on the problem of discontinuous shifts in
regression regimes at unknown points in the data series. This paper approaches this problem by
using non hierarchical clustering. So we tackle the most complex problem in which the different
regressions to be found may break off and start afresh at any instant. The aspects of generalized
inverse and ridge regression are also studied. We present also the program REGTYPO correspon-
ding to this new method. At last, we give a forecasting approach of the local linear models, which
may be more precise than the technics usually used to resolve a global linear model.
# SOMMAIRE

## INTRODUCTION

1. Rappels théoriques et inverses généralisées

2. L'algorithme

2.1. Les notations employées

2.2. Le critère

2.3. Les fonctions f et g

2.4. Les suites $u_n$ et $v_n$

2.5. La convergence de l'algorithme

3. Cas de variables non orthogonales

3.1. Comparaisons entre inverses généralisées et ridge régression

3.2. Schéma général de l'algorithme

3.3. Le programme informatique REGTYPO

4. Elimination de l'effet de cylindre

5. Choix optimal du nombre de classes

6. Méthode de prévision associée aux modèles linéaires locaux

## CONCLUSION

## BIBLIOGRAPHIE

## ANNEXES
# INTRODUCTION

Lorsqu'on utilise des séries chronologiques, on se pose souvent des problèmes de stabilité temporelle, c'est-à-dire qu'on désire tester à partir des observations si une même loi est bien adaptée pour représenter le même phénomène à des époques différentes. Par exemple : une fonction de consommation macroéconomique est-elle restée la même avant et après la guerre ?, les investissement industriels en France ont-ils obéi à la même loi avant et après la mise en place du marché commun ?, la fonction d'importation d'un pays anciennement colonisé s'est-elle modifiée depuis l'indépendance ?...

Dans tous les cas ci-dessus, le test de stabilité porte sur deux périodes consécutives. Mais il pourrait tout aussi bien porter sur des périodes entrenêlées. Par exemple : la relation de Philips (liant salaires prix, chômage) est-elle la même suivant les diverses phases de la conjoncture : expansion et récession ?

Dans ces exemples, les groupes d'observations sont définis a priori : consommations avant la guerre, consommations après la guerre, investissements avant le marché commun, investissements après le marché commun... La stabilité temporelle est donc mesurée grâce au test d'É Chow (1960) dans le cas de modèle linéaire.

Un problème plus difficile à résoudre est celui qui suppose que la série chronologique a été générée par plusieurs régimes de régression dont on se propose de rechercher les instants de discontinuités. Il s'agit alors de définir les différents groupes d'observations et les régimes de régression associés.

Pour reprendre l'exemple des investissements industriels, il ne s'agit plus de vérifier si le marché commun a modifié le comportement des entreprises, mais plutôt, sachant que ce comportement a été modifié au cours du temps, de détecter le ou les instants séparant les différents régimes (et si on
le désire les rattacher à une modification intervenue dans l'économie du pays) ainsi que ces régimes. Diverses méthodes se proposant de résoudre ce problème ont été proposées. On peut les regrouper en deux approches distinctes :

1) -QUANDT (1958) exprime la fonction de vraisemblance de l'échantillon comme une fonction du temps et la maximise par rapport aux moyennes $\beta_1$ et variances $\sigma_1^2, \sigma_2^2$ des deux régimes ainsi que par rapport au temps.

- FAIRLEY et HINICH (1970) supposent que le point de rupture est équiprobablement réparti sur l'ensemble de l'échantillon et dérive un test de vraisemblance.

- BROWN et DURBIN (1968) proposent de définir des écarts récursifs $u_r ; r = e+1, T$ (où e est le nombre de variables explicatives et T le nombre d'observations) obtenus en utilisant la première rangée de y et x (du modèle $y = x \beta + u$) et une estimation $\beta$ basée sur les (r-1) premières rangées de y et x. Ces écarts sont transformés en variables distribuées suivant une loi N(0, $\sigma^2$ I) et leurs sommes cumulées sont testées par rapport à zéro.

2) Une formulation plus complexe supposant que les régimes puissent se rompre et reprendre à différents instants du temps devient rapidement impraticable puisqu'elle nécessite l'évaluation de $2^T$ fois la fonction de vraisemblance au lieu de T fois.

Deux nouvelles approches ont été explorées :

- FAIR et JAFFEE (1972), GOLDFELD et QUANDT (1972) supposent l'existence d'une variable supplémentaire z qui permet de classer les individus $(y_r, x_r)$ suivant un critère du type :

$(y_r, x_r)$ est généré par le régime 1 si $z_r < z_o$

$(y_r, x_r)$ est généré par le régime 2 si $z_r > z_o$ où $z_o$ est connu.
MAC GEE et CARLETON (1970 )

utilisant un algorithme de classification hiérarchique et se propose de rechercher les $t^{(m)}$, m = 1,M (M connu) tels que dans le cas des moindres carrés ordinaires sur e variables explicatives, le modèle supposé approprié s'écrive :

$$
\left\{
\begin{array}{l}
y_t = \beta_0^{(1)} + \beta_1^{(1)} x_{1t} + \dots + \beta_e^{(1)} x_{et} + \varepsilon_t \\
\qquad \forall t < t^{(1)} \\[1em]
y_t = \beta_0^{(2)} + \beta_1^{(2)} x_{1t} + \dots + \beta_e^{(2)} x_{et} + \varepsilon_t \\
\qquad \forall t^{(1)} \le t < t^{(2)} \\[1em]
y_t = \beta_0^{(M)} + \beta_1^{(M)} x_{1t} + \dots + \beta_e^{(M)} x_{et} + \varepsilon_t \\
\qquad \forall t^{(M-1)} \le t \le t^{(M)}
\end{array}
\right.
\qquad \forall t^{(M)} = T
$$

avec $\epsilon_t \sim N(0, \sigma^2)$

La méthode que nous proposons réalise la synthèse des techniques de
régression et de classification automatique non hiérarchique. La classi-
fication non hiérarchique des observations, effectuée selon la méthode
des nuées dynamiques (E. DIDAY 1972 ' est réalisée simultanément
avec le calcul des hyperplans de régression. Nous obtenons ainsi une
typologie de l'ensemble des individus et un ensemble d'hyperplans de
régression locaux associés à ces types, d'ou le nom de régression
typologique ou régression locale (E. DIDAY 1976 )

La méthode proposée construit à partir d'une partition P<sup>(o)</sup> (estimée
ou tirée au hasard) en K classes (K fixé a priori) des observations une
suite L<sup>(n)</sup> de K samples de régression, et une
suite $P^{(n)}$ de K-partition de façon à minimiser un critère mesurant l'inadéquation des hyperplans aux groupes formés (on montre le lien entre ce critère et le test de Chow) en utilisant de façon alternative une fonction d'identification (constitution de P) et une fonction de représentation (constitution de L). Nous résolvons ainsi le problème indiqué en 2) par un algorithme moins sujet à caution que ceux supposant l'existence d'une variable supplémentaire $Z$ dont il faut fixer le seuil $Z_o$, et moins coûteux que celui de Mc Gee
# 1. Rappels théoriques et inverses généralisées

Dans ce paragraphe nous proposons une méthode permettant de détecter K régimes locaux de régression et les K hyperplans associés de façon à minimiser un critère mesurant l'inadéquation des sous échantillons à leurs représentations linéaires.

Nous rappelons tout d'abord les résultats classiques de la régression et introduisons les inverses généralisées.

Le modèle de régression linéaire multiple s'écrit habituellement :

$$ \forall n = 1, \dots, N $$

$$ y_n = b_1 x_{n1} + b_2 x_{n2} + \dots + b_m x_{nm} + \varepsilon_n $$

$$ \text{avec } E \varepsilon_n = 0 $$

$$ E \varepsilon_n \varepsilon_{n'} = \delta_{nn'} \sigma^2 \qquad 0 < \sigma^2 < \infty $$

ou matriciellement :

$$ Y = Xb + \varepsilon \qquad Y \text{ vecteur } (N,1) $$

$$ E \varepsilon = 0 \qquad \text{avec } X \text{ matrice } (N,m) $$

$$ V \varepsilon = \sigma^2 I_N \qquad \begin{array}{l} b \text{ vecteur } (m,1) \\ \varepsilon \text{ vecteur } (N,1) \end{array} $$

La méthode des moindres carrés ordinaires consiste à choisir pour estimateur $\hat{b}$ de $b$ le vecteur qui rend minimum dans l'échantillon la somme des carrés des écarts entre les valeurs observées $y_n$ et les valeurs ajustées $\hat{y}_n$ de la variable dépendante, ce qui s'écrit :

$$ \underset{b}{\operatorname{Min}} \quad \underset{n=1}{\sum} \quad \underset{e=1}{\sum} \quad \underset{b_e}{\underset{ne}{\sum}} \quad \begin{array}{l} N \\ m \\ 2 \end{array} \quad (y_n - \sum_{e=1}^{m} b_e x_{ne}) $$
ou matriciellement :

$$
\begin{array}{l}
\underset{b}{\operatorname{Min}} \quad (Y-Xb)'(Y-Xb) = \underset{b}{\operatorname{Min}} S(b) \\
b \qquad b
\end{array}
$$

Les conditions du premier ordre s'écrivent :

$$
\frac{\delta S}{\delta b} (\hat{b}) = -2 X'Y + 2 X'X \hat{b} = 0
$$

Les conditions du second ordre sont :

$$
\frac{\delta^2 S}{\delta b \delta b'} (\hat{b}) = 2 X'X
$$

et nous assurant que l'on obtient un minimum

d'où :

$$
\hat{X}'X \hat{b} = X'Y
$$

Deux cas se présentent alors :

• X'X inversible

• X'X non inversible

Dans le premier cas on obtient le résultat classique :

$$
\hat{b} = (X'X)^{-1} X'Y
$$

Dans le second (cas de multicollinéarité, variables non orthogonales...) nous
devons introduire la notion d'inverse généralisée.

D'après :

Soit A une matrice (m,n) de rang p. On appelle inverse généralisée de A
une matrice notée A⁻, de format (n,m) telle que :

$$
A \bar{A}^{-} A = A
$$
Nous allons montrer qu'une telle matrice existe. Considérons la matrice $A A'$ : elle est carrée $(m,m)$, symétrique de rang $p$ et définie non négative. On peut donc la diagonaliser et comme elle admet $p$ valeurs propres positives et $(m-p)$ valeurs propres nulles, il existe une matrice $P$ orthogonale telle que :

$$ P A A' P' = \begin{bmatrix} \Lambda_P & 0 \\ 0 & 0 \end{bmatrix} $$

$$ \text{avec } \mathbf{P}' = \mathbf{P}^{-1} \quad \text{et} \quad \Lambda_P = \begin{bmatrix} \lambda_1 & & & 0 \\ & \ddots & & \\ & & \lambda_p & 0 \end{bmatrix} $$

Partitionnons alors $P$

$$ P = \begin{bmatrix} P_1 \\ P_2 \end{bmatrix} \qquad \begin{array}{l} P_1 \text{ de format } (p,m) \\ P_2 \text{ de format } (m-p,m) \end{array} $$

Par suite de l'orthogonalité de $P$ on a

$$ PP' = \begin{bmatrix} P_1 P'_1 & P_1 P'_2 \\ P_2 P'_1 & P_2 P'_2 \end{bmatrix} = I = \begin{bmatrix} I & 0 \\ 0 & I \end{bmatrix} $$

ce qui implique

$$ P_1 P'_1 = I $$

Par ailleurs on a :

$$ PAA'P' = \begin{bmatrix} P_1 AA' P'_1 & P_1 AA' P'_2 \\ P_2 AA' P'_1 & P_2 AA' P'_2 \end{bmatrix} = \begin{bmatrix} \Lambda_P & 0 \\ 0 & 0 \end{bmatrix} $$

d'où

$$ P_1 AA' P'_1 = \Lambda_P $$

$$ \Rightarrow P'_1 P_1 AA' P'_1 = P'_1 \Lambda_P \text{ et puisque } P'_1 P_1 = I $$

$$ AA'P'_1 = P'_1 \Lambda_P $$
Comme

$$P'_{1}P_{1}A = A, \text{ si nous définissons la matrice G par } G = A'P'_{1} \Lambda^{-1}_{p} P_{1}$$

alors G est une inverse généralisée de A

En effet :

$$ 
\begin{align*}
AGA &= AA' P'_{1} \Lambda^{-1}_{p} P_{1} A = P'_{1} \Lambda_{p} \Lambda^{-1}_{p} P_{1} A \\
    &= P'_{1} P_{1} A \\
    &= A
\end{align*}
$$

En conclusion G existe toujours. Montrons qu'elle n'est pas l'unique inverse généralisée de A.

## Théorème

Soit G une inverse généralisée particulière de A.

Alors $\forall B$ matrice $(n,m)$ la matrice $G + B - GABAG$ est aussi une inverse généralisée de A,

## démonstration

Sachant que $AGA = A$ il suffit de vérifier que

$$A(G+B - GABAG)A = A$$

En effet en développant nous obtenons :

$$AGA + ABA - AGABAGA = A + ABA - ABA = A$$

## Théorème

Pour tout autre inverse généralisée $H$ de $A$ il existe une matrice $B$ $(n,m)$ telle que $H = G + B - GABAG$

## démonstration

Sachant que $AIBA = A$, choisissons $B = H - G$ et vérifions que

$$ 
\begin{align*}
G+B-GABAG &= G + H - G - GAHAG + GAGAG \\
           &= H - GAG + GAG \\
           &= H
\end{align*}
$$
Donc il existe toujours une et en général une infinité d'inverses généralisées de A.

## Application aux systèmes linéaires

$$ \text{Soit le système } A \begin{bmatrix} b \\ (m,n) \end{bmatrix} = g \begin{bmatrix} A \\ (n,1) \end{bmatrix} \quad (m,1) $$

$$ \text{avec rang } A = p < \text{Min}(m,n) $$

## Définition

On dit que le système ci-dessus est compatible si pour A et g donnés, il admet au moins une solutions en $\hat{b}$

### Théorème

$$ A \hat{b} = g \text{ compatible} \iff AA^{-} g = g \quad \forall A \neq 0 $$

### démonstration

$\Rightarrow$ S'il existe une solution $\hat{b}_0 / A \hat{b}_0 = g$ alors pour toute inverse généralisée $A^{-1}$ on a

$$ AA^{-}g = AA^{-}A \hat{b}_{0} = A \hat{b}_{0} = g $$

$\Leftarrow$ Soit G une inverse généralisée particulière de A telle que $AGg = g$. Alors $Gg$ est une solution du système et celui-ci est compatible.

### Théorème

Si le système $A \hat{b} = g$ est compatible il admet pour solution générale

$$ \hat{b} = A^{-1}g + (I - A^{-1}A)h \quad \forall A \neq 0 \quad \forall h \text{ vecteur de } \mathbb{R}^n $$

### démonstration

Il s'agit de montrer que si $AA^{-}g = g$ alors

1) $\forall A^{-1}$ et $h$ $\hat{b} = A^{-1}g + (I - A^{-1}A)h$ est tel que $A \hat{b} = g$, ce qui est inmédiat puisque

$$ A \hat{b} = AA^{-}g = g $$
2) $\hat{b}_0$ tel que $A \hat{b}_0 = g, \exists A^{-} \text{ et h tels que } \hat{b}_0 = A^{-}g + (I - A^{-}A) h$
Il suffit pour cela de prendre $h = \hat{b}_0$ et n'importe quelle matrice $A^{-}$.

Ainsi dans le cadre de la régression nous obtenons si

$$A = X'X, g = X'Y \\ (m,m) \quad (m,1)$$

$$\hat{b} = (X'X)^{-} X'Y + (I - (X'X)^{-} (X'X)) h$$

Pour évaluer $(X'X)^{-}$ nous utilisons la méthode due à GAUSS-DOOLITTLE
(cf.: Annexe 1).

Nous rappelons ci-dessous quelques résultats et tests statistiques
permettant de juger la validité de la régression.

- Dans le cas où le modèle contient un terme constant (une des variables explicatives constante et égale à 1 sur tout l'échantillon) on mesure l'adéquation du modèle à l'échantillon par le coefficient de détermination $R^2$

$$R^2 = 1 - \frac{\text{variance residuelle}}{\text{variance totale}} = 1 - \frac{\hat{\epsilon}' \hat{\epsilon}}{(y - \hat{y})' (y - \hat{y})}$$

$$\text{où } \hat{\epsilon} = \hat{y} - y, e = \begin{bmatrix} 1 \\ \vdots \\ 1 \end{bmatrix} \text{ vecteur } (N,1), \bar{y} = \frac{1}{N} \sum_{n=1}^{N} y_n$$

Remarque :

Dans le cas d'un modèle sans constante (cas rectificatif puisqu'il impose
à l'hyperplan de régression de passer par l'origine) on mesure cette adequati
par

$$\frac{\hat{y}'\hat{y}}{y'y} = 1 - \frac{\hat{\epsilon}'\hat{\epsilon}}{y'y}.$$

- L'estimateur de $\sigma^2$ est donné par $\hat{\sigma}^2 = \frac{\hat{\epsilon}'\hat{\epsilon}}{N-m}$ (estimateur sans biais)
La matrice de covariances de $\delta$ est donnée par

$$V(\delta) = \sigma^2 (X'X)^{-1}_r (X'X) (X'X)^{-1}_r$$

où r est le rang de $(X'X)$

$$\text{Donc } \hat{V}(\delta) = \hat{\sigma}^2 (X'X)^{-1}_r (X'X) (X'X)^{-1}_r$$

L'évaluation du Durbin-Watson (test d'indépendance des erreurs) permet
de vérifier l'hypothèse

$$E \epsilon_n \epsilon_{n'} = \delta_{nn'} \sigma^2 \quad \forall n, n' = 1, N$$

$$d = \frac{\sum_{n=1}^{N-1} (\hat{\epsilon}_{n+1} - \hat{\epsilon}_n)^2}{\sum_{n=1}^{N} \hat{\epsilon}_n^2}$$

(cf Annexe 2)

Les résultats rappelés ci-dessus seront utilisés sur chacun des
sous-échantillons, donc après que le découpage en K classes ait été effectué.

## 2. L'algorithme

L'algorithme de régression typologique a été proposé par E. DIDAY (1976 [6]) Nous le reprenons ici dans un cadre plus général en introduisant l'inverse généralisée et la ridge-régression. Il est basé sur la construction de trois fonctions f,g et w.

La fonction f fournit la partition, la fonction g, les hyperplans, et W mesure la validité des résultats. Une itération de l'algorithme consiste en l'application des trois fonctions f, g et w successivement.
## 2.1. - Les notations employées

Les données se présentent sous la forme de N lignes et (m+1) colonnes notées ($X^1$ ... $X^m$, y) où le ième élément est $z_i = (X_{i1}^1, ..., X_{i1}^m)$, $\forall i = 1,N$

1. ($X^1$... $X^m$) sont les m variables explicatives

2. y est la variable à expliquer

## Définitions

On appellera espace des noyaux $L_K = \{ (L_1, ..., L_K) \} \text{ l'ensemble des vecteurs (L_1, ..., L_K) où}$

$$L_j = \begin{bmatrix} b_j^1 \\ \vdots \\ b_j^m \end{bmatrix} = B_j = \text{vecteur des coefficients de régression}$$

On a ainsi $L_K = R^{mK}$

On appellera espace des partitions $P_K = \{ (P_1, ..., P_K) \} \text{ l'ensemble des des partitions en K classes au plus des observations}$

1. $n_j = \text{Card}(P_j)$, $\forall j = 1,K$

2. $z_i^j$ = ième élément de $P_j$

On a alors :

$$Z^j = \begin{bmatrix} Z^j_1 \\ \vdots \\ Z^j_{n_j} \end{bmatrix}, \quad X^j = \begin{bmatrix} x^j_1 & \cdots & x^j_{m_j} \\ \vdots & \ddots & \vdots \\ x^j_{n_j} & \cdots & x^j_{m_j} \end{bmatrix}, \quad y^j = \begin{bmatrix} y^j_1 \\ \vdots \\ y^j_{n_j} \end{bmatrix}$$
## 2.2. - Le critère

L'algorithme consiste en la recherche simultanée d'une partition en K classes de l'ensemble des observations et des K hyperplans de régression associés de façon à minimiser le critère :

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

$$ 
\begin{align*} 
W (L,P) &= \sum_{j=1}^{K} \left| \left| y^j - X^j B_j \right| \right|^2 n_j \\ 
&= \sum_{j=1}^{K} \sum_{i=1}^{n_j} (y_i^j - X_i^j B_j)^2 \\ 
&= \sum_{j=1}^{K} \sum_{i=1}^{n_j} (y_i^j - \sum_{e=1}^{m} x_i^{ej} b_j^{e})^2 
\end{align*} 
$$

On minimise la somme des carrés des résidus associés à la partition en K classes des N observations. Il s'agit de trouver le couple $(L^*, P^*)$ tel que :

$$W (L^*, P^*) = \min_{\mathbb{L}_K \times \mathbb{P}_K} W (L,P)$$

## 2.3. - Les fonctions f et g

La fonction d'identification f est une application de $\mathbb{L}_K$ dans $\mathbb{P}_K$

$$f : \mathbb{L}_K \longrightarrow \mathbb{P}_K$$

$$f : L \longrightarrow P = (P_1, \dots, P_K)$$

### Définition

Soit $d : E \times \mathbb{R}^m \to \mathbb{R}^+$ une fonction de proximité entre un individu et un noyau définie par :

$$d(Z, B_i) = (y - X B_i)^2$$
$$où X = (x^1 ... x^m)$$

$$Z = (X,y)$$

E = ensemble des observations

On définit alors $P_i$ ($\forall i = 1,K$) par :

$$P_i = \{ Z | d(Z, B_i) < d(Z, B_j) \forall j \neq i \}$$

Dans le cas où $\exists j / d(Z,B_i) = d(Z, B_j)$ on affecte Z à $P_{i\wedge j}$, c'est à-dire à la classe de plus petit indice.

La fonction de représentation g est une application de $\mathbb{P}_K$ dans $\mathbb{L}_K$

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

$$g : P \to L = (L_1 ... L_K) = (B_1 ... B_K)$$

où $B_j$ minimise $R(B,j,P) = \sum_{ZeP_j} d(Z,B)$ et

$$R(B_j,j,P) = \inf \{ R(B_j,j,P) \mid B \in \mathbb{R}^m \}$$

d'où par les équations normales et en introduisant l'inverse généralisée :

$$B_j = B_j(h) = ((x^j)'(x^j))^{-} (x^j)' y^j + (L_m - ((x^j)'(x^j))^{-})(x^j)'(x^j)). h$$

où h est un vecteur quelconque de $\mathbb{R}^m$.

Afin d'éviter que la reconstitution de la variable $y^j$ ne se fasse à l'aide de termes de valeurs élevées (qui se retranchant donnent une reconstitution illusoire de $y^i$) nous choisirons pour $B^j$ le $B^j(h)$ de norme minimale (cf. § 3 et annexe 3).
## 2.4. - Les suites $u_n$ et $v_n$

La suite $v_n$ est définie par :

• $v^{(0)} = (L^{(0)}, P^{(0)})_{où P^{(0)} (ou L^{(0)})}$ est tirée au hasard ou estimée

• $v^{(n+1)}$ est déduit de $v^{(n)}$ par

$$L^{(n+1)} = g(P^{(n)}, \text{et } P^{(n+1)} = f(L^{(n)}, \dots)$$

La suite $u_n$ est définie par $u_n = W(v_n)$

## 2.5. - La convergence de l'algorithme

**Théorème** : La suite $u_n$ converge en décroissant

démonstration

$$\text{Montrons } [1] \quad W(L^{(n)}, P^{(n)}) \ge W(L^{(n+1)}, P^{(n)})$$

$$[2] \quad W(L^{(n+1)}, P^{(n)}) \ge W(L^{(n+1)}, P^{(n+1)})$$

$[1]$ et $[2] \Rightarrow u_n$ non croissante.

Pour $[1]$

$$W(L^{(n)}, P^{(n)}) = \sum_{j=1}^{K} \sum_{i=1}^{n_j} (y_i^j - x_i^j)^2$$

$$W(L^{(n+1)}, P^{(n)}) = \sum_{j=1}^{K} \sum_{i=1}^{n_j} (y_i^j - x_i^j)^2$$

or $B_j^{(n+1)}$ est $\forall j = 1, K$ l'élément qui minimise

$$R(B_j, P^{(n)}) = \sum_{z \in P_j^{(n)}} d(z, B) \text{ avec } z = (x, y)$$

donc on a $[1]$
Pour [2]

$$
\text{On a } W(L^{(n+1)}, P^{(n+1)}) = \sum_{j=1}^{K} \sum_{i=1}^{n_j'} (y_i^j - x_i^j) (B_j^{(n+1)})^2 \\
\text{où } |P_j^{(n+1)}| = n_j'
$$

Dans $W(L^{(n+1)}, P^{(n)})$ et $W(L^{(n+1)}$, $P^{(n+1)}$) on ne retrouve chaque élément qu'une fois et une seule.

Soit i (X,y) l'indice de la classe à laquelle appartient (X,y) dans P<sup>(n)</sup>.
Soit i' (X,y) ce même indice dans P<sup>(n+1)</sup>.

On a alors

$$
W(L^{(n+1)}, P^{(n)}) = \sum_{z \in E} d(z, B_i^{(n+1)}) \quad \text{ou} \quad z = (X, y)
$$

$$
W(L^{(n+1)}, P^{(n+1)}) = \sum_{z \in E} d(z, B_i^{(n+1)}(X_j, y))
$$

Par construction de P<sup>(n+1)</sup> on a

$$
d(z, B_i^{(n+1)}(X_j, y)) \ge d(z, B_i^{(n+1)}(X_j, y)), \forall z \in E
$$

donc en sommant on a [2]

De [1] et [2] on a $u_n$ non croissante. Or $u_n$ est minorée par 0, donc $u_n$ est convergente.

La convergence s'effectue en un nombre fini d'étapes puisque $|P_k|$ est fini et que les éléments de la suite $L^{(n)}$ appartiennent à $g(P_k)$ avec $g : L_k \to g(P_k)$ fini.

Remarque : La suite $(L^{(n)}, P^{(n)})$ converge vers un élément non biaisé c'est-à-dire que $(L^{(n)}, P^{(n)}) \to (L, P) = (gof(L), f(L))$

L'indication

On sait que $u_n$ converge et atteint sa limite $\overline{3N/4n} > N$

$$
W(L^{(n)}, P^{(n)}) = W(L^{(N)}, P^{(N)})
$$
Mais on a $W(L^{(n+1)}, P^{(n+1)}) \le W(L^{(n+1)}, P^{(n)}) \le W(L^{(n)}, P^{(n)})$
d'où $W(L^{(n+1)}, P^{(n+1)}) = W(L^{(n+1)}, P^{(n)}) = W(L^{(n)}, P^{(n)})$

La dernière égalité s'écrit

$$ \sum_{j=1}^{k} \sum_{i=1}^{n_j} (y_i^j - x_i^j)^2 = \sum_{j=1}^{k} \sum_{i=1}^{n_j} (y_i^j - x_i^j)^2 \quad (n_j) $$

Par construction des $B_j^{(n+1)}$ on a

$$ \sum_{i=1}^{n_j} (y_i^j - x_i^j)^2 \le \sum_{i=1}^{n_j} (y_i^j - x_i^j)^2 \quad \forall j = 1, K $$

L'égalité n'est possible que si

$$ \sum_{i=1}^{n_j} (y_i^j - x_i^j)^2 = \sum_{i=1}^{n_j} (y_i^j - x_i^j)^2 \quad \forall j = 1, K \\ \Rightarrow \quad B_j^{(n+1)} = B_j^{(n)} \quad \forall j = 1, K \quad \Rightarrow \quad L^{(n+1)} = L^{(n)} $$

$$ \text{Comme } P^{(n+1)} = f(L^{(n+1)}) = f(L^{(n)}) \Rightarrow P^{(n)} = P^{(n+1)} $$

$$ \text{donc } \exists N \in \mathbb{N} \text{ tel } (L^{(n)}, P^{(n)}) = (L^{(N)}, P^{(N)}) $$

La suite $(L^{(n)}, P^{(n)})$ est une suite convergente. Soit $(L,P)$ la limite. On a
$(L,P) = (L^{(n+1)}, P^{(n+1)}) \quad \forall n > N$

$$ \begin{align*} \text{De plus } P^{(n)} &= f(L^{(n)}) \Rightarrow P = f(L) \\ \text{et } L^{(n+1)} &= g(P^{(n)}) \Rightarrow L = g(P) = \text{gof } (L) \\ \Rightarrow \quad (L,P) &= (\text{gof } (L), f(L)) \Rightarrow (L,P) \text{ non biaisé} \end{align*} $$

### Remarque :

L'utilisation des matrices inverses généralisées est particulièrement intéressante dans le cas où le nombre d'observations $n_j$ d'une partie $P_j$ est faible relativement au nombre de variables explicatives. En effet dans ce cas l'estimation par la méthode classique serait de mauvaise qualité, car dépendant fortement de l'échantillon. L'utilisation est également utile dans le cas de variables explicatives non orthogonales.
### 3. - Cas de variables non orthogonales : la ridge régression

Dans le cas de variables explicatives non orthogonales, la matrice de corrélations est "éloignée" de la matrice identité.

L'estimateur des moindres carrés tend à être de norme élevée et ses coefficients peuvent avoir un mauvais signe.

En effet soit le modèle :

$$Y = Xb + \varepsilon \quad E \varepsilon = 0, V\varepsilon = \sigma^2 I$$

$$On a \hat{b} = (X'X)^{-1} X'Y$$

Calculons la distance de $\hat{b}$ à $b$

$$E (\hat{b}-b)' (\hat{b}-b) = E (\text{trace} (\hat{b}-b)' (\hat{b}-b))$$

$$\text{puisque } (\hat{b}-b)' (\hat{b}-b) \in \mathbb{R}$$

$$\begin{cases} \text{Or } \hat{b} = (X'X)^{-1} X'Y \\ Y = Xb + \varepsilon \end{cases} \quad \Rightarrow \quad \begin{cases} \hat{b} = (X'X)^{-1} X' (Xb + \varepsilon) \\ \hat{b} = b + (X'X)^{-1} X' \varepsilon \end{cases}$$

d'où

$$\begin{aligned} (\hat{b}-b)' (\hat{b}-b) &= ((X'X)^{-1} X'\varepsilon)' ((X'X)^{-1} X' \varepsilon) \\ &= \text{trace}((X'X)^{-1} X'\varepsilon)' ((X'X)^{-1} X'\varepsilon) \\ &= \text{trace}((X'X)^{-1} X'\varepsilon)(\varepsilon'X (X'X)^{-1}) \end{aligned}$$

donc

$$\begin{aligned} E (\hat{b}-b)' (\hat{b}-b) &= \text{trace} \left((X'X)^{-1} X' \right) \sigma^2 (X(X'X)^{-1}) \\ &= \sigma^2 \text{trace} (X'X)^{-1} \end{aligned}$$

Soient $(\lambda_i)_{i=1,m}$ les m valeurs propres de $X'X$, donc $(1/\lambda_i)_{i=1,m}$

celles de $(X'X)^{-1}$. On a alors :

$$E (\hat{b}-b)' (\hat{b}-b) = \sigma^2 \sum_{i=1}^{m} \frac{1}{\lambda_i}$$
donc si $\lambda_{\max} = \lambda_1 > \lambda_2 \dots > \lambda_m = \lambda_{\min} > 0$

$$ E (\hat{b}-b)' (\hat{b}-b) > \sigma^2 \cdot \frac{1}{\lambda_m} $$

Si la matrice de corrélations est éloignée de la matrice identité alors $\lambda_m$ est voisin de zéro, donc la distance de $\hat{b}$ à $b$ est grande, ce qui implique une mauvaise estimation.

Par contre si elle est proche de la matrice identité alors cette distance est voisine de $m.\sigma^2$

Pour résoudre ce problème nous allons introduire la ridge régression.

Dans la ridge régression on impose au vecteur $\hat{b}$ d'avoir une norme pas trop grande de façon à éviter une fausse reconstitution de la variable à expliquer.

$$ \text{On cherche } \hat{b} / || Y - X \hat{b} ||^2 = \text{Min } || Y - X b ||^2 \\ \{ b / ||b||^2 = \alpha \} $$

où $\alpha$ est fixé a priori.

Par la méthode de Lagrange ce problème est équivalent à minimiser

$$ L (b, \lambda) = || Y - X b ||^2 + \lambda ( || b ||^2 - \alpha) $$

où $\lambda$ est un multiplicateur de Lagrange.

Par dérivation nous obtenons

$$ - 2 X'Y + 2 X'X \hat{b} + 2 \lambda \hat{b} = 0 $$

$$ \text{d'où } (X'X + \lambda I) \hat{b} = X'Y $$

$$ \text{et } \hat{b} = [(X'X) + \lambda I]^{-1} X'Y = \hat{b}(\lambda) $$

$$ \text{où } \lambda \text{ est choisi tel que } || \hat{b} ||^2 = \alpha $$

Pour déterminer la valeur de $\lambda$ on peut utiliser un processus itératif qui pour différentes valeurs de $\lambda$ ($\lambda = 0$ ; $0,1$ ; $0,2$ ; ...) évalue $\hat{b}$. On s'arrête dès que les valeurs des coefficients se sont stabilisées. Cette procédure est utile lorsqu'on ne connaît pas a priori la valeur à donner à $\alpha$.
### 3.1. - Comparaisons entre inverses généralisées et ridge régression

Ces deux méthodes sont supérieures à la méthode classique des moindres carrés.

D'après MARQUARDT, les inverses généralisées sont particulièrement intéressantes dans le cas où une ou plusieurs valeurs propres de la matrice de corrélations sont nulles ; la ridge régression étant plus adaptée au cas de valeurs propres voisines de zéro.

Nous allons introduire une nouvelle méthode de calcul de l'inverse généralisée, en utilisant le calcul des valeurs propres. Soit $A = X'X$. Soient $D$ la matrice diagonale des valeurs propres ordonnées $\lambda_1 > \lambda_2 > \dots > \lambda_m$ ; $S$ la matrice des vecteurs propres associés. On a $S'A S = D$ où $S'S = I$ d'où $A^{-1} = S D^{-1} S'$ si $A$ inversible.

Supposons $A$ de rang $r$, de telle sorte que les $(m-r)$ derniers éléments $d_i$ soient égaux à zéro (ou voisins).

On a ainsi $S = (S_r, S_{m-r})$ où $S_r$ est $(m,r)$

$S_{m-r}$ est $(m,m-r)$

$$ D = \begin{bmatrix} D_r & 0 \\ 0 & D_{m-r} \end{bmatrix} \quad \text{où } D_r \text{ est } (r,r) \\ D_{m-r} \text{ est } (m-r, m-r) $$

Si $D_{m-r}^{-1} = 0$ alors par définition $D_{m-r}^{-1} = 0$

$$ d'où A_r^{-} = S_r D_r^{-1} S'_r \quad \text{ou} \quad A_r^{-} = \sum_{i=1}^{r} \frac{1}{\lambda_i} S_i S'_i $$

où $S_i$ est le vecteur propre de $A$ associé à $\lambda_i$.
Pour choisir r on peut utiliser le critère suivant :

on teste $\frac{\lambda_m}{\text{Trace D}}$. Si $\frac{\lambda_m}{\text{Trace D}} > \omega$ où $\omega = 10^{-5}$ par exemple on admet $r = m$

sinon on teste $\frac{\lambda_m + \lambda_{m-1}}{\text{Trace D}}$ par rapport à $\omega$

Si $\frac{\lambda_m + \lambda_{m-1}}{\text{Trace D}} > \omega \Rightarrow r = m-1$. Sinon on itère.

A l'ordre k si $\frac{\sum_{i=m}^{m-k} \lambda_i}{\text{trace D}}$ est supérieur à $\omega$ alors $r = m-k$

sinon on passe à l'ordre $(k+1)$.

**Justification de cette procédure** : sachant que la trace est preservée par transformations orthogonales, et que $\lambda_j$ est la somme des carrés des projections des points représentés par les lignes de X sur le vecteur propre associé, le critère ci-dessus est tel que le rang choisi pour la matrice permet de représenter "à peu près" toutes les variations des points de X.

Suivant MARQUARDT nous allons introduire un estimateur combinant les avantages associés à l'inverse généralisée et à la ridge régression.

Soit $b = ((X'X)_r + \lambda I)^{-1}$ X'Y où r est choisi de façon à éliminer les valeurs propres nulles (ou voisines aux arrondis près) et $\lambda$ est tel qu'il réduit l'effet des valeurs propres proches de zéro.

$$ (avec (X'X)_r = S_r D_r S_r' = \sum_{j=1}^{r} \lambda_j S_j S_j') $$

**Remarque** La convergence de l'algorithme n'est pas modifiée. En effet le vecteur des coefficients de régression est obtenu par la méthode des moindres carrés sous contraintes.
### 3.2. - Schéma général de l'algorithme

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

Calcul de $L^{(n)}$

<figure><img src="image_3.png" /></figure>
### 3.3. - Programme informatique REGTYPO

L'algorithme précédemment décrit est programmé en Fortran.
L'organigramme est semblable à 3.2. Il comporte 3 modes de calcul des
estimateurs :

* MODE = 0 inversion classique

* MODE = 1 inversion généralisée et ridge regression
  (applicable aux modèles sans constante)

* MODE = 2 inversion généralisée

Il nécessite en entrée :

N : nombre d'observations

K : nombre de classes

m : nombre de variables explicatives

ITER : nombre maximum d'itérations de l'algorithme

(L<sub>1</sub><sup>(0)</sup> ... L<sub>K</sub><sup>(0)</sup>) : noyaux initiaux (ou P<sub>1</sub><sup>(0)</sup> ... P<sub>K</sub><sup>(0)</sup>)

Y,X : données exogènes et endogènes

En sortie nous obtenons :

* P = (P<sub>1</sub> ... P<sub>K</sub>), L = (L<sub>1</sub> ... L<sub>K</sub>), W

* les estimateurs $\hat{\sigma}_j^2$ : j = 1,K des variances

* les matrices de variances des paramètres

* les $R_j^2$ : j = 1,K carrés des coefficients de corrélation

* Les Durbin-Watson
  et un graphique des variables Y et ŶY
## 4. - Elimination de l'effet de "cylindre"

L'effet de cylindre remarqué en analyse factorielle typologique peut également s'introduire en régression typologique lorsque le tirage aléatoire de la partition initiale (ou des noyaux initiaux) est "malheureux", comme le montre l'exemple suivant :

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

Dans ce cas ($m = 2, K = 2$) le tirage aléatoire de $P = (P_1, P_2)$ ne permettra pas de reconnaître la structure en deux classes homogènes de la population :

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

Il s'agit donc de modifier l'algorithme au niveau de ses premières itérations afin d'éliminer cet effet de cylindre.

### 4.1. - Ellipsoïde d'un vecteur aléatoire

$$ \text{Soient } Z_i^j = (Y_i^j, X_i^j, \dots, X_i^{jm}) \quad \forall j = 1, \dots, K \\ \phantom{\text{Soient } Z_i^j = (} \forall i = 1, \dots, n_j $$

les vecteurs à $(m+1)$ dimensions définissant les individus.
Sur $P_j$ nous avons ainsi $n_j$ réalisations d'un vecteur $\overline{Z_j}^j$ à (m+1) dimensions

Soient $\overline{Z}^j$ moyenne de $Z_i^j$

$V_j$ matrice des covariances des $Z_i^j$ de format (m+1, m+1) de rang r

## Définition

Si $r = m+1$ on appelle ellipsoïde du vecteur à (m+1) dimensions l'ensemble des vecteurs $W$ à (m+1) dimensions tels que :

$W' V_j^{-1} W \le C_j^2$ avec $C_j$ constante positive

Si $C_j = 1$ alors on obtient l'ellipsoïde indicateur.

Si $r < m+1$ alors il existe des vecteurs $u$ non nuls tels que $V_j u = 0$ d'où $u' V_j u = 0$ donc par définition de $V_j$ :

$$ 
\begin{aligned}
& E \left[ u' (\overline{Z^j} - \overline{Z}^j)^2 \right] = 0 \\
\Rightarrow \quad & u' (\overline{Z^j} - \overline{Z}^j) = 0
\end{aligned}
$$

donc le vecteur $(\overline{Z^j} - \overline{Z}^j)$ appartient au sous espace vectoriel orthogonal à tous les vecteurs $u$ vérifiant $V_j u = 0$.

On appelle support de $\overline{Z}^j - \overline{Z}^j$ cet espace de dimension r

Soit $R_j$ la matrice (m+1, m-r+1) de rang (m-r+1) telle que $V_j R_j = 0$, définissant ce support. Ainsi $R_j' Z_j = 0$

## Définition

Si $r < m+1$ alors l'ellipsoïde est l'ensemble des vecteurs $W$ tels que

$$ 
\begin{gathered}
W' (V_j + R_j R_j')^{-1} W \le C_j^2 \\
\text{vérifiant } R_j' W = 0
\end{gathered}
$$
## 4.2. - Les modifications apportées à l'algorithme

### Définition

On appelle partie indicatrice de $P_j$ l'ensemble des éléments de $P_j$ intérieurs à l'ellipsoïde indicateur.

Soit $Q_j$ cet ensemble défini par :

$$ egin{cases} 
\text{ze } P_j \text{ (ze } Q_j \leftarrow \left( Z - \bar{Z}_j \right)' V_j^{-1} (Z - \bar{Z}_j) \le 1) & \text{si } r=m+1 \\ 
Z \in P_j \text{ ( Z \in Q_j \leftarrow \left( Z - \bar{Z}_j \right)' (V_j + R_j R_j')^{-1} (Z - \bar{Z}_j) \le 1) } & \text{si } r<m+1 
\end{cases} $$

$$ \text{Ainsi, } Q_j \subsetneq P_j \quad \forall j = 1, K $$

($Q_1$ ... $Q_K$) n'est pas une partition de $E$

Pour les premières itérations l'algorithme aura la démarche suivante :

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

avec :

$$ \bullet \ h : \mathbb{P}_K \longrightarrow (\mathcal{S}(E))^K $$

$$ \bullet \ P_j \longrightarrow Q_j = \{ Z \in P_j \text{ intérieur à son ellipsoïde indicatrice. } \} $$

$$ \bullet \ g' : (\mathcal{S}(E))^K \longrightarrow \mathbb{P}_K \text{ extension de } g \text{ à } (\mathcal{S}(E))^K $$

$$ \bullet \ f : \mathbb{P}_K \longrightarrow \mathbb{P}_K \text{ définie comme précédemment} $$
La convergence de l'algorithme n'est plus assurée lorsqu'on applique cette procédure (ceci est dû à la fonction h).

Toutefois n'étant appliquée que pendant quelques itérations la convergence est de nouveau acquise par la suite.

### 4.3. - Taille des ellipsoïdes

Si on suppose que le vecteur $\bar{z}^j$ à (m+1) composantes suit une loi de probabilité normale de moyenne $\bar{z}_j$ et de matrice de covariances $V_j$ non singulière alors :

$$ (\bar{z}^j - \bar{z}_j)' V_j^{-1} (\bar{z}^j - \bar{z}_j) \rightsquigarrow \chi_{m+1}^2 $$

En effet si $V_j$ inversible, il existe une matrice D telle que $V_j = DD'$. Alors $\alpha = D^{-1}(\bar{z}^j - \bar{z}_j)$ suit une loi normale de matrice de covariances

$$ (D^{-1}) V_j (D^{-1})' = I $$

La somme des carrés de ses composantes suit un $\chi_{m+1}^2$ et s'écrit :

$$ \begin{aligned} \alpha' \alpha &= (\bar{z}^j - \bar{z}_j)' (D^{-1})' (D^{-1}) (\bar{z}^j - \bar{z}_j) \\ &= (\bar{z}^j - \bar{z}_j)' V_j^{-1} (\bar{z}^j - \bar{z}_j) \end{aligned} $$

On peut ainsi calculer en fonction de (m+1) et du pourcentage $\tau$ de la population désirée à l'intérieur de l'ellipsoïde la valeur de $\chi_{m+1}^2$ telle que $P(\chi^2 > \chi_{m+1}^2) = 1 - \tau$ et ainsi fixer $C_j$.
4.4. - Exemple d'application m = 2, K = 2

Soit le nuage suivant :

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

P₁⁽⁰⁾, P₂⁽⁰⁾ classes initiales

Les hyperplans associés sont :

H 1 y = -0,057 x + 16,5

H 2 y = -0,05 x + 13,4
Calculons les ellipsoïdes indicateurs et les parties indicatrices.

$$ \text{P 1} \left| \begin{array}{l} V(x) = 32,06 \\ \text{Cov}(xy) = -1,8375 \\ V(y) = 0,822 \end{array} \right. $$

$$ \text{P 2} \left| \begin{array}{l} V(x) = 59,16 \\ \text{Cov}(xy) = -2,84 \\ V(y) = 0,96 \end{array} \right. $$

d'où

$$ V_1^{-1} = \begin{bmatrix} 0,0537 & 0,08 \\ 0,08 & 1,395 \end{bmatrix} $$

$$ V_2^{-1} = \begin{bmatrix} 0,02 & 0,058 \\ 0,058 & 1,214 \end{bmatrix} $$

On teste alors si pour $j = 1$ ou $2$

$$ (x - \bar{x}, y - \bar{y}) \cdot V_j^{-1} \begin{pmatrix} x - \bar{x} \\ y - \bar{y} \end{pmatrix} \le 1 $$

d'où

$$ Q_1 = \{5,6\} \quad Q_2 = \{12, 13\} $$

qui donneront les hyperplans $H'_1$ et $H'_2$ passant par les points (5,6) et (12, 13) respectivement.

Autour de ces points s'agrégeront les classes

$$ P'_1 = \{1,2,3..., 10\} $$

$$ P'_2 = \{11,12, ..., 20\} $$

ainsi reconnues pour lesquelles on évaluera les hyperplans associés.
## 5. Choix optimal du nombre de classes

Quelle que soit la méthode de classification utilisée le problème du choix du nombre de classes se pose. En pratique on calcule un (ou plusieurs) indice permettant d'évaluer ce nombre de classes.

Nous en présenterons deux. Le premier est:

$$c(K) = \frac{\text{variance résiduelle inter-classes}}{\text{variance résiduelle totale}}$$

où K est le nombre de classes.

L'inertie résiduelle totale s'écrit:

$$\sum_{j=1}^{K} \sum_{i \in P_j} (\hat{Y}_j^i - \hat{Y}_j)^2$$

où $\hat{Y}_j^i$ est la valeur calculée de $\hat{Y}_j^i$ par une régression globale.

Soit $\hat{Y}_j^{iL}$ la valeur calculée de $\hat{Y}_j^i$ par la méthode de régression typologique en K classes. Montrons que l'on a:

$$\sum_{j=1}^{K} \sum_{i \in P_j} (\hat{Y}_j^i - \hat{Y}_j^{iL})^2 = \sum_{j=1}^{K} \sum_{i \in P_j} (\hat{Y}_j^i - \hat{Y}_j^{iL})^2 + \sum_{j=1}^{K} \sum_{i \in P_j} (\hat{Y}_j^{iL} - \hat{Y}_j^{i})^2$$

c'est à dire :

$$\text{variance résiduelle totale} = \text{variance résiduelle intra-classes} + \text{variance résiduelle inter-classes}$$

Nous avons $\forall j, \forall i$

$$\hat{Y}_j^i - \hat{Y}_j^{iL} = (\hat{Y}_j^i - \hat{Y}_j^{iL}) + (\hat{Y}_j^{iL} - \hat{Y}_j^i)$$

d'où sous forme vectorielle :

$$\|\hat{Y} - \hat{Y}^L\|^2 = \|\hat{Y} - \hat{Y}^L\|^2 + \|\hat{Y}^L - \hat{Y}\|^2 + 2 \langle \hat{Y} - \hat{Y}^L, \hat{Y}^L - \hat{Y} \rangle$$

Il s'agit de montrer que ce dernier terme est nul. Or pour tout j nous avons:

$$\langle \hat{Y}_j - \hat{Y}_j^L, \hat{Y}_j - \hat{Y}_j^L \rangle = 0$$

car $\hat{Y}_j^L$ se projette en $\hat{Y}_j^L$ sur la variété linéaire engendrée par les variables explicatives dans chaque classe $P_j$.

$$\text{Il reste alors à montrer } \langle \hat{Y} - \hat{Y}^L, \hat{Y} \rangle = 0$$
<figure><img src="image_8.png" /></figure>

**Proposition**

$$ \left( \sum_{i=1}^{K} A_i \right)^{-1} = \left( K I + \sum_{i=1}^{K} \sum_{j \neq i}^{K} A_i^{-1} A_j \right)^{-1} \left( \sum_{i=1}^{K} A_i^{-1} \right) $$

**démonstration:** En effet soit M la matrice inversible telle que:

$$ \left( \sum_{i=1}^{K} A_i \right)^{-1} = M \left( \sum_{i=1}^{K} A_i^{-1} \right) $$

On a :

$$ \begin{aligned} I &= M \left( \sum_{i=1}^{K} A_i^{-1} \right) \left( \sum_{i=1}^{K} A_i \right) \\ &= M \left( K I + \sum_{i=1}^{K} \sum_{\substack{j=1 \\ j \neq i}}^{K} A_i^{-1} A_j \right) \end{aligned} $$

$$ \text{d'où } M = \left( K I + \sum_{i=1}^{K} \sum_{\substack{j=1 \\ j \neq i}}^{K} A_i^{-1} A_j \right)^{-1} $$

$$ \text{et } \left( \sum_{i=1}^{K} A_i \right)^{-1} = \left( M \right) \left( \sum_{i=1}^{K} A_i^{-1} \right) $$

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

**Proposition:** Le vecteur des écarts résiduels locaux est orthogonal au vecteur des valeurs calculées globales.

$$ \langle Y - \hat{Y}^L, \hat{Y} \rangle = 0 $$

**démonstration:**
Soit $\mathbf{X} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_K \end{bmatrix}$ et $\mathbf{Y} = \begin{bmatrix} \mathbf{Y}_1 \\ \mathbf{Y}_2 \\ \vdots \\ \mathbf{Y}_K \end{bmatrix}$ les vecteurs des valeurs observées sur $P_1 \cdots P_K$.

On a : $\hat{\mathbf{Y}} = \mathbf{X} ( \mathbf{X}'\mathbf{X} )^{-1} \mathbf{X}'\mathbf{Y}$ et $\hat{\mathbf{Y}}_j^L = \mathbf{X}_j ( \mathbf{X}_j' \mathbf{X}_j )^{-1} \mathbf{X}_j' \mathbf{Y}_j \quad \forall j=1,K$

d'où $\mathbf{Y}_j - \hat{\mathbf{Y}}_j^L = ( \mathbf{I} - \mathbf{X}_j ( \mathbf{X}_j' \mathbf{X}_j )^{-1} \mathbf{X}_j' ) \mathbf{Y}_j \quad \forall j=1,K$

Décomposons le vecteur $\hat{\mathbf{Y}}$ en:

$$ \hat{\mathbf{Y}} = \begin{bmatrix} \hat{\mathbf{Y}}_1 \\ \hat{\mathbf{Y}}_2 \\ \vdots \\ \hat{\mathbf{Y}}_K \end{bmatrix}, $$

$$ \hat{\mathbf{Y}} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_K \end{bmatrix} ( \mathbf{K} \mathbf{I} + \sum_{i=1}^{K} \sum_{\substack{j=1 \\ j \neq i}}^{K} (x_i' x_i)^{-1} (x_j' x_j)^{-1} (\sum_{i=1}^{K} (x_i' x_i)^{-1}) \cdot \begin{bmatrix} x_1' \\ x_2' \\ \vdots \\ x_K' \end{bmatrix} \begin{bmatrix} \mathbf{Y}_1 \\ \mathbf{Y}_2 \\ \vdots \\ \mathbf{Y}_K \end{bmatrix} ) $$

avec $x'x = \sum_{i=1}^{K} x_i' x_i$ et la proposition précédente, d'où :

$$ \begin{bmatrix} \hat{\mathbf{Y}}_1 \\ \hat{\mathbf{Y}}_2 \\ \vdots \\ \hat{\mathbf{Y}}_K \end{bmatrix} = \begin{bmatrix} x_1 \mathbf{M} ( \sum_{i=1}^{K} (x_i' x_i)^{-1} ) ( \sum_{i=1}^{K} x_i' \mathbf{Y}_i ) \\ \vdots \\ x_K \mathbf{M} ( \sum_{i=1}^{K} (x_i' x_i)^{-1} ) ( \sum_{i=1}^{K} x_i' \mathbf{Y}_i ) \end{bmatrix} $$

avec $\mathbf{M} = ( \mathbf{K} \mathbf{I} + \sum_{i=1}^{K} \sum_{\substack{j=1 \\ j \neq i}}^{K} (x_i' x_i)^{-1} (x_j' x_j)^{-1} )^{-1}$

Ainsi avec $\langle \mathbf{Y} - \hat{\mathbf{Y}}^L, \hat{\mathbf{Y}} \rangle = \sum_{j=1}^{K} \langle \mathbf{Y}_j - \hat{\mathbf{Y}}_j^L, \hat{\mathbf{Y}}_j \rangle$
$$ \langle Y - \hat{Y}^L, \hat{Y} \rangle = \sum_{j=1}^{K} \langle (I - X_j(X_j'X_j)^{-1}X_j')Y_j, X_j M(\sum_{i=1}^{K} (X_i'X_i)^{-1}) (\sum_{i=1}^{K} X_i'Y_i) \rangle $$

$$ = \sum_{j=1}^{K} Y_j' (I - X_j(X_j'X_j)^{-1}X_j') X_j M(\sum_{i=1}^{K} (X_i'X_i)^{-1}) (\sum_{i=1}^{K} X_i'Y_i) $$

$$ = \sum_{j=1}^{K} Y_j' (X_j - X_j(X_j'X_j)^{-1}X_j'X_j) M(\sum_{i=1}^{K} (X_i'X_i)^{-1}) (\sum_{i=1}^{K} X_i'Y_i) $$

= 0.

On posera alors:

$$ c(K) = \frac{\| \hat{Y}^L - \hat{Y} \|^2}{\| Y - \hat{Y} \|^2} \quad \forall K \ge 1 $$

On a $c(K) = 0$ si $K = 1$

$c(K) = 1$ si la variance intra-classes est nulle,c'est à dire quand
le critère da l'algorithme de régression typologique est nul donc lorsque

$\forall j=1,K \mid P_j \mid \le m$ où m est le nombre de variables explicatives.

Pour évaluer K on trace la courbe $c(K)$ d'où par exemple:

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

En analysant cette courbe on peut conclure à une valeur optimale de K, où
serait pour notre exemple 4.

Le second indice est défini par analogie avec le F de Fisher

$$ F(K) = \frac{\| \hat{Y}^L - \hat{Y} \|^2}{\| Y - \hat{Y} \|^2} \quad \begin{array}{l} \text{N-K} \\ K-1 \end{array} $$

où N est le nombre d'observations.
# 6. Méthode de prévision associée aux modèles linéaires locaux

## 6.1 Le problème

Nous désirons résoudre le problème suivant: soit X un nouvel individu $X = (x_1^1, ..., x_m^m)$

On suppose que X appartient à l'un des groupes mais on ignore lequel. On cherche alors à définir une règle de décision permettant d'affecter cet individu supplémentaire à l'une des classes $P_1^*, ..., P_K^*$ connaissant uniquement X et à en déduire une valeur de la variable à expliquer Y pour cet individu.

## 6.2 La régression par boules

La méthode de régression par boules (BENZECRI, 1974; CAZES, 1975) traite dans le cas empirique le problème de la régression fonctionnelle.

Aux n réalisations $(Y_i, x_i^1, ..., x_i^m)_{i=1,n}$ on associe le nuage des n points de $IR^m$ (un point étant repéré par ses m coordonnées sur les variables explicatives). A chaque point est affecté une valeur $Y_i$. (Nous supposerons qu'il y a une probabilité nulle pour que deux points soient confondus).

Pour prédire la variable Y d'un nouvel individu repéré par $(x_1^1, ..., x_m^m)$ nous cherchons les s points (s fixé) du nuage des individus les plus proches de $(x_1^1, ..., x_m^m)$, la moyenne des valeurs de la variable à expliquer de ces s points nous donne la valeur de Y.

## 6.3 La règle d'affectation d'un nouvel individu

Soient $P_1^*, ..., P_K^*$ les K classes obtenues à la convergence de l'algorithme.

Soit $X = (x_1^1, ..., x_m^m)$ les m coordonnées d'un nouvel individu.

Soit $S = \{ s \in N / s_{Max} \ge s \ge s_{Min} \}$ à $s_{Max}$ et $s_{Min}$ fixés.

Soit $\forall s \in S, \forall j \in \{1, ..., K\} \rho_j(s)$ le nombre de ces points appartenant à $P_j^*$.
Soit $\forall s \in S \quad J_s = \{ j \in \{1, \dots, K\} : \ell_j(s) = \max_{\ell=1,K} \ell \in \{1\} \}$

Soit $\forall j \in \{1, \dots, K\} \ \exists j = \sum_{s=s_{\min}}^{s=s_{\max}} \chi_{J_s}^{(j)} \quad \text{où}$

$$ \chi_{J_s}^{(j)} = \begin{cases} 1 & \text{si } j \in J_s \\ 0 & \text{sinon} \end{cases} $$

$$ \text{Soit } J^* = \{ j^* : \exists j_* = \max_{j=1,K} \ell_j \} $$

Si $|J^*|=1$ alors nous affectons $X$ à $P_j^{*}$. Sinon nous modifions les bornes de variations de $s$. Ainsi nous obtenons la valeur prévue de $Y_i$:

$$ \hat{Y} = \sum_{i=1}^{m} b_{ji}^{i} X^{i} $$

6.4 **Un exemple** Nous supposerons $K=2, m=2, S_{\min} = 1, S_{\max} = 5$ $X = (5,2;1)$

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

Nous obtenons:

$J_1 = (1)$

$J_2 = (1,2)$

$J_3 = (2)$

$J_4 = (2)$

$J_5 = (2)$

d'où $X \in P_2^*$
## CONCLUSION

Dans [3] nous appliquons cet algorithme de régression typologique aux modèles des moindres carrés généralisés, des moindres valeurs absolues, de régression orthogonale, à erreurs sur les variables, autorégressifs.

D'autre part dans le cadre de la régression typologique nous abordons les problèmes de **reconnaissance des formes**. En utilisant les méthodes de régression sous contraintes linéaires locales et la classification non hiérarchique nous donnons des algorithmes de mémorisation condensée de courbes, (courbes de niveaux d'altitudes en cartographie automatique, par exemple). Ces algorithmes permettent à la reconstitution d'être continue, à dérivées continues ou spline. (Les raccords en pente et en courbure, dans le cas de courbes d'altitudes, semblent nécessaires à une bonne représentation du profil).

Enfin nous étendons les modèles de régression typologique (utilisables dans le cas de variables **quantitatives**) aux situations dans lesquelles nous désirons expliquer une variable **ordinale**, **qualitative** ou **quantitative** par un groupe de variable **hétérogènes** (mélange de variables ordinales, qualitatives, quantitatives). Nous donnons un algorithme qui de façon itérative fournit simultanément la partition des observations, les hyperplans locaux de régression et les codages associés aux variables.

Nous donnons également diverses applications de ces algorithmes notamment en géologie et en cartographie automatique.
# ANNEXES

## ANNEXE N° I

Nous donnons dans cette annexe une description, sur un exemple, de l'algorithme de GAUSS-DOLITTLE permettant de calculer l'inverse généralisée d'une matrice carrée symétrique. Soit à résoudre le système:

$$ (X'X) \hat{b} = X'Y $$

$$ \text{avec: } X'X = \begin{bmatrix} 1 & 2 & 3 \\ 2 & 4 & 6 \\ 3 & 6 & 12 \end{bmatrix} \qquad X'Y = \begin{bmatrix} 4 \\ 8 \\ 18 \end{bmatrix} $$

Nous avons le tableau:

<table><thead><tr><th>Ligne</th><th>I</th><th>2</th><th>3</th><th>0'</th><th>I'</th><th>2'</th><th>3'</th></tr></thead><tbody><tr><td>I</td><td>I</td><td>2</td><td>3</td><td>4</td><td>I</td><td>0</td><td>0</td></tr><tr><td>2</td><td>2</td><td>4</td><td>6</td><td>8</td><td>0</td><td>I</td><td>0</td></tr><tr><td>3</td><td>3</td><td>6</td><td>12</td><td>18</td><td>0</td><td>0</td><td>I</td></tr><tr><td>1.0</td><td>I</td><td>2</td><td>3</td><td>4</td><td>I</td><td>0</td><td>0</td></tr><tr><td>1.1</td><td>I</td><td>2</td><td>3</td><td>4</td><td>I</td><td>0</td><td>0</td></tr><tr><td>2.0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-2</td><td>I</td><td>0</td></tr><tr><td>2.1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>3.0</td><td>0</td><td>0</td><td>0</td><td>6</td><td>-3</td><td>0</td><td>I</td></tr><tr><td>3.1</td><td>0</td><td>0</td><td>I</td><td>2</td><td>-I</td><td>0</td><td>I/3</td></tr></tbody></table>

### pivots

Les lignes (i,0) sont les lignes non normalisées.

Lorsque le pivot i est nul, on rend nul toute la ligne (i,I)

On définit la matrice $(X'X)^{-}$ comme le produit de la matrice constituée des $(i,0)_i$ transposée par celle des $(i,I)_i$ d'où:

$$ (X'X)^{-} = \begin{bmatrix} I & -2 & -3 \\ 0 & I & 0 \\ 0 & 0 & I \end{bmatrix} \begin{bmatrix} I & 0 & 0 \\ 0 & 0 & 0 \\ -I & 0 & I/3 \end{bmatrix} = \begin{bmatrix} 4 & 0 & -I \\ 0 & 0 & 0 \\ -I & 0 & I/3 \end{bmatrix} $$

On obtient alors le vecteur b, solution de $X'X \hat{b} = X'Y$ comme le produit de la matrice $(i,I)_i$ transposée par le vecteur Y des rangées $(i,0)_i$.

$$ \hat{b} = \begin{bmatrix} I & 0 & -I \\ 0 & 0 & 0 \\ 0 & 0 & I/3 \end{bmatrix} \begin{bmatrix} 4 \\ 0 \\ 6 \end{bmatrix} = \begin{bmatrix} -2 \\ 0 \\ 2 \end{bmatrix} $$
ANNEXE N° 2

L'évaluation du DURBIN-WATSON (test d'indépendance des erreurs) permet de vérifier l'hypothèse:

$$E \epsilon_n \epsilon_{n'} = \hat{\delta}_{n n'} \sigma^2$$

Il s'écrit:

$$d = \frac{\sum_{n=1}^{N-1} ( \hat{\epsilon}_{n+1} - \hat{\epsilon}_n )^2}{\sum_{n=1}^{N} \hat{\epsilon}_n^2}$$

Nous savons que $d_L \le d \le d_U$. Soit $d^*$ la valeur calculée sur l'échantillon. Pour N observations sur m variables explicatives on compare $d^*$ à $d_L$ et $d_U$:

Soit $H_0$ l'hypothèse: il n'y a pas de corrélation $\rho=0$

$H_1$ l'hypothèse alternative $\rho \ne 0$

on a alors:

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

$AH_i$ zone d'acceptation de l'hypothèse i $i=0,1$

? zone d'indétermination

Remarque: Le calcul du Durbin-Watson suppose la normalité des perturbations et la présence d'une constante comme variable explicative.
ANNEXE N° 3

Nous recherchons le vecteur $\hat{a} = \begin{bmatrix} a_1 \\ \vdots \\ a_m \end{bmatrix}$ tel que $\hat{a}$ soit de norme minimale

et vérifie le système de n équations linéaires :

$$ \begin{bmatrix} I \end{bmatrix} \left\{ \begin{array}{l} y_1 = \sum_{i=1}^{m} a_i x_1^i \\ \vdots \\ y_n = \sum_{i=1}^{m} a_i x_n^i \end{array} \right. $$

d'où le problème de minimisation :

$$ \text{Min} \sum_{i=1}^{m} a_i^2 $$

sous les contraintes $[I]$

Par la méthode de Lagrange et puisque la fonction à minimiser est convexe et les contraintes linéaires les conditions du premier ordre sont nécessaires et suffisantes d'où :

$$ L(a_1, \dots, a_m, \lambda_1, \dots, \lambda_n) = \sum_{j=1}^{m} a_j^2 - \sum_{j=1}^{n} \lambda_j \left( y_j - \sum_{i=1}^{m} a_i x_j^i \right) $$

$$ \frac{\delta L}{\delta a_i} = 2 a_i + \sum_{j=1}^{n} \lambda_j x_j^i = 0 \quad \forall i=1,m $$

$$ -\frac{\delta L}{\delta \lambda_j} = y_j - \sum_{i=1}^{m} a_i x_j^i = 0 \quad \forall j=1,n $$

d'où matriciellement :

$$ \begin{array}{l@{\quad}c} \text{m} & \left\{ \begin{bmatrix} 2 & \cdots & 2 \\ \vdots & \ddots & \vdots \\ 2 & \cdots & 2 \\ \hline 1 & \cdots & 1 \\ \vdots & \ddots & \vdots \\ n & \cdots & n \end{bmatrix} \begin{bmatrix} x_1^1 & \cdots & x_n^1 \\ \vdots & \ddots & \vdots \\ x_1^m & \cdots & x_n^m \\ \hline 0 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & 0 \end{bmatrix} \begin{bmatrix} a_1 \\ \vdots \\ a_m \\ \hline \lambda_1 \\ \vdots \\ \lambda_n \end{bmatrix} \right\} = \begin{bmatrix} 0 \\ \vdots \\ 0 \\ \hline y_1 \\ \vdots \\ y_n \end{bmatrix} \\ \text{n} & \left\{ \begin{bmatrix} 1 & \cdots & 1 \\ \vdots & \ddots & \vdots \\ 1 & \cdots & 1 \\ \hline 1 & \cdots & 1 \\ \vdots & \ddots & \vdots \\ n & \cdots & n \end{bmatrix} \begin{bmatrix} a_1 \\ \vdots \\ a_m \\ \hline \lambda_1 \\ \vdots \\ \lambda_n \end{bmatrix} \right\} \end{array} $$

m

n
ce qui est équivalent à :

$$ \begin{bmatrix} 2 & I \\ X & 0 \end{bmatrix} \begin{bmatrix} a \\ \Lambda \end{bmatrix} = \begin{bmatrix} 0 \\ Y \end{bmatrix} $$

En utilisant la méthode d'inversion par blocs nous obtenons :

$$ \begin{bmatrix} a \\ \Lambda \end{bmatrix} = \begin{bmatrix} \frac{1}{2} & (I - X^{-1}(XX')^{-1}X) & X^{-1}(XX')^{-1} \\ & X^{-1}(XX')^{-1} & -2(XX')^{-1} \end{bmatrix} \begin{bmatrix} 0 \\ Y \end{bmatrix} $$

d'où : $\hat{a} = X^{-1}(XX')^{-1}Y$
BIBLIOGRAPHIE

[1] BROWN, DURBIN - "Methods for investigating whether a regression relationship is constant over time",
Congrès européen de statistique - Amsterdam 1968.

[2] BENZECRI, CAZES - "Méthodes statistiques de la taxinomie",
LSM Paris VI 1974.

[3] CHARLES C - " Regression typologique et reconnaissance des formes",
Thèse de 3ème cycle - Paris IX 1977.

[4] CHOW - "Test of equality between two sets of coefficients in two linear regressions",
Econometrica 28 (Juillet 1960) p. 561-605.

[5] DIDAY - "Nouvelles méthodes et nouveaux concepts en classification auto-
matique et reconnaissance des formes",
Thèse d'Etat 1972.

[6] DIDAY - "Classification et sélection de paramètres sous contraintes"
Rapport de Recherche n° 188 - IRIA-LABORIA 1976.

[7] FAIR, JAFFEE - "Methods of estimation for markets in desequilibrium",
Econometrica.

[8] FAIRLEY, HINICH - "A test for a shifting slope coefficient in a linear model",
JASA 65 (Sept 1970) 1320-29.

[9] GOLDFELD, QUANDT - "Non linear methods in econometrics",
Amsterdam - North Holland Publishing Co. (1972).

[10] Mc. GEE, CARLETON - "Piecewise regression",
JASA, 65, Sept 70, 1109-24.

[11] MARQUARDT - "Generalized inverses, ridge regression. Biased linear estimation and non linear estimation",
Technometrics, vol 12, N° 3, p. 591-611.

[12] OK Y. - "L'analyse factorielle typologique",
Thèse de 3ème cycle 1975.

[13] QUANDT - "The estimation of the parameters of a linear regression system obeying two separate regimes",
JASA 53 (Dec.58) 873-880.