Mathématice, intégration des Tice dans l'enseignement des mathématiques  
Sommaire > N°41 - septembre 2014 > Regards croisés sur l’algorithmique et la (...)

Regards croisés sur l’algorithmique et la programmation (6)
Les chaînes de Markov au bac S 2014
Moteur de recherche
Mis en ligne le 27 août 2014, par Alain Busser, Angelo Laplace, Guillaume Connan, Hubert Raymondaud, Pierre-Marc Mazat, Stéphan Manganelli

Cet article peut être librement diffusé et son contenu réutilisé pour une utilisation non commerciale (contacter l’auteur pour une utilisation commerciale) suivant la licence CC-by-nc-sa (http://creativecommons.org/licenses/by-nc-sa/3.0/fr/legalcode)

Cette sixième édition de la rubrique Regards croisés sur l’algorithmique et la programmation est consacrée à un exercice du bac 2014 : Le sujet de spécialité du bac S Pondichery.

Le sujet mettait en place une chaîne de Markov en dimension 3 :

Noter qu’alors que des chaînes de Markov sont évoquées dans des sujets du bac ES depuis une quinzaine d’années, le sujet devient une mode au bac S puisqu’il a été vu également

  • au bac S Liban (exercice spécialité)
  • au bac S Amérique du Nord (exercice 4 "obligatoire", avec un terme constant, il ne s’agit donc pas stricto sensu de chaîne de Markov)
  • Un exercice sur une chaîne de Markov avait déjà été donné à Pondichery en 2013
  • Le sujet de spé maths du bac Métropole/Réunion portait, comme celui de l’Amérique du Nord, sur un "système d’évolution avec terme constant". Les idées développées ici lui sont donc applicables. Voir par exemple la version CoffeeScript...
  • Plus récent, cet article de Philippe Gay évoque l’historique de la notion, ainsi que des exemples intéressants et une bibliographie indispensable.

Contribution d’Alain Busser

Alors que le calcul est en général effectué dans le CPU (processeur principal de l’ordinateur), le GPU (processeur graphique) est plus spécialisé dans le calcul matriciel. Or toute tablette sous ARM est dotée d’un GPU, et son gestionnaire openGL permet de faire du calcul matriciel. Dans le cas particulier des matrices d’ordre 4, Blender 3D possède lui aussi des fonctions de calcul matriciel...

Sur un navigateur gérant webGL (Chrome recommandé), le calcul matriciel peut être délégué au processeur graphique via webGL. Ce qui permet de modifier en live les éléments de la matrice [1] et le vecteur initial, et voir l’effet immédiat sur la figure (le nuage de points est suggéré par des petites sphères de taille et de couleur différentes pour mieux les distinguer).

le fichier avec rendu webGL
HTML - 10.4 ko
corrigé en webGL
ouvrir de préférence avec Chrome ou Safari (qualité du graphisme)

webGL sans three.js

webGL permet de faire du calcul parallèle dans le GPU à partir du JavaScript de navigateur internet. Le principe est le suivant :

  • les données sont transformées par JavaScript en float32Array() où chaque élément est stocké dans 4 octets ; par exemple un vecteur de dimension 3 devient un tableau de 12 octets consécutifs.
  • les données d’entrée sont transférées, via un buffer d’openGLES, dans le GPU ; par exemple, un transfert de 12 octets peut être considéré par le GPU comme un vec3, un objet typique d’openGL.
  • les données concernant la géométrie d’une scène 3D (par exemple, des points représentés par des vec3) sont traités par un "vertex shader" qui réalise la projection sur un écran 2D ;
  • puis le résultat est passé à un "fragment shader" qui modifie les couleurs par algorithmes (par exemple pour des textures).
  • les résultats de ces calculs (il y en a plusieurs menés en même temps sur plusieurs cœurs du GPU) sont ensuite envoyés à l’écran de l’ordinateur.

Il y a donc un problème technologique : On sait comment entrer des données dans un programme webGL (un vertex -> un fragment) mais on ne sait pas comment sortir des données du GPU vers le JavaScript [2]. Alors on va transformer ces données en couleur et lire cette couleur à l’écran. Le codage des couleurs étant de la forme vec4(r,v,b,a) où r est la quantité de rouge (entre 0 et 1), v la quantité de vert et b la quantité de bleu (a est la composante alpha qui concerne la transparence et qu’on laissera à 1), cette technique s’appliquera bien ici, parce que le vecteur V représentant une distribution de probabilité, a bien ses composantes entre 0 et 1.

Première étape : Aller sur cet interpréteur webGL en ligne :

On voit le code en langage webGL d’un "fragement shader" (ça tombe bien, c’est lui qui s’occupe des couleurs), et à côté, l’effet qu’il produit sur Suzanne. Explications sur le shader initial (ça servira pour la suite) :

  1. ligne 1, la précision est positionnée sur la plus haute valeur possible (en général 32 bits) pour les "float" (les réels) ;
  2. la ligne 2 définit une variable "time" de type "float" (un réel) ; c’est cette variable qui fait tourner la tête à Suzanne !
  3. la ligne 3 définit une autre variable appelée "resolution", de type vec2 soit un vecteur de dimension 2 ;
  4. la ligne 4 définit un vec3 appelé fPosition (il vient de l’autre shader) ;
  5. la ligne 5 définit une variable "fNormal" de type vec3 (le vecteur normal à la face courante) ; ces variables existent en plusieurs exemplaires puisqu’on fait du calcul parallèle...
  6. la ligne 6 est vide, pour marquer une séparation entre les déclarations de variables et le programme proprement dit.
  7. la ligne 7 dit que le programme s’appelle "main" et ne renvoie aucune valeur.
  8. les lignes 8 et 10 délimitent le programme en question.
  9. La ligne 9 est donc l’unique ligne du programme : Elle affecte à la couleur de Suzanne ("gl_FragColor") la valeur obtenue en concaténant fNormal (on se rappelle que c’est un vecteur de dimension 3) et le nombre 1 : Le tout est donc un vec4, un vecteur de l’espace de dimension 4 (le "1" final est la transparence, on déclare donc que Suzanne est opaque).

Un exemple de modification aisée est de remplacer dans l’unique ligne du programme, le vecteur fNormal par le vecteur fPosition : La couleur dépend maintenant de l’endroit qu’on regarde et plus du vecteur normal. En emplaçant Suzanne par un tore (dans le menu en haut à droite) on obtient alors un silmaril du plus bel effet :

On a maintenant tout ce qu’il faut pour rapidement calculer la distribution de probabilité invariante du problème de bac :

  • déclarer une variable V de type vec3 et initialisée avec vec3(0.5,0.3,0.2)
  • déclarer une matrice M de type mat3 et initialisée de façon analogue.
  • remplacer le vecteur donnant la couleur par V :

Ce vecteur est donné par une couleur maron plutôt foncée. Ensuite, on rajoute au début du programme une boucle dans laquelle on multiplie U par M avant de le représenter par une couleur :

Une seule multiplication par M donne déjà à Suzanne un teint tout gris : La convergence est rapide.

Puis il suffit de modifier le nombre de passages dans la boucle pour voir la distribution limite représentée par un gris foncé :

À titre de comparaison, avec 200 itérations, la couleur ne change plus :

Voici le programme à copier-coller :

  1. precision highp float;
  2. uniform float time;
  3. uniform vec2 resolution;
  4. varying vec3 fPosition;
  5. varying vec3 fNormal;
  6. vec3 U = vec3(0.5,0.3,0.2);
  7. mat3 M = mat3(0.5,0.5,0.1,0.4,0.3,0.2,0.1,0.2,0.7);
  8.  
  9. void main()
  10. {
  11. for(int n=0;n<=8;n++){
  12. U = M*U;
  13. }
  14. gl_FragColor = vec4(U, 1.0);
  15. }

Télécharger

Mais quelle est cette couleur ? Les copies d’écran ci-dessus ont été faites avec Gimp, qui a un outil "pipette" permettant de connaître la couleur de Suzanne ; il donne ceci :

C’est vraiment du gris, les composantes rouge, verte et bleue étant égales à 87 sur 256 soit environ 0,34 chacune : Résultat peu précis si on le compare au résultat exact :

  • x ≈ 0,37 (et non 0,34)
  • y≈ 0,30 (et non 0,34)
  • z ≈ 0,33 (plutôt bon)

Il existe une bibliothèque JavaScript permettant de faire des statistiques à l’aide de webGL].

