Category Archives: Maths

Discussions géométriques du soir

J’avais raconté l’autre jour l’aventure de mon fils qui souhaitait faire des triangles avec des bouts de bois (ici). Il a récidivé ce soir. Sauf qu’il n’avait qu’un grand baton, et que j’avais fermé le portillon qui l’empêchait d’aller dans la forêt en chercher (la nuance entre romancer et devenir mythomane doit être mince). Je lui ai suggéré de le casser en trois portions, puis de constituer un triangle avec les trois morceaux restant. Manque de bol, l’expérience rate (comme l’autre fois). Et me souvenant de mes pénibles calculs de la dernière fois, je me suis rappelé qu’il avait une chance sur deux de pouvoir effectivement construire un triangle… sauf que la manière d’obtenir les trois bouts de bois me semble différente… donc peut être que ce 1/2 n’est plus valide ?
En fait, je crois qu’il faut définir clairement comment on va couper les bouts de bois.

  • méthode 1: le coupage simultané

Je suppose que je peux couper d’un coup le bout de bois en deux. Formellement, je tire deux points sur l’intevalle correspondant à mon bout de bois, les tirages étant uniformes, et indépendants. On arrive alors à trois morceaux, alatex2png.2.php.pngalatex2png-1.2.php.png et alatex2png-2.2.php.pngen sachant que 

alatex2png-3.2.php.png

la longueur du bout de bois dont je diposais initialement. On supposera que aalatex2png.2.png=1 pour simplifier (ce qui ne changera rien au résultat final, on s’en doute). L’ensemble des points que je peux obtenir peuvent alors se visualiser comme des points du simplexe, 

alatex2png-4.2.php.png

autrement dit, on peut les représenter dans le triangle équilatéral ci-dessous.
Or on ne peut constituer un triangle que si la plus grande des longeurs n’excède pas la somme des deux autres. Supposons que alatex2png-2.2.php.png soit la plus grande longueur (et qu’elle soit visualisée graphiquement par rapport à la distance au côté en bas), le point rouge ci-dessous (correpondant à un triplet de longueurs possibles) ne convient pas. En fait, dès lors que l’on entre dans l’aire rosée, on ne peut plus constituer un triangle…

Moralité, à moins de tomber dans le triangle jaune, on ne peut pas constituer de triangle…. Ce qui veut dire qu’on n’a plus qu’une chance sur quatre !

  • méthode 2: le coupage séquentiel

Une autre méthodde consister à couper d’abord en deux (en tirant un point au hasard sur mon intervalle), puis de recouper en deux le plus grand morceau… Damned, mon fils a de bonne idée, mais ça sent les calculs à n’en plus finir ça…Ou alors va falloir continuer à essayer de faire des dessins, car mine de rien, c’est plus simple que les calculs de la dernière fois. Bon, la longueur de mon bout de bois est la hauteur de mon triangle. Je commence par le couper en deux, la distance au côté en bas devant être la plus petite des deux.

Une fois fixé la taille de ce petit bout, il me reste à couper le plus grand morceaux en deux. Pour cela, je vais me déplacer le long du segment rouge ci-dessous,
La position sur ce petit segment me donnant ensuite les deux autres longueurs des morceaux. Or rappelons que pour que nos bouts de bois forments un triangle, il faut forcément appartenir au triangle jaune,
Pour résumer, nous avions exclus la partie supérieure du triangle, et conditionnellement au fait que nous sommes dans le lozange inférieur, nous souhaitons être dans le triangle jaune. Arg, cette fois-ci la probabilité passe à une chance sur trois !

  • couper d’un coup comme découpage séquentiel aléatoire ?

Mais attends deux minutes mon bonhomme (oui, je parle à mon fils). Si j’avais choisi de découper le plus petit morceau, on sait qu’il est impossible de constituer un triangle car par définition alatex2png-5.2.php.png Supposons que l’on casse le baton en deux, et que l’on tire ensuite au hasard (avec probabilité 1/2 et 1/2) quelle partie sera alors coupée en deux, on doit se ramener au cas précédant non ? La probabilité de constituer un triangle doit alors être

alatex2png-6.2.php.png

Damned, mais ça ne fait pas un quart ça ! Bon, je crois que je vais prendre une aspirine et aller demander à mon fils d’aller ranger ses bâtons !

Zénon, sommes infinies et loi géométrique

Ayant parlé plusieurs fois de Luc Chatel sur mon blog je n’ai pas pu m’enpêcher de mettre l’image qui est à la une du Canard de cette semaine.
Pour ceux qui se souviennent de leurs cours de philo au lycée, c’est le paradoxe de Zénon (ou presque). Pour cela il faut fermer les yeux et imaginer une flèche en vol (en fait non, je viens de faire l’expérience, si on ferme les yeux, on ne peut pas lire la suite, ça ne marche qu’en cours cette histoire). La flèche est lancée à 10 mètre de sa cible. Il faut un certain temps, non nul, à cette pierre pour parcourir la moitié de la distance qui la sépare de la cible. Ensuite, il lui reste encore 5 mètres à parcourir, dont elle accomplit d’abord la moitié, 2,5 mètres, ce qui lui prend un certain temps. Puis elle parcourt la moitié, puis la moitié, et cela – à chaque fois – avec un temps non nul. Zénon (Ζήνων, parfois appelé aussi Zénon d’Elée)  en conclut que la pierre ne pourra frapper l’arbre qu’au bout d’un temps infini, c’est-à-dire jamais (ou pour la version illustrée, je peux renvoyer à un extrait de Kid Paddle, ici). En fait, Zénon introduit  la notion d’infini divisibilité (on dira qu’un élément est indéfiniment divisible, c’est d’ailleurs surprenant ce lien en français entre l’infini et l’indéfini), qui philosophiquement, est un complexe qui a mis du temps à s’imposer: une somme infinie d’éléments non nuls peut être finie (je renvois ici par exemple).

Formellement, l’histoire de la flèche revient à se demander si on peut calculer quelquechose qui ressemble à

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/zenon-01.png

soit, de manière plus formelle,

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/zenon-02.png

on retrouve la somme d’une série géométrique (au premier terme près), dont on sait que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/zenon-03.png

Un application intéressante est que la loi géométrique, vérifiant

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/zenon-06.png

est effectivement une loi de probabilité,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/zenon-05.png

(les termes sont positifs, et la somme vaut 1). Il existe une démonstration élégante de ce dernière résultat, qui tient en un dessin (à condition d’accepter ce principe un peu fractal d’infini divisibilité):

Je reviendrais un jour sur les lois de probas infiniment divisibles, car elles sont très très intéressantes de par certaines propriétés….

Fibonacci: des lapins à l’actualisation

Depuis qu’on peut facilement commander des livres d’occasion outre-atlantique, j’avoue en abuser un peu… La dernière commande que j’ai faite est le “origins of value“. Cet ouvrage, qui est une compilation d’articles d’histoire de la finance, contient un chapitre passionnant sur Leonardo Pisano alias Fibonacci. Pour remettre dans le contexte, Fibonacci vivait il y a fort fort longtemps, autour de 1200 (300 ans avant la renaissance et Léonard de Vinci). Je peux renvoyer ici par exemple pour une analyse de quelques uns de ses ouvrages.
Pour moi, Fibonnaci, c’était l’histoire de la reproduction des lapins, c’est à dire, l’étude
des suites récurentes pour ceux qui se souviennent de leurs cours de lycée, avec des suites de la forme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fibo-001.png