L’outil MathsOntologie n’est pas le seul traité dans cet article, la manip est presque la même en html, en utilisant SVG ; les fichiers html obtenus sont aussi donnés en téléchargement :

l’article en pdf son source en odt l’algorithme en JavaScript La résolution du système en JavaScript
PDF - 510.2 ko
OpenDocument Text - 293.6 ko
HTML - 954 octets
HTML - 839 octets

Bonus : Des fichiers html, évoqués dans l’article ci-dessus, permettant de résoudre des systèmes linéaires sans passer trop de temps à entrer les coefficients :

Système 2×2 Système 3×3
HTML - 1.4 ko
HTML - 2.2 ko

L’outil alcoffeethmique a encore une fois été utilisé ici, mais il s’agit d’une version enrichie du calcul matriciel, téléchargeable ci-dessous :

Zip - 390.7 ko
l’article en pdf son source en odt
PDF - 378.7 ko
OpenDocument Text - 351.3 ko
  • La Ti 82 Stat fr

Matrices

La Ti 82 Stat Fr possède un mystérieux bouton

(vers le haut à gauche, juste à côté du bouton "math"). On va voir comment ce bouton permet de rapidement tester l’algorithme matriciel du sujet.

En dimension 3

En appuyant sur ce bouton, on découvre les noms possibles pour les matrices, qui vont de A à I entre crochets. Les lettres A, B et C étant réservées pour des variables vues plus tard dans l’énoncé, on va utiliser la lettre D pour désigner la matrice stochastique. On va la créer en allant dans le menu "édition" qui permet de modifier les matrices (et la création d’une matrice est une modification !), et, dans ce menu, on choisit D :

Par défaut, D est une matrice carré d’ordre 1 (c’est-à-dire un nombre...) :

On remplace alors les "1" par des "3" au clavier :

Alors il ne reste plus, toujours à l’aide du clavier, qu’à entrer les 9 coefficients de D :

Alors le menu des noms de matrices montre que D existe puisque ses dimensions sont affichées :

Pour calculer la puissance quatrième de D, on utilise la notation habituelle de la Ti :

Même les deux premières lignes montrent que lorsque l’exposant est suffisamment grand, les lignes de la puissance sont constantes, et représentent la distribution d’équilibre de la chaîne de Markov :

En effet, en multipliant une matrice dont les lignes sont constantes par un vecteur colonne dont la somme des éléments vaut 1, le produit est égal à l’une quelconque des colonnes de la matrice.

Algorithme du sujet

A est une matrice à 2 lignes et 2 colonnes, B et U sont des matrices à 2 lignes et 1 colonne (soit, des vecteurs colonne), comme on le voit dans l’énoncé. On va donc commencer par entrer A, en allant dans le menu d’édition des matrices, en entrant 2 et 2 pour ses dimensions et en entrant ses 4 coefficients :

Idem pour B :

Puis pour U, rebaptisé D puisque la lettre U n’existe pas dans l’espace de noms des matrices :

Pour programmer l’algorithme, il suffit, comme à l’accoutumée, d’aller dans le menu "programme" de la Ti, et d’y créer un programme :

Le nombre de lettres disponibles oblige à appeler quelque peu familièrement "Pondiche" le programme en question :

Le programme est alors concis, en utilisant les opérations de la calculatrice, et en cherchant les matrices dans le menu des noms de matrices :

Pour tester le programme, comme d’habitude on quitte l’éditeur de programme et on lance le programme :

La suite des valeurs successives prises par le vecteur U (ici, D) apparaît alors à l’écran :

Calcul de N

Dans le menu "math" du menu des matrices, se trouve l’instruction "identité" :

Elle permet d’entrer la matrice I de l’énoncé (en fournissant sa dimension comme argument), à laquelle il suffit ensuite de soustraire A pour avoir N :

Calcul de C

Pour voir comment a été obtenue N, voir l’onglet précédent. Il suffit de multiplier l’inverse de N par B pour avoir C :

On peut également trouver C en résolvant le système

  • 0,6x-0,4y=0,1
  • -0,2x+0,9y=0,2

Pour cela, on entre les 6 coefficients (en comptant le second membre) dans une matrice à 3 lignes et 2 colonnes (ici C) :

Pour résoudre le système, on applique une élimination de Gauss-Jordan à cette matrice. Celle-ci se trouve dans le menu mathématique des matrices :

La troisième colonne de C contient alors la solution du système :

Le logiciel Scilab étant spécialisé dans le calcul matriciel, devrait facilement résoudre ce genre d’exercice.

Voici une méthode possible.

Pour entrer M, pas besoin de mettre des crochets dans les crochets, mais ce sont des points-virgule qui délimitent les lignes de la matrice :

  1. M=[0.5,0.5,0.1;0.4,0.3,0.2;0.1,0.1,0.7]

Pour entrer V, on peut faire pareil, ou séparer ses éléments par des virgules et transposer le vecteur ligne ainsi obtenu pour avoir un vecteur colonne. La transposition en Scilab se note par un "prime" après le vecteur ligne :

  1. V=[0.5,0.3,0.2]'

Pour avoir directement le vecteur de probabilité au bout du quatrième mois, on peut passer par la puissance quatrième de M :

  1. M**4*V

ou alors, boucler en multipliant à chaque fois V par M pour avoir le nouveau V :

  1. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">for n=1:4, V=M*V, scilab.org/product/dic-mat-sci/M2SCI_doc.htm">end

Pour vérifier que les puissances de M tendent vers une limite, on peut afficher les valeurs propres de M (ou son "spectre") par

  1. spec(M)

ou faire calculer en plus la matrice des vecteurs propres par

  1. [vp,ep]=spec(M)

vp est une matrice diagonale, les valeurs propres étant sur la diagonale, et ep est la matrice de passage, dont les colonnes sont les vecteurs propres correspondant aux valeurs propres de vp.

L’algorithme de l’énoncé s’écrit alors

  1. A=[0.4,0.4;0.2,0.1]
  2. B=[0.1;0.2]
  3. U=[0.5;0.3]
  4. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">for n=1:5, U=A*U+B, scilab.org/product/dic-mat-sci/M2SCI_doc.htm">end

Télécharger

C’est assez concis ! On remarque que tout est matriciel dans ce script, puisque la boucle porte sur n qui parcourt le vecteur ligne 1:5 qui est formé des nombres entiers allant de 1 à 5...

Pour finir, on va résoudre le système NC=B avec Scilab. Avec N=I-A, on crée la matrice identité I avec

  1. I=eye(2,2)

Puis on effectue la soustraction pour avoir N, et on résout le système avec linsolve :

  1. I=eye(2,2)
  2. N=I-A
  3. [C,kerN]=linsolve(N,-B)

Télécharger

On obtient alors le noyau de N (kerN) qui ici est vide, et la solution C cherchée.

Pour en savoir plus sur la résolution par Scilab, voir la contribution d’Angelo Laplace ci-dessous.

Géométrie dynamique

La transformation itérée dans l’algorithme de l’énoncé est une affinité plane. Pour voir l’effet qu’elle produit, bouger les sommets du quadrilatère bleu ci-dessous, et observer l’effet obtenu sur son image (en marron) :

On peut également voir comment le cercle bleu est transformé en l’ellipse marron :

Outre le fait que la transformation ci-dessus n’est pas une isométrie, on constate son côté "contractante" : elle rapproche tout d’une droite ; ce qui laisse conjecturer la convergence de la suite de points :

Contribution de Guillaume Connan

Guillaume Connan propose un générateur de texte
aléatoire selon une distribution "markovienne" : on étudie les
successeurs de chaque couple de lettres successives apparaissant dans un
texte et on en fait une table de transition pour générer un texte.

Le dossier ci-dessous comprend le script Python pour

  • télécharger des livres au format "txt"
  • engendrer la table de transition (statistiques !)
  • Simuler la chaîne de Markov pour écrire un texte aléatoire "dans le style de..."
Le fichier Python, zippé
Zip - 1.1 ko

Le script Python en ligne