Le problème que l’on cherche à résoudre – posé par Fibonacci – est “partant d’un couple, combien de couples de lapins obtiendrons-nous après un nombre donné de mois sachant que chaque couple produit chaque mois un nouveau couple, lequel ne devient productif qu’après deux mois ?”. Pour ceux qui se souviennent de la résolution, l’équation charactéristique de cette relation de récurence est un polynôme de degré 2,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fibo-002.png

dont les solutions sont

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fibo-003.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fibo-004.png

Le premier nombre est souvent connu comme étant le nombre d’or (ici). Bref, Fibonacci a été recyclé abondamment par l’utilisation ésotérique de cette racine de polynôme de degré 2. Aussi, lorsque j’ai vu évoqué le nom de Fibonacci dans cet ouvrage de finance, j’ai eu peu peur, d’autant plus que certaines personnes faisant de l’analyse chartiste utilisent ce nombre d’or (comme l’évoque rapidement certains site, comme http://www.forexfibonacci.com/ ou http://www.fibonaccitrader.com/). Aussi, parler de Fibonacci en finance, c’est un peu comme parler de Nostradamus dans les modèles de prévisions.

Que nenni mon ami ! Fibonacci est l’auteur d’un  Liber Abaci, correspondant à un livre des calculs. En fait, Leonardo s’est énormément inspiré des mathématiques arabes, alors que la norme était plutôt l’arithmétique romaine. Pour ceux qui ont essayé un jour dans leur vie de lire des nombres latins (même dans Astérix), la logique est assez déroutante ! Par exemple 888 s’écrit DCCCLXXXVIII… Essayez avec cette écriture de poser une addition ! Les nombres arabes – eux – utilisaient le système décimale, ce qui a fait gagné énormément de temps de calcul (on peut aussi remonter à l’arithmétique indienne, mais on sort largement du cadre de mon billet). Fibonacci évoque d’ailleurs constamment l’utilisation de la règle de trois (sans l’appeler ainsi, j’étais revenu ici sur l’origine de cette règle).
Dans son livre d’arithmétique des nombres arabes, Fibonacci montre une application intéressante: le calculs d’actualisation et des marges commerciales. William Goetzmann reprend tout cela dans un article publié en 2003, ” Fibonacci and the Financial Revolution” en ligne ici. Il faut aller pour cela dans le chapitre 12 du Liber Abaci, où Fibonacci parle de partage du profit (ce que nous appellerions l’allocation de capital, et que l’on retrouve dans les travaux sur les jeux coopératifs, en particulier les travaux de Shapley), mais aussi d’actualisationde flux futurs, et surtout, de réflexions sur la différence entre un taux annuel et un taux trimestriel. La version en anglaise est la suivante

ou encore

Dans la version originale, cela donne

Quand on pense qu’auparavant l’arithmétique se faisait en chiffre romain, et n’autorisait pas de raisonnements aussi subtils, on peut légitimement penser que Fibonacci est révolutionné le calcul actuariel (certes, en se contentant de traduire des textes arabes). Si on cherche à remonter encore dans le temps, on peut également citer un mathématicien indien, Brahmagupta, ब्रह्मगुप्त, vers 650, qui semble également avoir beaucoup inspiré Fibonacci.

De la vitesse de convergence, partie 1

Samedi après midi, à l’heure de la sieste pour certains, nous discutions avec Magali, Jean Baptiste et Florent des machines qui pouvaient trier les cartes, de là, la discussion a dévié sur les algorithmes de tri, et sur la difficulté d’expliquer l’importance de la vitesse d’un algorithme, ou plus généralement, la vitesse de convergence. Bref, avant de parler de vitesse de convergence en statistique – ou peut être avant sur les méthodes de calcul d’intégrale (ou d’espérance) – je voulais faire un petit billet rapide sur les algorithmes de tri. Je commence par des choses simples. De mémoire c’est des choses que je me souviens avoir vu au lycée dans les années 90, dans les cours d’informatique… Mais l’exemple peut être intéressant pour comprendre la différence, par exemple, entre https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-02.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-03.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-04.png.
Le premier algorithme est appelé tri à bulles. C’est probablement l’algorithme le plus simple, et sa vitesse est en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-01.png. On va balayer autant de fois la base que l’on a d’observations, et à chaque passage, on récupère le maximum (ou le minimum, peu importe). L’algorithme ressemble à ça,

> n=25; s=1
> X=sample(1:n,size=n,replace=FALSE)
>          for(i in 1:(n-1)){
+          for(j in 1:(n-1)){
+          s=s+1
+          if(X[j] > X[j + 1]){X[c(j,j+1)]=X[c(j+1,j)]}
+          }}

Si on l’implémente, voilà ce que ça donne (tous les graphiques partent de la même permutation initiale, et tous contiennent autant d’images, de l’ordre de https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-01.png)

Bref, un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png.Bon, j’ai été un peu brutal dans mon algorithme. En fait, on voit que les derniers points sont rapidement bien classés, donc inutile que la seconde boucle aille jusqu’au bout… on peut donc faire du https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-06.png, en faisant

 

+          for(j in 1:(n-i)){

Bref, la constante dans les https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png peut être relativement importante comme on le voit ici.

Mais on peut faire  mieux, voire beaucoup mieux comme on le verra ensuite.
Le second algorithme est le tri par insertion. C’est aussi un algorithme très simple:  on balaye l’ensemble en plaçant les valeurs rencontrées à leur place, un peu comme quand on tri ses cartes… L’algorithme est là aussi en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png. Ici l’algorithme converge un peu plus rapidement que le précédant, mais cela vient du fait que les données étaient partiellement ordonnées sur la fin. Autrement dit, un autre échantillon aurait pu donner autre chose… Mais généralement, c’est le cas.

>          for(i in 1:n){
+           x = X[i]
+          for(j in (i+1):2){
+          s=s+1
+         if(X[j - 1] < x){X[j] = X[j - 1];  j = j - 1;  X[j] = x}
+         }
+         }

(le tri est fait à l’envers ici…. mais c’est un détail).

Enfin, le tri fusion est un autre algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-07.png. L’idée est de couper, de trier, puis de fusionner. L’idée étant que trier un tableau dont la première moitié des éléments sont ordonnés, et où la dernière moitié est ordonné doit être plus rapide que quand tout est dans le désordre. A chaque étape, l’algorithme est simple, on compare le minimum des deux listes, puis on le vire à l’étape suivante. Sous R, ça donne, quand on commence

 

> for(i in 0:(trunc(n/2)-1)){
+ Y=X[c(1,2)+2*i]
+ if(Y[1]>Y[2]){X[c(1,2)+2*i]=X[c(2,1)+2*i]}
+ }

Bref, en commençant par classer les paires (1,2) puis (3,4), etc, puis les listes (1,2) versus (3,4). Le code est alors de la forme

> k=4
> for(i in 0:(trunc(n/4)-1)){
+ Y=X[1:k+k*i]
+ Y1=c(Y[1:(k/2)],rep(999999,k/2));Y2=c(Y[(k/2+1):k],rep(999999,k/2));
+ for(j in 1:k+k*i){
+     X[j]=min(c(Y1[1],Y2[1]),na.rm=TRUE)
+    if(X[j]==Y1[1]){Y1=Y1[
-1]}
+    if(X[j]==Y2[1]){Y2=Y2[-1]}
+ }}

(avec une petite correction pour le dernier bloc, qui ne compare pas forcément un bloc de taille 2 contre un autre de taille 2 (mais l’adaptation du code est infime). Puis on continue avec (1,2,3,4) versus (5,6,7,8), c’est relativement rapide de tout classer, comme le montre l’animation ci-dessous.

Il existe de très nombreux autres algorithmes de tri, le but étant de monter qu’un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-07.png va effectivement beaucoup plus vite qu’un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png. J’ai tout fait ici en R, mais pour les adeptes d’autres langages, il existe plein de code en ligne, ici par exemple. Pour conclure, le graphique ci-dessous montre les différences entre les vitesses de convergence en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-08.png (en noir), en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-09.png (en rouge) et en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-10..png (en bleu),

La suite sur la convergence des méthodes de monte carlo, puis des théorèmes limites en stats….

Discussions géométriques du matin (partie 2)

Suite rapide de mon billet sur la construction de triangles avec des bouts de bois (ici). Adrian me faisait remarquer () qu’il y avait plus simple pour calculer la probabilité de construire un triangle (je suis d’ailleurs affreusement jaloux de ne pas y avoir pensé tout seul). Mon principal soucis que j’ai en lisant le commentaire est “ça semble juste, mais est-ce que c’est juste ?” (mon expérience en TD de proba est que l’on peut souvent trouver un paquet de solutions qui “semblent” justes). Bref, faisons preuve de rigueur quelques minutes (même si France Inter continue à nous dire que le mot est toujours tabou ces temps-ci).
Rappelons que si on dispose d’un échantillon L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-01.png” ne peut être affichée car elle contient des erreurs. i.i.d. de loi https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-02.png, et de densité https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-03.png, et que l’on définie les statistiques d’ordre

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-04.png

la loi du couple https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-05.png, avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-06.png,est alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-07.png

à condition que https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-08.png (sinon la densité sera nulle). C’est jsute des histoires de compter les classements possibles, et puis de calculer quelques probas simples. Un cas particulier connu de tous est la loi du couple https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-09.png, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-10.png

Dans le cas qui nous intéresse, supposons que l’on tire deux bâtons, alors la loi du couple s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-11.png

alors que la densité du maximum était

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-12.png

donc la loi du minimum conditionnellement au maximum doit d’écrire, par la formule de Bayes

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-13.png

qui est effectivement la densité de la loi uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-14.png. Bref, ça marche pour deux bâtons.
Dans le cas de trois bâtons, la loi du triplet s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-15.png

alors que la loi du maximum devient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-30.png

aussi la loi conditionnelle du couple des deux premiers s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-31.png

Or si l’on avait deux lois uniformes sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-16.png, compte tenu de ce que l’on a écrit auparavant, la loi du couple des deux statistiques d’ordre seraient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-17.png

soit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-18.png

Bref, il semble qu’effectivement, le résultat soit correct: conditionnellement au maximum, les autres variables sont indépendantes, et suivent un loi uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-19.png.
Moralité, effectivement, les non-maximum sont indépendants, et leur loi est uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-16.png. On peut alors sommer ces deux lois uniformes et indépendantes, ce qui donne une loi triangulaire, sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-20.png, avec le sommet en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-21.png. Et par des raisons de symmétrie, la probabilité de dépasser https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-21.png est alors https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-22.png.Merci Adrian. C’est effectivement plus élégant que mon calcul moche…

De là à en faire un papier…

mais mon blog est l’endroit rêvé pour exposer mes interrogations récurentes. Je ne sais pas si vous avez remarqué, mais au toilettes, au bureau, il y a (généralement) du papier toilette dans un distributeur, et souvent, un rouleau qui traîne à côté….  et régulièrement, je me demande quelle peut être la stratégie optimale: prendre la papier du distributeur ou le papier du rouleau qui traîne…? car mine de rien, la personne chargée du réapprovisionement ne semble intervenir que si l’un des deux est vide.. autrement dit, si les deux diminuent trop vite ensemble, on risque d’avoir un problème un jour. Bon, ce qui me rassure, c’est que Donald Knuth avait remarqué des choses similaires, et avait eu la prétention de publier ses observations dans le American Mathematical Monthly. En fait, pour formaliser le problème, on note https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ02.png la probabilité que les gens utilisent du papier dans le rouleau le plus plein (on suppose connaître la quantité restant dans le distributeur qui dispose d’une fenêtre transparente), et https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ04.png la probabilité de piocher dans le plus petit rouleau. Et quand les deux font la même taille, on choisit au hasard (avec probabilités 1/2). On suppose que les rouleaux sont valables pour https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ60.png utilisateurs (je n’introduis pas de d’aléa, voire de dépendance en fonction du menu à la cantine). On parlera de portions. En gardant les notation de Knuth, on note https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ01.png l’espérance du nombre de portions restant sur le rouleau restant dès que l’un est terminé. Le but du papier de Knuth était d’étudier le comportement de https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ01.png quand https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ61.png. Mais comme ça n’a que peu d’intérêt, je trouve, on va regarder ce qui se passe à https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ60.png fixé.
Knuth note que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ05.png et que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ06.png. En fait, histoire d’épater les filles, Donald prétend que on peut obtenir facilement les résultats suivants

Pour le premier, c’est trivial, et pour le second, on peut faire un dessin (oui, je suis très dessin en ce moment),

A partir de la diagonale (points verts), il y a égalité entre les paquets, et donc on choisit avec probabilité 1/2. Pour calculer https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ62.png, on note qu’il y a deux issues possibles, soit il reste 1 portion sur le rouleau non vide (les points rouges) soit il reste 2 portions sur le rouleau non vide (les points bleus). Les probabilités de passer d’un point à l’autre sont indiquées sur le dessin. Bref, on a https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ03.png chances d’atteindre un point bleu, et https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ02.png chances d’atteindre un point rouge, dont l’espérance vaut ici

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ07.png

En continuant un peu, on obtient que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ08.png

En effet, sur notre petit dession, on rajoute juste les points autours, les valeurs à l’intérieur restant inchangées, et pour faire simple, je m’éloigne de la diagonale avec probabilité https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ03.png, et je m’en rapproche avec probabilité https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ02.png.

Comme tout à l’heure, on compte les probas d’arriver sur les différents points. On a https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ09.png chances d’atteindre un point violet, https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ10.png chances d’atteindre un point bleu, et enfin https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ11.png chances d’atteindre un point rouge (il suffit d’étudier toutes les trajectoires possibles). Au total, on arrive à

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ12.png

qui se simplifie en

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ13.png

Quant aux conditions de bord, on notera que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ14.png et que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ1.png. En fait, on peut étudier un problème plus général, mais plus réaliste: les paquets dans les distributeurs sont souvent beaucoup plus gros que celui qui traîne. On peut donc noter https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ16.png la même quantité, avec m le nombre de portions dans le distributeur et https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ60.png dans le rouleau qui traîne. Si on a compris le dessin ci-dessus, on se doute que ça ne change pas grande, chose, ormis que l’on ne commence pas sur la diagonale…
En fait, on note que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ17.png, avec les conditions de bord de la forme https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ18.png.
On peut aussi noter (il suffit de regarde sur le dessin en haut) que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ19.png

La preuve plus propre est que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ20.png

or comme https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ21.png par une symmétrie du problème, on a l’égalité. En utilisant l’idée de la preuve du résultat précédant (et le dessin, comme toujours), si https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ22.png, alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ23.png

Bref, on est capable de faire tous les calculs de manière récursive et relativement simple (pour ceux qui compte le faire la prochaine fois qu’ils vont aux toilettes, il faudra peut être emmener quand même leur ordinateur). On note que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ24.png

or https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ25.png (comme on l’a montré) et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ26.png

On a tout cette fois ci pour vérifier le calcul “facile” évoqué par Donald,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ27.png

Bref, on retombe sur le résultat précédant. On peut alors faire le calcul en fonction des valeurs de https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ02.png (que l’on fixera, puis de m et https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ60.png),

> p=0.4
> q=1-p
> m=50
> n=50
> M=matrix(NA,m+1,n+1)
> M[,1]=seq(0,m)
> M[1,]=seq(0,n)
> M[2,2]=2-p
> for(i in 3:(m+1)){
+ for(j in 2:(i-1)){
+ M[i,j]=p*M[i-1,j]+q*M[i,j-1]
+ }
+ M[i,i]=M[i,i-1]
+ }
> for(i in 2:(m)){
+ for(j in (i+1):(n+1)){
+ M[i,j]=M[j,i]
+ }}
> colnames(M)=seq(0,n)
> rownames(M)=seq(0,m)
> image(seq(0,n),seq(0,n),M,)
> contour(seq(0,n),seq(0,n),M,add=TRUE,
+ levels=seq(0,50,by=5))

Voilà ce que l’on obtient si https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ28.png. le nombre de portions qui reste dès que l’un des rouleaux est fini peut se visualiser sur le graphique suivant, la couleur rouge correspond aux cas un peu dangereux où il reste peu de feuilles dans le rouleau, les zones jaunes indiquant au contraire des niveaux plus élevés,

Avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ29.png, on obtient

Autrement dit, il reste toujours du papier sur l’autre dès lors que les gens pensent à alterner suffisement (ici 95% des utilisateurs) alors que pour https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ30.png,

autrement dit, on risque constament de se retrouver à court… Ce qui est marrant (?) c’est le cas où https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ31.png: s’il y autant de personnes qui piochent dans le gros paquet ou dans le petit, alors il y a de fortes chances pour que les deux rouleaux soient vides dès que l’un l’est…

[ce billet a été initialement posté sur freakonometrics.blog.free.fr/… en mai 2010]

Moyenne géométrique versus moyenne arithmétique…

L’autre jour, j’étais de surveillance d’examen. Pour passer le temps, j’ai survolé l’énoncé, et je me souviens avoir lu “montrez que la moyenne arithmétique est plus grande que la moyenne géométrique“. Ayant généralement tendance à préférer les copies où il y a peu à lire, mais avec de jolis dessins, j’ai eu envie de proposer une correction sans écrire un mot aux étudiants (qui semblaient souffrir dans leur démonstration). Il suffit de monter que

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/moyennes-latex.png

Ben voilà,

quod erat demonstrandum. Non ? On note en plus que l’on a égalité si et seulement si les nombres sont égaux…. Allez, pour convaincre les sceptiques, je rajoute une deuxième preuve… Bon, en fait je peux prouver que

mais c’est la même chose, non ?

Discussions géométriques du matin (partie 1)

Petit résumé du débat du matin avec mon fils …
– papa, papa j’ai trouvé trois bouts de bois,
– ben… fais une figure géométrique avec *
– ça marche pas !
– comment ça, “ça marche pas” ?!?
– j’arrive pas à faire un triangle !
– … hein ? attends, j’arrive
genre “mais qu’est-ce que tu racontes ? t’es pas réveillé ce matin ou quoi ?“. Bref, après avoir vu qu’effectivement ça n’était pas possible (un des bouts de bois était beaucoup plus long que les deux autres),  je me suis demandé si mon fils n’avait vraiment pas eu de chance… ou si au contraire, c’était plutôt “normal” de ne pas pouvoir faire un triangle à partir de trois bouts de bois….
Soient https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-01.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-02.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-03.png les longueurs des trois morceaux de bois (oui, une fois le café avalé, je peux prendre un peu de temps pour formaliser). Ces trois bouts de bois forment  un triangle si et seulement si la somme des deux plus petits est plus grande que le plus grand. Autrement dit, si je m’intéresse aux statistiques d’ordre

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-04.png

il faut que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-05.png

Pour faire simple, on supposera que les https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-50.png suivent une loi uniforme sur [0,1].

  • Approche par simulation

Par simulation, ça ne devrait pas être trop compliqué d’avoir un ordre de grandeur de la probabilité

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-06.png

Le code est tout simplement

> n=1000000
> U=matrix(runif(3*n),n,3)
> U1=apply(U,1,max)
> U2=apply(U,1,sum)-U1
> U3=U1<U2
> mean(U3)
[1] 0.501067

Finalement la probabilité semble être proche de 1/2… ça doit pouvoir se montrer….

  • Approche par des calculs de proba

Il doit être possible de calculer directement

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-06.png

ou bien, on peut aussi noter que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-07.png

ce qui ne fait intervenir que deux variables (certes, toujours non indépendantes, mais on devrait s’en sortir malgré tout): le maximum des trois variables, et la somme des trois variables. On cherche alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-08.png

Je pense qu’il est plus simple de travailler sur la première expression, en notant qu’il existe 6 permutations de ces 3 valeurs, et cette probabilité s’écrit simplement

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-10.png

(qui  est un des 6 arrangements possibles, tous équiprobables). Cette probabilité est alors tout simplement, compte tenu de l’indépendance des tailles des bouts de bois, ce qui est peut être une hypothèse un peu forte, certes,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-11.png

qui doit pouvoir se réécrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-12.png

L’intégrale au centre se réécrit alors simplement, à condition de savoir comment https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-13.png se positionne par rapport à 1. Deux cas assez différents sont à distinguer: comme on veut https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-134.png, soit https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-15.png, soit https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-16.png. Autrement dit on écrira l’ensemble sur lequel on calcule la première intégrale comme la réunion de deux ensembles (disjoints),

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-18.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-19.png

Sur le premier ensemble,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-20.png

et sur le second

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-21.pngBref, l’intégrale se retrouve sous la forme de deux intégrales

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-22.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-23.png

La première intégrale correspond à un cas où https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-24.png est petit, ou plutôt https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-25.png. La seconde intégrale est un cas où https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-26.png est grand, ou plutôt https://perso.univ-rennes1.fr/arthur.charpentier/latex/b<br />
ois-27.png. Bref, la première intégrale peut s’écrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-28.png

et la seconde

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-29.png

(on aura en mémoire le facteur 6, que je ne mes pas pour alléger un peu). Bref, on a deux fois la même intégrale, qui s’écrire (dans le second cas à une translation des indices près)

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-30.png

Notre probabilité est alors 6 fois le double de cette intégrale, soit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-31.png

Bon, finalement c’est beaucoup plus improbable de pouvoir construire un triangle que je ne le croyais au départ…. Damned ! Moralité, mieux vaut ne pas trop réfléchir avant d’avoir pris un café…. Et comme “papa” est très fort en mauvaise foi, “papa” expliquera que si la longueur des bois de bois n’était pas distribuée suivant une loi uniforme, l’intuition initiale n’était pas complètement stupide… Non mais !
* pour ma défense, c’était ça ou la révision de la table de multiplication de 5… en gros, “papa a pas encore pris son café, c’est pas le moment“.

Acheter des pochettes Pokemon ? et puis quoi encore !

Il y a quelque temps, Olivier Gergaud était venu à Rennes présenter son papier Stardom and Perceived Talent.Nous en avions profité pour discuter un peu d’un ancien papier qu’il avait fait sur les Pokemon (ici). Il faut avouer que je suis un peu jaloux car il essaye de répondre à une question que je me posais: le prix des cartes Pokemon a-t-il un lien avec les “pouvoirs” de la carte ou sa rareté.

Mais aujourd’hui je voulais répondre à une question simple que je me pose depuis que mon fils me dit (sans grand succès) “Papa, achète moi des cartes Pokemon” (ou “Naruto“, ou je ne sais quoi encore). Mon argument, en tant que père rationnel, est que ces cartes couchent très chers. Aujourd’hui, je vais le prouver. Supposons qu’il existe 200 cartes distinctes, et que ces cartes sont vendus 5 euros le paquet de 10. Pour simplifier, supposons que ces cartes soient vendues individuellement dans des cachés hermétiques à 50 centimes pièce. On va supposer que les cartes soient équi-réparties. Supposons que je puisse acheter une vignette spécifique (par exemple en ligne), dont je connais le personnage, pour 1 euro (soit le double du prix si le personnage est inconnu à l’achat). Quelle stratégie doit-je adopter ? et combien ça va me coûter ?
Je note https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-17.png le nombre de carte total (ici 200 car je suppose qu’il existe 200 personnages). Je note https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-01.png le coût d’une carte cachée (50 centimes ici) et https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-02.png le coût d’une carte non-cachée (1 euro). Supposons que j’ai à en ma possession les cartes https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-03.png (et que je n’ai pas les suivantes). Si je veux récupérer une nouvelle carte (que je noterais https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-05.png), en achetant au hasard i.e. dans des sachets hermétiques, le nombre de cartes à acheter est https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-06.png, qui suit une loi géométrique, de paramètre https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-07.png, où

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-08.png

Par exemple, si j’ai https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-09.pngcartes, et qu’il m’en manque simplement une, alors la probabilité de la tirer est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-10.png

Pour revenir au tirage d’une nouvelle carte que je n’avais pas (alors que j’en https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-04.png), le coût moyen de ma stratégie est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-11.png

J’ai donc intérêt à acheter cette carte directement dès lors que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-12.png

ce qui devient intéressant dès lors que

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-13.png

Dans mon exemple numérique, https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-02.png étant le double de https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-01.png, je commence à acheter directement dès lors que j’ai en ma possession 100 cartes distinctes. Le coût moyen est alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-14.png

Aussi, dans le cas que j’évoquais, on obtient 171€,
> .5*sum(200/(200-(1:101)-1))+1*(200-101+1)
[1] 171.0933

Je peux comparer à deux stratégies que l’on connaît. La première est de tout acheter à un collectionneur (i.e. à découvert, ce qui limite l’intérêt du “jeu“, j’en conviens). Dans ce cas, le coût total sera

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-21.png

soit 200€. L’autre stratégie est d’acheter continuellement des paquets clos chez mon buraliste. Dans ce cas, le coût moyen serait

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-16.png

Le calcul exact donne ici 587,80€, comme  le montre le calcul ci-dessous,
> 0.5*200*sum(1/1:200)
[1] 587.8031

mais on peut aussi utiliser une approximation de la forme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-22.png

(ce qui donnerait 529€ environ). Pour ceux qui ont besoin de voir pour croire, j’ai représenté ci-dessous le prix en fonction du moment où je décide de passer des cartes fermées aux cartes visibles.

Ça y est, j’ai enfin un argument rationnel à faire valoir à mon fils…. même si je ne suis pas certain de réussir à le convaincre aussi facilement…

Optimal control, part 2

In the first part (here), we introduced Bellman’s idea of backward induction. But what if we consider now infinite time horizon ? Actually, the maths will be even more simple… and we will be able to use fixed pointed theorem to derive solutions.

  • The mathematical framework

Here, consider the following value function

https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-01.png

and define

https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-02.png

A sequence https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-04.png is said to be an admissible solution for starting point x,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-05.png

if https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-06.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-07.png. If we reformulate the dynamic programming idea, we obtain that if https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-04.png is a solution to problem https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-08.png, then for all https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-09.png, sequence https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-10.png is a solution to problem https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-11.png. It comes that function v is a solution of Bellman’s equation

https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-12.png

Note that is can be ssen as some fixed point resul, since

https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-13.png

i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/opt-contr-14.png

So far, it shouldn’t be so hard….

  • Frank Ramsey’s model (discrete version)

In 1928, Frank Ramsey wanted to understand the amount of savings in a dynamic perspective (in how much of its income should a nation save). Consider the following infinite horizon problem, where some planifier wants to maximize

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-01.png

subject to constraints https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-02.png https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-03.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-04.png.
Before looking at dynamic programming answers, we might start with standard Lagrangian optimization techniques.
Assuming concavity of utility function and production function, we should look only for interior solutions. Define the Lagrangian as

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-05.png

Thus, the first order conditions are then given by

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-06.png

and

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-088.png

Assume further some terminal condition, e.g.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-08.png

(also called transversality condition). If we combine those two conditions, and assume that the first constraint is saturated, we obtain the so-called Euler equation,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-09.png

It is also possible to use Bellman’s equation: given the dynamic of the capital

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-10.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-11.png

The first order condition states

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-12.png

But since v is unknown, so is its derivative. But from the enveloppe theroem, we obtain something like

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-13.png

where

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-14.png

We can then write

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-15.png

i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-16.png

and finally

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-17.png

which is, Euler’s equation.

  • A specified model, with calculations

As in the previous post (here), consider a log utility function, and a power production function,https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-18.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-19.png. The dynamic is then

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-20.pngand https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-21.png

Note that fixed points are here

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-22.png

and

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-23.png

Recall that the value function is defined as

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-24.png

A natural idea to derive the value function can be to iterate, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-25.png

starting with a simple function, e.g the null function, at step 0. At step n=1

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-26.png

thus

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-27.png

At step n=2,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-28.png

i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-29.png

The first order condition is then

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-30.png

and thus, we obtain

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-31.png

that can be plugged in the previous equation, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-32.png

At step 3, we start from that new expression, derive the first order condition, and we get

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-33.png

and

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-34.png

and so on…
And finally, we can prove that

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-35.png

i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-36.png. Assuming that

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-37.png

actually, we can prove that

https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-38.png

(and https://perso.univ-rennes1.fr/arthur.charpentier/latex/co-opt-39.png has a form that can be explicited).

Optimal control, part 1

Since I had recently a request (from Benoît) about optimal control, I will start here a series of posts on that topic. But first, let us start with a simple problem, with discrete time, no randomness, and with a finite horizon. This might be a too simple framework to model complex problem, but that should be interesting to derive heuristic intuitions (I will skip here the mathematical problems, that can be found in very good books… references will come soon).

  • An introduction to backward induction

Before starting seriously, let us consider the following example: we want to reach the red city on the right from the red city on the left, as fast as possible. There are some roads, and the number is the number of hours it takes. Let us prove that the optimal way is the red one,

The first idea can be to calculate all possible trajectories, but with a large number of roads, the number of possible ways can soon be extremely large. An alternative can be to look backwards (like any students facing a question where the answer is given: start from the end, and try to find a possible way to reach it).


Numbers are greens are the number of hours that we still need one we’ve reached that point. Let us move again one step backward, and consider the orange points,

In the middle, we had to chose either to go up (it will still take 9 hours) or go down (and then 14 hours are necessary). Thus, the optimal strategy, once we’ve reach that point, it to take the 9 hour road. This will be idea the idea proposed by Bellman. Let us go backward again, to the purple cities. Again, we have to chose the shorter way,

From the top, the fastest road will take 13 hours, and from the bottom, it will take 16 hours. Thus, since it takes 10 hours to reach the top city, this road is necessarily faster than taking the one below (since it takes 8 hours to reach the nearest city).
Any this is it. We now have an intuitive idea of what should be done to find optimal strategies.

  • The optimization problem, discrete with finite horizon

Let us consider the following optimization problem: we want to find the optimal strategy https://freakonometrics.hypotheses.org/files/2022/04/control-01.png that maximizes the following function

https://freakonometrics.hypotheses.org/files/2022/04/control-02.png

with simple constraints, such as https://freakonometrics.hypotheses.org/files/2022/04/control-03.png (a state space) and https://freakonometrics.hypotheses.org/files/2022/04/control-04.png (i.e. some dynamic constraints). Assume further that the starting point is given, i.e.
https://freakonometrics.hypotheses.org/files/2022/04/control-05.png. In economic application, note that frequently https://freakonometrics.hypotheses.org/files/2022/04/control-06.png i.e. we consider a discounted version of the value.

  • The idea of dynamic programming

The intuitive idea of dynamic programming is that if the optimal path from A to C goes throught B, then the path is optimal from B to C. Thus, it will be natural to consider backward induction techniques.
Thus, define
https://freakonometrics.hypotheses.org/files/2022/04/control-07.png
https://freakonometrics.hypotheses.org/files/2022/04/control-08.png
https://freakonometrics.hypotheses.org/files/2022/04/control-09.png

https://freakonometrics.hypotheses.org/files/2022/04/control-10.png
https://freakonometrics.hypotheses.org/files/2022/04/control-11.png
Then Bellman’s principle can be used to link those problems: if https://freakonometrics.hypotheses.org/files/2022/04/control-12.png is a solution of the problem https://freakonometrics.hypotheses.org/files/2022/04/control-13.png then, for all
https://freakonometrics.hypotheses.org/files/2022/04/control-14.png, https://freakonometrics.hypotheses.org/files/2022/04/control-15.png is a solution of problem https://freakonometrics.hypotheses.org/files/2022/04/control-16.png.
Note that, so far, we assume that such an optimal sequence does exist. Thus, we get that for all x

https://freakonometrics.hypotheses.org/files/2022/04/control-17.png

and more generally,

https://freakonometrics.hypotheses.org/files/2022/04/control-18.png

Hence, from a practical point of view, we solve those equation using a backward approch. I.e. first,

https://freakonometrics.hypotheses.org/files/2022/04/control-20.png

and then

https://freakonometrics.hypotheses.org/files/2022/04/control-21.png

and so on

https://freakonometrics.hypotheses.org/files/2022/04/control-22.png

… etc. It can be proved that the sequence https://freakonometrics.hypotheses.org/files/2022/04/control-23.png is solution of https://freakonometrics.hypotheses.org/files/2022/04/control-24.png if and only if for all https://freakonometrics.hypotheses.org/files/2022/04/control-25.png, https://freakonometrics.hypotheses.org/files/2022/04/control-15.png
is solution of

https://freakonometrics.hypotheses.org/files/2022/04/control-27.png

So far, it does not look so difficult….

  • A simple example

Let https://freakonometrics.hypotheses.org/files/2022/04/control-30.png denote the consumption at period t, and assume consumption yields utility

https://freakonometrics.hypotheses.org/files/2022/04/control-31.png

as long as the consumer lives. Assume the consumer is impatient, or has a stronger preference for present, so that he discounts future utility by a factor https://freakonometrics.hypotheses.org/files/2022/04/control-32.png. Let https://freakonometrics.hypotheses.org/files/2022/04/control-33.png be capital he got at time t. Assume that his initial capital is a given amount https://freakonometrics.hypotheses.org/files/2022/04/control-34.png, and suppose that this period’s capital and consumption determine next period’s capital as

https://freakonometrics.hypotheses.org/files/2022/04/control-35.png

where https://freakonometrics.hypotheses.org/files/2022/04/control-36.png is a positive constant and https://freakonometrics.hypotheses.org/files/2022/04/control-37.png. Assume further that capital cannot be negative. Then the consumer’s problem is simply

https://freakonometrics.hypotheses.org/files/2022/04/control-38.png

given https://freakonometrics.hypotheses.org/files/2022/04/control-39.png. Bellman’s equation is then

https://freakonometrics.hypotheses.org/files/2022/04/control-40B.png

whih leads us to a simplier problem than the initial one, since only two variables are involved here https://freakonometrics.hypotheses.org/files/2022/04/control-41.png and https://freakonometrics.hypotheses.org/files/2022/04/control-42.png. And to solve that problem, we use backwards induction techniques.
Since https://freakonometrics.hypotheses.org/files/2022/04/control-44.png is known, we can derive easily https://freakonometrics.hypotheses.org/files/2022/04/control-45.png, and so on until https://freakonometrics.hypotheses.org/files/2022/04/control-46.png. More precisely, given https://freakonometrics.hypotheses.org/files/2022/04/control-47.png, we can get https://freakonometrics.hypotheses.org/files/2022/04/control-49.png which is the maximum of function

https://freakonometrics.hypotheses.org/files/2022/04/control-50.png

with https://freakonometrics.hypotheses.org/files/2022/04/control-51.png. One can see that the following function is a possible solution

https://freakonometrics.hypotheses.org/files/2022/04/control-52.png

where each https://freakonometrics.hypotheses.org/files/2022/04/control-53.png is a constant. Further, the optimal amount to consume at time https://freakonometrics.hypotheses.org/files/2022/04/control-54.png is

https://freakonometrics.hypotheses.org/files/2022/04/control-55.png

i.e., if we explicit those expressions

https://freakonometrics.hypotheses.org/files/2022/04/control-60.png
https://freakonometrics.hypotheses.org/files/2022/04/control-61.png
https://freakonometrics.hypotheses.org/files/2022/04/control-62.png

…etc.

  • A reformulation of the optimisation problem

Another way of expressing the optimization problem is the following: https://freakonometrics.hypotheses.org/files/2022/04/control-70.png was the variable of interest, https://perso.univ-rennes1.fr/arthur.charpenti<br /> er/latex/control-71.png was the control variable, and https://freakonometrics.hypotheses.org/files/2022/04/control-72.png. Thus, the programm

https://freakonometrics.hypotheses.org/files/2022/04/control-73b.png

can be expressed as

https://freakonometrics.hypotheses.org/files/2022/04/control-74.png

Several extension can be considered,

  • consider an infinite sum
https://freakonometrics.hypotheses.org/files/2022/04/control-80.png
  • consider a random component
https://freakonometrics.hypotheses.org/files/2022/04/control-81.png
  • consider a continuous version
https://freakonometrics.hypotheses.org/files/2022/04/control-82.png

But those items will be for posts that I still have to write down….

Quand peu d’information donne énormément d’information

Une petite histoire empruntée à Alfred (et Nicolas1) suite au repas de mardi, à l’X. Un explorateur débarque un jour sur une petite île, perdue au milieu des mers. Cette île est un peu particulière: elle est peuplée de personnes ayant soit les yeux bleu, soit les yeux noirs (et tout le monde sait qu’il n’y a pas d’autre alternative). Et sous aucun prétexte personne ne doit connaître la couleur de ses propres yeux, sinon la légende dit qu’il doit se suicider avant le coucher du soleil (du jour où il l’apprend). Et s’ils ne le font pas, ils se retrouvent maudit jusqu’à la huitième génération. Cette île vivait en parfaite harmonie depuis plusieurs siècles, lorsque un explorateur débarqua. Il y avait alors 60 habitants sur cette île.

Bref, comme le veut la coutume, notre explorateur doit tenir un discours devant ses hôtes. Ne sachant trop que dire, il ne trouve rien d’autre à dire que de s’étonner de voir à l’autre bout du monde un île avec des habitants aux yeux bleus comme les siens (la légende lui avait racontée auparavant, de manière à éviter toute gaffe).
Deux mois plus tard, tous les habitants de l’île sont morts !
Quoi de plus étrange n’est-ce pas… ? Il n’a pourtant pas dit grand chose, si ce n’est une trivialité, à savoir qu’il y a des personnes aux yeux bleus sur l’île (ce que tout le monde devait savoir car s’ils ne connaissaient pas leur propre couleur d’yeux, ils pouvaient observer celle des autres). Pourtant cette micro-information semble avoir suffit à ce que les habitants en déduisent la couleur de leurs yeux !
En effet, le lendemain du discours, tout le monde s’est fait la réflexion suivante “d’après le discours de l’explorateur, il y a au moins une personne aux yeux bleus sur l’île. Or s’il était seul, il se serait rendu compte que tout le monde a des yeux noirs, et donc il aurait du se suicider dans la nuit. S’il ne s’est rien passé, c’est qu’il y a au moins deux personnes dans l’île avec des yeux bleu. Si quelqu’un s’est suicidé, c’est que j’ai les yeux noirs, et je dois me suicider avant le coucher du soleil“. Le deuxième jour, si tout le monde est encore en vie, tout le monde se fait la réflexion suivante “comme personne ne s’est suicidé, il y a au moins deux personnes aux yeux bleus sur l’île. Mais s’ils n’étaient que deux, chacun devrait penser que l’autre est le seul, et donc maintenant il sait qu’il est le second. Ils aurait du se suicider dans la nuit. S’il ne s’est rien passé, c’est qu’il y a au moins trois personnes dans l’île avec des yeux bleu. Si deux personnes se sont suicidé, c’est que j’ai les yeux noirs, et je dois me suicider avant le coucher du soleil“. Etc (je dois continuer ?). Bref, au plus tard, on arrive 60 jours après le discours du navigateur, et si tout le monde est encore en vie, c’est que  tout le monde se fait la réflexion suivante “comme personne ne s’est suicidé, il y a au moins soixante personnes aux yeux bleus sur l’île, c’est à dire que tout le monde a les yeux bleu. Je dois donc me suicider ce soir“.
https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso/desproges.jpgPourtant l’explorateur n’a pas franchement donné d’information… et pourtant cette information a suffi à décimer l’île ! Étonnant non ?
1 Nicolas avait raconté une version beaucoup plus connue, que l’on peut retrouvé dans tous les classiques sur les énigmes sous le titre des cocus de Bagdad. L’histoire est assez proche, avec un calife à Bagdad qui, jugeant immoral l’adultère dicte une loi qui stipule que toute femme convaincue d’adultère sera étranglée par son propre mari le lendemain matin même de la découverte par le  mari de la faute. Il faut ajouter que tous les maris voient si leurs voisins sont cocus ou pas, mais ils sont incapables de voir si leur femme les trompe… Et le calife a eu la mauvaise idée de dire que certains hommes étaient cocus… moralité, au bout de quelques jours toutes les femmes qui ont trompé leur mari finissent égorgées. Moralité, il n’y a pas besoin de beaucoup d’information pour agir de manière efficiente…

Théorème de l’enveloppe

Un tout petit complément sur l’utilité inverse, et l’exercice 2 de la PC de jeudi après midi.
Tout  d’abord un petit rappel sur le théorème de l’enveloppe. Supposons que l’on cherche à maximiser une fonction (dite objectif) f(bold x,bold r) qui dépend de paramètres bold r, i.e. on cherche

f^*(bold r) = max_{bold x} f(bold x,bold r),

La fonction f^*(bold r) est alors la fonction optimale pour les paramètres bold r donnés. Soit  bold x^*(bold r) la valeur optimale en bold x, fonction du paramètre bold r , de telle sorte que

f^*(bold r) = f(bold x^*(bold r), bold r)

Le théorème de l’enveloppe nous permet de voir comment f^*(bold r) varie lorsque le paramètre change, et plus particulièrement, on a

frac{partial f^*(bold r)}{partial r_i} = frac{partial f(bold x,bold r)}{ partial r_i} Bigg|_{bold x =

Autrement dit, la dérivée de f^*(bold x,bold r) suivant  ri est donnée par la dérivée partielle de f(bold x,bold r)suivant ri, lorsqu’on l’évalue en bold x =.

Notons qu’il existe aussi une version un peu plus générale qui autorise un programme d’optimisation avec contrainte.
Supposons que l’on s’intéresse à

max_{bold x} f(bold x,bold r) ;; s.t. ;; bold g(bold x,bold r) = bold 0

dont le Lagrangien s’écrit

mathcal{L}(bold x,bold r) = f(bold x,bold r) - boldsymbol{lambda} cdot bold g(bold x,bold r)

où de manière classique boldsymbol{lambda} = (lambda_{1},dots,lambda_{n}) est l’ensemble de multiplicateurs de Lagrange où

bold g(bold x,bold r) = (g_{1}(bold x,bold r),dots,g_{n}(bold x,bold r))

Le théorème de l’enveloppe permet d’écrire

frac{partial f^*(bold r)}{partial r_i} = frac{partial mathcal{L}(bold x,bold r)}{partial r_i} Bigg|_{ bold x =

L’autre point important de la PC était le principe de dualité dans les programme d’optimisation, que j’avais abordé ici.
Bref, tout ça permet d’introduire la notion d’utilité inverse et de demande Hicksienne. Pour reprendre calmement cette partie (en s’inspirant d’une page de wikiversité sur le sujet), le programme classique visait à chercher

 max_{mathbf{x}in Bleft(mathbf{p},y<br />
ight)}uleft(mathbf{x}<br />
ight)

où la contrainte budgétaire s’écrit classiquement

B=left{mathbf{x}in X|mathbf{p}mathbf{x}=sum_{i=1}^Np_ix_i le y<br />
ight}

La solution est alors notée mathbf{x}^mleft(mathbf{p},y<br />
ight) (on parle de demande Marshalienne). Je passe rapidement sur le programme d’optimisation de l’agent, mais nous avions vu que le taux marginal de substitution était égal au rapport des prix. En effet, le Lagrangien s’écrit

 mathcal{L}equiv u(mathbf{x})-lambda(mathbf{p}mathbf{x}-y)

ce qui donne de conditions du premier ordre de la forme

 frac{partial u(mathbf{x}^m)}{partial x_i}=u'_i=lambda p_i forall i.
L’idée maintenant est d’injecter les fonction de demande dans la fonction d’utilité, ce qui permet d’obtenir la fonction d’utilité inverse

 v(mathbf{p},y)=u(mathbf{x}^m(mathbf{p},y))

Une manière alternative d’étudier ce problème est de regarder le programme dual (on a supposé la fonction d’utilité convexe), i.e on cherche à trouver le panier de bien le moins cher possible nous procurant une utilité au moins égale à une utilité minimale, i.e. on cherche à résoudre

 min_{mathbf{x}|u(mathbf{x})gebar u}mathbf{p}mathbf{x}

La solution de ce programme d’optimisation mathbf{x}^h(mathbf{p},bar u) est appelée demande hicksienne. On peut alors considérer la fonction de dépense minimale,

e(mathbf{p},bar u)=mathbf{p}mathbf{x}^h(mathbf{p},bar u)

Le lemme de Shepard permet d’écrire

: x^h_i(mathbf{p},bar u)=frac{partial e(mathbf{p},bar u)}{partial p_i}

Ce résultat se montre finalement assez simplement: en effet, par construction

e(mathbf{p},bar u)=sum_{i=1}^Np_ix_i^h(mathbf{p},bar u)

dont un calcul de dérivée rapide permet d’écrite

 frac{partial e(mathbf{p},bar u)}{partial p_k}={color{Red}sum_{i=1}^Np_ifrac{partial x_i^h(mathbf{p},bar u)}{partial p_k}}+x_k^h(mathbf{p},bar u)

Rappelons que la condition de premier ordre s’écrit

u'_i=frac{p_i}{lambda}

Comme en u(mathbf{x})=bar u, nous sommes à un optimum, on en déduit

 sum_{i=1}^Nu'_ifrac{partial x_i^h(mathbf{p},bar u)}{partial p_k}=0

et donc

frac{1}{lambda}{color{Red}sum_{i=1}^Np_ifrac{partial x_i^h(mathbf{p},bar u)}{partial p_k}}=0

et comme lambda e 0 c’est que  la somme rouge est nulle. D’où le lemme dit de Shepard.
Les théorèmes de dualité permettent d’obtenir un certain nombre de résultats, parmi lesquels

  • e[mathbf{p},v(mathbf{p},y)]equiv y
  • v[mathbf{p},e(mathbf{p},bar u)]equivbar u
  • x^m_i(mathbf{p},y)equiv x^h_i[mathbf{p},v(mathbf{p},y)]
  • x^h_i(mathbf{p},bar u)equiv x^m_i[mathbf{p},e(mathbf{p},bar u)]

Enfin, en appliquant le théorème de Slutsky, on écrit

ce qui s’écrit encore

frac{partial x_k^m}{partial p_j}=frac{partial x_k^h}{partial p_j}-frac{partial x_k^m}{partial y}{color{Red}frac{partial e}{partial p_j}}

Bref, on retrouve ici ce qu’on avait fait en PC.
Mais si ça n’est toujours pas plus clair, il y a des notes de cours de Michael Manove, à Boston University, sur “indirect utility” et “duality“, dont les slides sont ici. Il y a également des compléments mathématiques, dont le théorème de l’enveloppe (Prop. 21.2) ou le théorème de Kuhn-Tucker (Prop.21.13), ici.
On pourra aussi trouver en ligne l’article de Lawrence Blume sur la dualité pour le dictionnaire Palgrave d’économie, ici. Larry a aussi écrit l’article sur la programmation convexe, diponible .
Sinon je mets ici les exercices des deux prochaines PC.

Economie, et optimisation sous contraintes

Premier cours de microéconomie à l’X. Pour les PC (salle PC10), la liste des exercices est disponible ici, avec au programme la consommation des ménages. Je mets ici les exos de la prochaine PC (et une autre feuille d’exercices – que nous ne traiterons pas).

Petit complément sur la résolution du programme d’optimisation avec contrainte, via le Lagrangien. Comme je l’ai évoqué (très) rapidement, c’est lié à des résultats connus sous le nom de théorème (ou condition) de Kuhn et Tucker. De manière assez générale, le but est de trouver X^* in mathbb{R}^n qui maximise f(X) , sous contrainte g(X) ge 0, c’est-à-dire résoudre le programme
X^*in arg max_{X in mathbb{R}^n} f(X) sous contrainte g(X) ge 0g(X)=(g_i(X))_{1 le i le m} (correspondant à un ensemble de m contraintes).

Comme précisé en cours, l’idée de la résolution du problème (sous des hypothèses – fortes – de concavité des fonctions) est d’introduire des multiplicateurs de Lagrange (autant que de contraintes) lambda = (lambda_i)_{1 le i le m} in mathbb{R}_+^m , et on définit le Lagrangien comme la fonction 
L(X, lambda) = f(X) + sum_{i=1}^m lambda_i g_i(X)

On considère alors le problème d’optimisation (sans contrainte) du Lagrangien. Les conditions du premier ordre donne n premières conditions, à savoir

frac{partial f}{partial X}(X) + sum_{i=1}^m lambda_i frac{partial g_i}{partial X}(X)=0

et les m dernières correspondent à une saturation des contraintes (en fait  g_i(X)>0 Rightarrow lambda_i = 0, donc si la contrainte n’est pas saturée, le multiplicateur de Lagrange est nul). Ces n+m conditions sont appelées “conditions de Kuhn-Tucker“. L’intérêt de cette formalisation est que l’on a le résultat suivant:

X^*in arg max_{X in mathbb{R}^n} f(X) sous les contraintes g(X) ge 0 si et seulement si (X^*, lambda^*) , est solution des conditions de Kuhn-Tucker.

Pour plus plus d’information, je peux renvoyer aux slides de Christophe Gonzales àParis 6, ici, ou aux notes de Julien Grenet, de l’ENS, . Sinon j’avais fait un billet il y a quelque temps sur un sujet connexe (et qu’on verra apparaître dans les cours suivants normalement), i.e. le principe de dualité en optimisation (ici).

Paradoxe et enseignement des probabilités

L’été semble être propice aux billets sur les paradoxes et autres bizarreries intellectuelles (je renvoie par exemple au très bon billet de Laurent, ici). Je vais en rajouter une couche (j’avais déjà parlé du paradoxe des enveloppes, ici ou de Simpson, ). En fait, je voulais revenir sur un paradoxe que j’avais lu en Terminale (à l’époque il y avait une introduction aux probabilités). J’y repense car je suis en train de finir de préparer mon cours de rappels de probabilité, et je cherche des exercices “simples” évitant toute ambiguité afin d’être un peu éclairant… pas comme ici….

Lançons trois pièces de monnaie. Quelle est la probabilité que les 3 pièces tombent du même côté (peu importe lequel) ?

Réponse 1: facile, faisons un peu de dénombrement, il y a 8 cas possibles (PPP, FPP, PFP, PPF, FFP, FPF, PFF, FFF), qui sont équiprobables, et seuls deux nous intéressent (FFF et PPP), on a alors 1 chance sur 4.

Réponse 2: facile, raisonnons deux minutes et essayons de trouver une réponse élégante. Comme j’ai 3 pièces pour seulement deux faces, il y en a forcément deux qui sont tombées du même côté, et pour la troisième, il me reste une chance sur deux, et la réponse est alors 1 chance sur 2.

Etonnant, non?“.

Bon, c’est un exemple assez classique en théorie des probabilités, je laisse aux lecteurs le temps de la réflexion…. le but étant de montrer qu’il y a une faute dans le raisonnement de la personne ayant donné la seconde réponse… Oui, c’est toute la difficulté des TD de probas: comment expliquer qu’un raisonnement qui semble juste est falacieux…