Pour consultation en ligne, voici l’intégralité du script Python qui produit des textes aléatoires :

  1. from codecs import open
  2. from urllib.request import urlretrieve
  3. from random import choice
  4. from itertools import chain
  5.  
  6.  
  7. class Markov(object):
  8.  
  9. def __init__(self,fichier):
  10. self.dicAssoc = {}
  11. self.fichier = fichier
  12. self.lettres = self.fichierVersLettres()
  13. self.baseDeTuplets()
  14.  
  15. def fichierVersLettres(self):
  16. """
  17. Renvoie les lettres du texte sous forme d'un itérateur
  18. """
  19. data = self.fichier.read() # le fichier entier de type string
  20. lettres = map(lambda lettre: ','.join(lettre),data) # on sépare chaque caractère du texte par une virgule
  21. return list(chain.from_iterable(lettres)) # on en fait un itérateur
  22.  
  23. def tuplets(self):
  24. """
  25. crée un générateur de triplets de lettres successives
  26.  
  27. """
  28. for i in range(len(self.lettres) - 2):
  29. yield (self.lettres[i], self.lettres[i+1], self.lettres[i+2])
  30.  
  31. def baseDeTuplets(self):
  32. """
  33. Crée un dictionnaire dont les clés sont des couples de lettres successives et la valeur la liste des successeurs observés.
  34.  
  35. """
  36. for l1,l2,l3 in self.tuplets():
  37. key = l1,l2
  38. if key in self.dicAssoc:
  39. self.dicAssoc[key].append(l3)
  40. else:
  41. self.dicAssoc[key] = [l3]
  42.  
  43. def genereMarkovText(self, nbLettres = 100):
  44. """
  45. génère un texte selon la distribution du dictionnaire d'association
  46. """
  47. (l1,l2) = choice(list(self.dicAssoc.keys())) # on choisit un couple de lettres au hasard dans le dico
  48. gen_lettres = "" # on initialise le texte généré
  49. for i in range(nbLettres):
  50. gen_lettres += l1 # on écrit l1 dans letexte
  51. l1,l2 = l2,choice(self.dicAssoc[(l1,l2)]) # on avance d'un cran
  52. print(gen_lettres)
  53.  
  54. # Notre-Dame de Paris
  55. urlretrieve('http://www.gutenberg.org/cache/epub/19657/pg19657.txt','ndp.txt')
  56. # Fahrenheit 451
  57. urlretrieve('http://sami.is.free.fr/Oeuvres/Bradbury%20Ray%20-%20Fahrenheit%20451-FR.txt','451.txt')
  58. # Gone with the wind
  59. urlretrieve('http://gutenberg.net.au/ebooks02/0200161.txt','autant.txt')
  60.  
  61. ndp = Markov(open('ndp.txt','r'))
  62. # Attention : celui-ci n'est pas en UTF-8 mais en latin-1
  63. fahr = Markov(open('451.txt','r','iso-8859-1'))
  64.  
  65. autant = Markov(open('autant.txt','r'))
  66.  
  67. # exemple : ndp.genereMarkovText(1000)

Télécharger

Par exemple avec le texte Fahrenheit 451, une production possible est celle-ci :

« Qui qu’un mait qu’occontre. Une res la vierpe.
« .........
« vous d’aper ses pas à per vie, de.
- Je pas rer de lactement vait pro des, des il éti, mous ré debont vid ond et un en nonts ! Oh. Je cons passimps, aut de quissez vour poncomme qui le dés suparquelquant sus de mons taiser, les détrevitre ça leurnivrengte ler le le ous verchoi. Cler eu deureil si trut une qu’illes sildrepte lesse touclaçongt-elé avouelles d’ennant cole. Au ! »

« Eh bous jouccifé landit seu main imin pâleillest et brés un d’ouls.. Il ricule ?
- Aus que nous l’on. Que pareille cycheffroi. « Voyager la s’ait et Commer lesses de proffou sur ses yeux bêtes l’étres ques quellibris gue mard, un a port. Et à tren conges la plusemps tres lever ! » d’étaisse larieur quatu ? Et gée.
- Mongez à cour aude vintie, ispar dont ge sa pare des côté ; commoigis per réva ser, n’avous permes saine Mon mace par s’air ne sallettelomment le que chos vis, lagerrente mait une dis puis de nourqueluier que d’houvrant à nouvenstre

Contribution d’Angelo Laplace

Angelo Laplace utilise Scilab pour faire du calcul matriciel, mais aussi Maple avec lequel il applique les chaînes de Markov à d’autres exemples.

Introduction

Les chaînes de Markov servent à modéliser les phénomènes qui se déroulent dans le temps de manière à ce que le passé et le futur à un instant donné soient conditionnellement indépendants par rapport au présent. On dit que c’est un processus "sans mémoire". Ainsi, pour faire une prédiction à l’instant n sur le comportement futur de la chaîne, il suffit de connaître la valeur $X_n$ et il est inutile de connaître les valeurs antérieures. Ceci est réalisé, en première approximation au moins, dans de nombreuses situations. On rappelle qu’une chaîne de Markov est homogène si la matrice de transition qui définit le passage de l’instant n à son successeur n+1 ne dépend pas de la valeur de n.

Dans le sujet de bac que nous étudions dans cette chronique, nous avons affaire à un processus markovien homogène où l’ensemble des états est fini (3 états, les marques de petits pots pour bébés). Je ne me propose pas de reprendre les calculs matriciels 2X2 suggérés dans l’énoncé de ce sujet, Alain Busser notamment l’ayant fait dans sa contribution. Je vais plutôt essayer de montrer comment les chaînes de Markov peuvent apporter des informations marketing précieuses en simulant notamment mois par mois l’évolution du marché de manière à en trouver les tendances à long terme. Pour les programmes, je me suis orienté vers deux logiciels qui gèrent les calculs matriciels 3X3 et qui m’étaient familiers : Scilab et Maple. Je me suis ensuite écarté du sujet de Bac et j’ai essayé d’en imaginer un prolongement qui pourra, je l’espère, donner des idées aux fans de programmation.

Quelques éléments théoriques

Le théorème de Perron-Frobenius s’applique aux matrices carrées A = $(a_{ij})$ d’ordre n qui vérifient les hypothèses suivantes :

  • tous les éléments $a_{ij}$ sont positifs ou nuls ;
  • il existe au moins une puissance à exposant entier de A tous les éléments sont strictement positifs.

Sous ces conditions, le théorème affirme que A possède une valeur propre réelle dominante $\lambda$ associée à un vecteur propre dominant V. C’est à dire que toutes les autres valeurs propres ont un module strictement inférieur à celui de $\lambda$ et de plus toutes les composantes de V sont strictement positives et de somme égale à 1.

Ce théorème permet de montrer qu’une chaîne de Markov à espace d’états fini converge en loi vers son unique distribution stationnaire (ou invariante). C’est alors $\lambda$ = 1 qui est la valeur propre dominante. Cette convergence est indépendante de la loi initiale choisie pour un processus markovien donné. On peut donc rien qu’en itérant la matrice atteindre une très bonne approximation de la distribution stationnaire, ceci simulant d’ailleurs l’évolution du système markovien au fil du temps.

L’exercice posé au bac : recherche de la mesure stationnaire

Dans l’exercice de bac S proposé en 2014 à Pondichéry, la matrice de transition entre les différentes marques obéit aux hypothèses du théorème de Perron-Frobenius puisque tous ses coefficients sont strictement positifs :

$A=\left(\begin{array}{ccc} 0.5 & 0.5 & 0.1 \\ 0.4 & 0.3 & 0.2 \\ 0.1 & 0.2 & 0.7 \end{array}\right)$

Dans un premier temps, il est facile de chercher la distribution stationnaire : il suffit pour cela de chercher un élément du noyau de A-Id et de le normaliser de manière à ce que la somme de ses composantes soit égale à 1.

Pour réaliser cette tache, j’ai utilisé le logiciel de calcul numérique gratuit scilab : celui-ci peut être téléchargé à cette adresse. Lors d’une utilisation pédagogique avec les élèves de terminale en enseignement de spécialité, l’installation du "Module Lycée" de Claude Gomez et le téléchargement d’un guide pratique peut s’avérer utile. La particularité de logiciel est d’être particulièrement bien adapté au calcul matriciel puisque les nombres en tant que tels sont pour scilab des matrices 1X1 et qu’ainsi tous les calculs sont définis selon un procédé matriciel. Cela en fait un outil pertinent pour l’étude des chaînes de Markov, a fortiori dans les établissements scolaires (lycée) et universitaires où la gratuité est un atout non négligeable. Notons tout de suite que la commande "afficher" que nous utiliserons plus bas est propre au module lycée et qu’il faut impérativement l’installer pour pouvoir l’utiliser.

Une fois tout installé, on clique sur le menu Applications puis dans le menu déroulant sur SciNotes, une fenêtre d’édition s’ouvre pour saisir la suite d’opérations algorithmique à réaliser. La prise en mains du logiciel est assez simple lorsque l’on connaît la syntaxe propre à l’édition des matrices.

Pour créer la matrice A et le vecteur indiquant la distribution initiale des clients, il suffit de saisir :

JPEG - 30 ko

Notons que le "double slash" est utilisé pour écrire un commentaire qui n’est pas pris en compte par le logiciel lors de l’exécution.

La recherche du noyau dans scilab s’effectue assez naturellement en utilisant la fonction kernel.

Le programme scilab suivant réalise donc la recherche de la distribution stationnaire :

  1. // saisie de la matrice de transition
  2. A=[0.5,0.5,0.1;0.4,0.3,0.2;0.1,0.2,0.7];
  3. // saisie de la matrice identité
  4. I=[1,0,0;0,1,0;0,0,1];
  5. W=A-I;
  6. // calcul du noyau de A-I
  7. X=kernel(W);
  8. // procédure de normalisation du vecteur
  9. norme=X(1)+X(2)+X(3);
  10. stationnaire=X/norme;
  11. // affichage des composantes de la distribution stationnaire
  12. afficher(stationnaire);

Télécharger

Voici ce que la console de scilab affiche comme résultat à la suite de l’exécution de ce programme :

JPEG - 18.6 ko

Si vous êtes un peu déçus par les nombreuses décimales affichées et vous auriez préféré une valeur exacte. Alors, il faut se tourner vers un logiciel de calcul formel. J’y reviendrai un peu plus bas. En réalité, les composantes de la probabilité stationnaire sont exactement $17 \over 46$, $7 \over 23$ et $15 \over 46$. Le logiciel scilab en donne donc une très bonne approximation.

En tout cas, après quelques secondes de travail, nous pouvons déjà prédire que le marché va se stabiliser ainsi : au bout d’un certain temps, environ 37% des clients vont se fournir chez la marque X, 30% des gens s’orienteront vers la marque Y et 33% des parents prendront des petits pots de la marque Z. De mois en mois, même si certains clients décident ponctuellement de changer de marque, le marché restera stable autour de ces proportions.

Convergence en loi vers la mesure invariante

La seconde question qui peut nous occuper est de déterminer à partir de combien de mois cette position d’équilibre est atteinte ou plutôt pratiquement atteinte. Il s’agit donc d’étudier la vitesse de convergence vers la distribution stationnaire, convergence qui ,doit-on le rappeler, est assurée par la théorie.

Étant donné que la distribution initiale est connue, pour étudier la convergence, je me fixe une précision et je regarde le nombre d’itérations nécessaires pour atteindre cette précision prédéterminée. J’utilise donc une boucle "while" pour gérer ce processus itératif qui doit s’arrêter naturellement dès que la précision voulue est atteinte.
Une fois la sortie de la boucle effective, il convient de faire afficher le nombre de mois nécessaires pour réaliser l’objectif et les valeurs obtenues de la répartition probable des clients à l’issue du processus.

Voici l’algorithme implémenté en langage scilab :

  1. // saisie de la distribution initiale en janvier
  2. U=[0.5;0.3;0.2];
  3. //saisie de la matrice de transition 3X3
  4. A=[0.5,0.5,0.1;0.4,0.3,0.2;0.1,0.2,0.7];
  5. //saisie des composantes de la distribution stationnaire
  6. f=17/46; g=7/23; h=15/46;
  7. //calcul de la distribution à la fin du premier mois
  8. V=A*U;
  9. //initialisation du compteur d'itérations(nombre de mois écoulés)
  10. i=1;
  11. //boucle for sur la précision attendue des calculs
  12. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">for k=1:20
  13. //boucle while pour la condition d'arrêt
  14. //tant que nous sommes éloignés de la distribution stationnaire, on applique A
  15. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">while abs(V(1)-f)>10^(-k)
  16. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">while abs(V(2)-g)>10^(-k)
  17. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">while abs(V(3)-h)>10^-k i=i+1; V=A^i*U;
  18. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">end
  19. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">end
  20. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">end
  21. //Affichage du nombre de mois écoulés, de la précision et de la distribution obtenue
  22. afficher("En "+string(i)+" itérations, le système est à 10^"+string(-k)+" près de la ..
  23. distribution stationnaire.")
  24. afficher(V(1)); afficher(V(2)); afficher(V(3));
  25. scilab.org/product/dic-mat-sci/M2SCI_doc.htm">end

Télécharger

Nous pouvons tirer de cette étude l’enseignement suivant : au bout de 9 mois (9 itérations pour une précision de 0,001), soit au mois d’octobre, le marché est déjà très stabilisé entre les différentes marques. Le vecteur issu de cette simulation du marché jusqu’en octobre a pour composantes :

0.3701577252,

0.3045998733,

0.3252424015.

Cela signifie que 37% des acheteurs environ s’orientent vers la marque X, 30,4% achètent la marque Y et 32,6% préfèrent la marque Z. Ensuite, même si des consommateurs prennent ponctuellement et individuellement la décision de changer de marque, le marché est globalement stable à 0,1% près, ceci tant que la matrice de transition qui répartit les échanges de clients reste inchangée. Au bout de 13 mois, le marché est même stable à 0,01% près. La copie d’écran suivante reprend d’ailleurs les résultats obtenus avec scilab :

JPEG - 89.3 ko

Par contre, dans scilab, tous les calculs sont numériques et sont effectués en double précision. Cela signifie que la précision maximale est d’environ 16 chiffres décimaux. Un utilisateur novice pourrait penser que la saisie de format(20) permettra d’obtenir une meilleure précision. En fait, il n’en est rien comme le prouve cet exemple :

JPEG - 11.1 ko

En réalité, les 18 premières décimales de $\pi$ sont :

3,14159265358979323

On voit donc que scilab présente 2 décimales "terminales" erronées. D’ailleurs, en laissant la commande format(20), dès la fin de la première itération dans mon programme, le résultat est plutôt surprenant puisque nous attendions notamment 0.42 exactement pour la première composante (voir énoncé).

JPEG - 24.2 ko

L’erreur est certes minime mais l’erreur est réelle.

On peut d’ailleurs voir que scilab ne peut aller plus loin qu’une précision à $10^{-15}$. Il s’arrête sans message d’erreur d’ailleurs.

En conclusion, nous voyons que scilab est un logiciel qui permet d’obtenir rapidement une très bonne valeur approchée à $10^{-3}$ près de la distribution stationnaire. Mais pour aller au-delà, il faudra immanquablement se tourner vers un logiciel de calcul formel.

Dans la famille des logiciels de calcul formel, il existe des logiciels gratuits comme Maxima ou Derive...et des logiciels propriétaires comme Maple et Mathematica. Pour la suite de ce propos, je me suis tourné vers Maple car c’est un logiciel qui a fait fureur en classe préparatoire et fut utilisé également dans l’épreuve de modélisation de l’agrégation externe.

La première vertu que l’on peut trouver à Maple est de pouvoir donner les valeurs exactes des composantes de la distribution stationnaire. Le programme suivant se charge de donner un vecteur propre de la matrice A-Id pour la valeur propre 1.

>restart ;
>with(linalg):
>A:= matrix(3,3,[5/10,5/10,1/10,4/10,3/10,2/10,1/10,2/10,7/10]);
>J:=array(identity,1..3,1..3);
>K:=kernel(A-J);

Le programme sort un vecteur dont la somme des composantes n’est pas 1, mais $23 \over 7$, il convient donc de la normaliser pour obtenir la distribution dont j’ai déjà parlé.

>P:=vector([17/14,1,15/14]):
>P[1]+P[2]+P[3];

sortie : $23 \over 7$

>evalm(7/23*P);

Le package linalg est indispensable pour gérer les matrices et utiliser la fonction kernel dans Maple. On peut observer quelques différences dans la manière de saisir une matrice par rapport à scilab.

Avec Maple, on peut également réaliser l’’itération du processus markovien de manière à retrouver les résultats relatifs à la convergence de la chaîne :

>U[0]:=vector([5/10,3/10,2/10]);
>for i from 1 to 4 do U[i]:=A&*U[i-1] ; od :

Cette boucle "for" limitée à 4 itérations permet de vérifier les valeurs numériques demandées dans la dernière question du sujet de Bac après avoir complété avec cette commande :

>V:=evalm(U[4]);  for k from 1 to 3 do evalf(V[k]); od;

En programmant cette boucle jusqu’à i=9, on retrouve les résultats obtenus par scilab précédemment.

Et si on change la distribution initiale ?

Dans le paragraphe précédent , nous avons vu que la convergence vers la loi stationnaire est extrêmement rapide. On pourrait croire que le choix de la distribution initiale y est pour beaucoup. Pour réfuter cette hypothèse, je propose d’examiner l’exemple suivant.

Supposons que le nombre de clients qui achètent la marque X (celle qui est la plus populaire) est tombé à 0 (par exemple à cause d’une campagne publicitaire désastreuse mettant en scène une célébrité qui n’est plus du tout populaire ou à la suite d’un scandale alimentaire.) Prenons donc pour distribution initiale le vecteur U(0 ; 0.6 ; 0.4) et supposons que nous sommes au moment où les consommateurs commencent à oublier les problèmes rencontrés par la marque X et que la matrice de transition A gère de nouveau les échanges entre les différentes marques. On se souvient, que sous ces nouvelles conditions, la distribution stationnaire reste inchangée. Les résultats fournis par scilab sont sans appel :

JPEG - 94.3 ko

Même si après une itération, la première composante associée à X est à 0.36 au lieu de 0.42, nous constatons la convergence attendue vers la distribution stationnaire (c’est indépendant de la distribution initiale, rappelons-le) et même un retour proche à $10^{-3}$ de l’équilibre en 9 itérations également. Tant que le système est régi par la matrice A, il est convenu qu’en 9 mois, le marché est "stabilisé". Le seul moyen pour les marques de faire changer cet équilibre du marché est d’agir pour que les coefficients de la matrice leur soit plus favorable (soit faire en sorte que moins de clients les quittent, soit faire en sorte que plus de clients les rejoignent.)

Une quatrième marque ?

a) Garder les clients conquis

Élargissons un peu le champ d’application. Oublions un instant les petits pots pour bébé et leurs trois marques X,Y et Z. Les consommateurs grandissent vite !
Les voici maintenant dans un environnement où 3 opérateurs de téléphonie mobile X,Y et Z se partagent le marché quand un quatrième larron W vient entrer dans la danse (toute ressemblance avec des faits réels ne serait que pure et fortuite coïncidence). Nous supposons que la distribution initiale est (0,5 ; 0,3 ; 0,2 ; 0), ce qui est la même distribution que pour les petits pots de bébé mais avec une quatrième composante nulle. L’opérateur X est en tête et possède 1 client sur 2, l’opérateur Y a une courte avance sur l’opérateur Z tandis que W n’a aucun client puisqu’il vient d’être créé. Mais W a une stratégie agressive et pratique des tarifs alléchants (pure coïncidence fortuite, on vous dit ! ) Par contre, les nouveaux clients de W sont parfois déçus par le service et retournent vers la concurrence. Nous supposons en outre que tous les pourcentages concernant les transitions de clientèle entre les opérateurs X, Y et Z baissent de 5% par rapport à l’exercice des petits pots pour bébé. Ainsi 15 % des clients de chacun des opérateurs X,Y et Z cèdent au chant des sirènes et décident de les quitter pour rejoindre W.

Un client de l’opérateur X a le mois suivant :

  • 45% de chance de rester fidèle à X ;
  • 35% de chance de changer pour s’engager envers Y ;
  • 5% de chance de changer pour s’engager envers Z ;
  • 15% de chance de rejoindre le nouvel opérateur W.

Un client de l’opérateur Y a le mois suivant :

  • 25% de chance de rester fidèle à Y ;
  • 45% de chance de changer pour s’engager envers X ;
  • 15% de chance de changer pour s’engager envers Z ;
  • 15% de chance de rejoindre le nouvel opérateur W.

Un client de l’opérateur Z a le mois suivant :

  • 65% de chance de rester fidèle à Z ;
  • 5% de chance de changer pour s’engager envers X ;
  • 15% de chance de changer pour s’engager envers Y ;
  • 15% de chance de rejoindre le nouvel opérateur W.

Quant aux clients de l’opérateur W, nous supposons qu’ils quittent W, insatisfaits du service en égale proportion vers chacun des concurrents.

Les dirigeants de l’entreprise W doivent faire une étude précise du marché pour savoir si limiter les pertes de clients en proposant un service de qualité est une stratégie intéressante pour gagner le leadership. Je me propose d’aider les dirigeants de W en rédigeant un programme avec Maple.

Le paramètre $a$ sert à mesurer le niveau de satisfaction des clients de W. Par commodité, je choisis de dire que $3a+1$% des clients de W se réengagent le moins suivant. En conséquence $33-a$% des clients de W repartent vers chacun des 3 concurrents.

On réalise deux procédures Maple :

  • La première "Vecteurpropre" se charge d’écrire la matrice de transition A, dépendante du paramètre $a$, puis elle cherche le noyau de A - Id. Elle normalise le vecteur propre de manière à en faire une distribution de probabilité.
  • La seconde procédure "Meilleursparametres" est chargée de renvoyer les paramètres $a$ pour lesquels le l’opérateur nouvellement créé devient le premier en nombre d’utilisateurs à la distribution d’équilibre.

Voici le programme complet (la distribution initiale n’intervient pas) :

>restart;
>with(linalg):

>J:=array(identity,1..4,1..4):

>Vecteurpropre:=proc(a)
local A,V,W,X,N,k;
A:= matrix(4,4,[45/100,45/100,5/100,(33-a)/100,35/100,25/100,15/100,(33-a)/100,5/100,15/100,65/100,(33-a)/100,3/20,3/20,3/20,(1+3*a)/100]):
V:=kernel(A-J): W:=op(1,V):
           for k from 1 to 4 do evalf(W[k]) od :
N:=sum(evalf(W[o]),o=1..4):  X:=evalm(1/N*W);
end;

>Meilleursparametres:=proc()
local a;
for a from 1 to 33 do Vecteurpropre(a); if  Vecteurpropre(a)[4]> Vecteurpropre(a)[1] and Vecteurpropre(a)[4]> Vecteurpropre(a)[2] and Vecteurpropre(a)[4]> Vecteurpropre(a)[3] then print(a); fi; od;
end;

>Meilleursparametres();

Maple donne en sortie de la procédure "Meilleursparametres" tous les entiers compris entre 20 et 33. En analysant les résultats, on s’aperçoit que l’opérateur qui part avec 0% du marché devient l’opérateur prépondérant à terme lorsque son taux de réengagement atteint 61%. En dessous, il restera dominé par au moins un concurrent.

On peut même obtenir la distribution stationnaire dans le cas où le taux de réengagement est de 61% et estimer l’écart par rapport aux autres opérateurs en saisissant :
>B :=Vecteurpropre(20) ;
>for i from 1 to 4 do evalf(B[i]) od ;

L’opérateur W devance X d’une courte tête 27,8% des clients contre 26,7% à l’équilibre du marché.

Voici même une procédure Maple supplémentaire "Leadership" qui permet de déterminer au bout de combien de mois, le quatrième opérateur devient le leader du marché :

>Leadership:=proc()
local A,U,i,V;
U[0]:=vector([5/10,3/10,2/10,0]);
A:= matrix(4,4,[45/100,45/100,5/100,13/100,35/100,25/100,15/100,13/100,5/100,15/100,65/100,13/100,3/20,3/20,3/20,61/100]);
> for i from 1 to 100 do U[i]:=A&*U[i-1]; V:=evalm(U[i]); if V[4]> V[1] then RETURN(i); fi; od;
> end;

Celle-ci retourne la première valeur pour laquelle le nombre de clients du quatrième opérateur dépasse celui du premier, à savoir 6. Il faut donc seulement 6 mois pour être le numéro 1 à cet opérateur sous ces conditions. On peut même s’apercevoir en calculant U[6] que ce mois-ci, cet opérateur possède déjà 27,5% des clients et est tout proche de son maximum.

Les dirigeants de W ont donc tout intérêt à proposer un service de qualité.

b) Casser les prix quitte à faire des insatisfaits ?

On peut quand même se demander si W a intérêt à casser les prix quitte à ne pas satisfaire ses nouveaux clients : on suppose ici que les clients se W se partagent le mois suivant en :

  • 25% qui se réengagent vers W ;
  • 25% qui quittent W pour X ;
  • 25% qui quittent W pour Y ;
  • 25% qui quittent W pour Z.

Cette fois-ci, nous baissons tous les pourcentages de transition entre opérateurs de a% par rapport à l’exercice des petits pots pour bébés. Ainsi l’opérateur W récupère automatiquement 3a% des clients de chaque opérateur. Le paramètre a sera limité à 10 car sinon certains éléments de la matrice de transition A deviendraient négatifs.

La matrice A s’écrit donc en multipliant celle-ci par le scalaire $1 \over 100$ :

$$\left(\begin{array}{cccc} 50-a & 50-a & 10-a & 25 \\ 40-a & 30-a & 20-a & 25 \\ 10-a & 20-a & 70-a & 25 \\ 3a & 3a & 3a &25 \end{array}\right)$$

La procédure "Vecteurpropre" devient :

Vecteurpropre:=proc(a)
> local A,V,W,X,N,k;
> A:= matrix(4,4,[(50-a)/100,(50-a)/100,(10-a)/100,25/100,(40-a)/100,(30-a)/100,(20-a)/100,25/100,(10-a)/100,(20-a)/100,(70-a)/100,25/100,3*a/100,3*a/100,3*a/100,25/100]): V:=kernel(A-J):
W:=op(1,V):
           for k from 1 to 4 do evalf(W[k]) od :
N:=sum(evalf(W[o]),o=1..4):   X:=evalm(1/N*W);
end;

La procédure "Meilleursparametres" donne en sortie seulement le nombre 10 : ainsi l’opérateur de téléphonie W doit impérativement "voler" 30% des clients de chaque concurrent pour devenir le leader du marché avec 28,6% contre 26,4% à X.

Une stratégie agressive de baisse des prix s’avère ainsi également efficace même dans le cas où 75% des clients récupérés s’avèrent finalement insatisfaits.

Conclusion

Les chaînes de Markov que nous avons croisées ici servent à modéliser en première approximation des situations que la vie courante génère (diffusion des gaz, files d’attente, modèles proies-prédateurs...) Dans cet exercice de bac de 2014 à Pondichéry ainsi que dans le prolongement que nous avons examiné, celles-ci sont étudiées à des fins de marketing. Bien entendu, la vie réelle n’est jamais aussi idéale et l’homogénéité de la chaîne qui maintient stable la matrice A de mois en mois est une vue de l’esprit que le mathématicien s’accorde à accepter en première approximation.

Je me souviens que lorsque j’étais en passe de tenter mon bac, il n’y avait aucune matrice, aucune chaîne de Markov au menu. Leur apparition récente dans les programmes de terminale doit sûrement quelque chose à l’économie de marché dans laquelle nous vivons. Quelques années plus tard, lorsque je rédigeais mon mémoire de maîtrise sur l’étude du Monopoly avec ces mêmes processus markoviens, j’étais bien loin de penser que les bacheliers auraient 15 ans plus tard à se pencher sur des phénomènes stochastiques. Je m’en félicite en tout cas.

Je parlais à l’instant de Monopoly.... on peut déterminer la probabilité d’occupation de chacune des cases du plateau de jeu. Essayez, vous verrez bien, c’est assez drôle : commencez par déterminer les probabilités de passer d’une case à une autre avec les dés, puis d’une case à une autre avec des cartes et vous aurez la matrice de transition, 40X40 quand même. Elle remplit les conditions du théorème de Perron-Frobenius et donc vous obtiendrez un vecteur propre qu’il faudra normaliser en utilisant un logiciel adapté. Après, vous verrez, vous aimerez les quartiers rouges et oranges ! Sacrées chaînes de Markov !

Contribution de Pierre-Marc Mazat

L’article au format pdf
PDF - 341.4 ko

Le document pdf ci-dessus prolonge et complète l’intervention de son auteur dans cet article : Celui-ci portait en effet déjà sur une chaîne de Markov...

Contribution de Brigitte Chaput, Stephan Manganelli et Hubert Raymondaud

Le problème des "chaînes de 6" (6 "pile" de suite à pile ou face) peut s’expliquer par une chaîne de Markov (dessinée dans le document de Stephan Manganelli ci-dessous). Il est analysé conjointement par LARP, R et un tableur :

La contribution de Stephan Manganelli Le fichier Excel (Brigitte Chaput) La contribution d’Hubert Raymondaud
PDF - 756.9 ko
Excel - 119.5 ko
PDF - 149.7 ko

notes

[1seul le cas où la matrice est d’ordre 3 est traité ; en effet three.js ne possède pas de matrices d’ordre 2, le processeur graphique étant spécialisé dans le calcul quaternionnique...

[2ce problème n’existe pas avec webCL qui est en cours de développement

Documents associés à l'article
     |   (HTML - 954 octets)
     |   (HTML - 839 octets)
     |   (HTML - 954 octets)
Réagir à cet article
Vous souhaitez compléter cet article pour un numéro futur, réagir à son contenu, demander des précisions à l'auteur ou au comité de rédaction...
À lire aussi ici
MathémaTICE est un projet
en collaboration avec
Suivre la vie du site Flux RSS 2.0  |  Espace de rédaction  |  Nous contacter  |  Site réalisé avec: SPIP  |  N° ISSN 2109-9197