Recently, one of my students asked me about optimization routines in R. He told me he that R performed well on the estimation of a time series model with different regimes, while he had trouble with a (simple) GARCH process, and he was wondering if R was good in optimization routines. Actually, I always thought that mixtures (and regimes) was something difficult to estimate, so I was a bit surprised…

Indeed, it reminded me some trouble I experienced once, while I was talking about maximum likelihooh estimation, for non standard distribution, i.e. when optimization had to be done on the log likelihood function. And even when generating nice samples, giving appropriate initial values (actually the true value used in random generation), each time I tried to optimize my log likelihood, it failed. So I decided to play a little bit with standard optimization functions, to see which one performed better when trying to estimate mixture parameter (from a mixture based sample). Here, I generate a mixture of two gaussian distributions, and I would like to see how different the mean should be to have a high probability to estimate properly the parameters of the mixture.

The density is here
proportional to

The true model is
, and
being a parameter that will change, from 0 to 4.
The log likelihood (actually, I add a minus since most of the optimization functions actually minimize functions) is
> logvraineg <- function(param, obs) {
+ p <- param[1]
+ m1 <- param[2]
+ sd1 <- param[3]
+ m2 <- param[4]
+ sd2 <- param[5]
+ -sum(log(p * dnorm(x = obs, mean = m1, sd = sd1) + (1 – p) *
+ dnorm(x = obs, mean = m2, sd = sd2)))
+ }
The code to generate my samples is the following,
>X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z = sample(c(1,2,2),size=n,replace=TRUE)
> X2=m+X20
> X = c(X1[Z==1],X2[Z==2])
Then I use two functions to optimize my log likelihood, with identical intial values,
> O1=nlm(f = logvraineg, p = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)), obs = X)
> logvrainegX <- function(param) {logvraineg(param,X)}
> O2=optim( par = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)),
+ fn = logvrainegX)
Actually, since I might have identification problems, I take either
or
, depending whether
or
is the smallest parameter.

On the graph above, the x-axis is the difference between means of the mixture (as on the animated grap above). Then, the red point is the median of estimated parameter I have (here
), and I have included something that can be interpreted as a confidence interval, i.e. where I have been in 90% of my scenarios: theblack vertical segments. Obviously, when the sample is not enough heterogeneous (i.e.
and
rather different), I cannot estimate properly my parameters, I might even have a probability that exceed 1 (I did not add any constraint). The blue plain horizontal line is the true value of the parameter, while the blue dotted horizontal line is the initial value of the parameter in the optimization algorithm (I started assuming that the mixture probability was around 0.2).
The graph below is based on the second optimization routine (with identical starting values, and of course on the same generated samples),

(just to be honest, in many cases, it did not converge, so the loop stopped, and I had to run it again… so finally, my study is based on a bit less than 500 samples (times 15 since I considered several values for the mean of my second underlying distribution), with 200 generated observations from a mixture).
The graph below compares the two (empty circles are the first algorithm, while plain circles the second one),

On average, it is not so bad…. but the probability to be far away from the tru value is not small at all… except when the difference between the two means exceeds 3…
If I change starting values for the optimization algorithm (previously, I assumed that the mixture probability was 1/5, here I start from 1/2), we have the following graph

which look like the previous one, except for small differences between the two underlying distributions (just as if initial values had not impact on the optimization, but it might come from the fact that the surface is nice, and we are not trapped in regions of local minimum).
Thus, I am far from being an expert in optimization routines in R (see here for further information), but so far, it looks like R is not doing so bad… and the two algorithm perform similarly (maybe the first one being a bit closer to the trueparameter).










J’avais raconté l’autre jour l’aventure de mon fils qui souhaitait faire des triangles avec des bouts de bois (
,
, 


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
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






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 







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
,
et
.
. 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,
.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
, en faisant

. 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
(en noir), en
(en rouge) et en
(en bleu),
Suite rapide de mon billet sur la construction de triangles avec des bouts de bois (
i.i.d. de loi
, et de densité
, et que l’on définie les statistiques d’ordre
, avec
,est alors
(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
, i.e.



. Bref, ça marche pour deux bâtons.


, compte tenu de ce que l’on a écrit auparavant, la loi du couple des deux statistiques d’ordre seraient

.
, avec le sommet en
. Et par des raisons de symmétrie, la probabilité de dépasser
.Merci Adrian. C’est effectivement plus élégant que mon calcul moche…
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
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
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
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
. Mais comme ça n’a que peu d’intérêt, je trouve, on va regarder ce qui se passe à
et que
. En fait, histoire d’épater les filles, Donald prétend que on peut obtenir facilement les résultats suivants

, 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
chances d’atteindre un point bleu, et 


chances d’atteindre un point violet,
chances d’atteindre un point bleu, et enfin
chances d’atteindre un point rouge (il suffit d’étudier toutes les trajectoires possibles). Au total, on arrive à

et que
. 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
la même quantité, avec m le nombre de portions dans le distributeur et
, avec les conditions de bord de la forme
.

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
, alors

(comme on l’a montré) et

. 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,
, on obtient
,
: 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…




Petit résumé du débat du matin avec mon fils …
,
et
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

suivent une loi uniforme sur [0,1].





se positionne par rapport à 1. Deux cas assez différents sont à distinguer: comme on veut
, soit
, soit
. Autrement dit on écrira l’ensemble sur lequel on calcule la première intégrale comme la réunion de deux ensembles (disjoints),


Bref, l’intégrale se retrouve sous la forme de deux intégrales


est petit, ou plutôt
. La seconde intégrale est un cas où
est grand, ou plutôt
. Bref, la première intégrale peut s’écrire




le nombre de carte total (ici 200 car je suppose qu’il existe 200 personnages). Je note
le coût d’une carte cachée (50 centimes ici) et
le coût d’une carte non-cachée (1 euro). Supposons que j’ai à en ma possession les cartes
(et que je n’ai pas les suivantes). Si je veux récupérer une nouvelle carte (que je noterais
), en achetant au hasard i.e. dans des sachets hermétiques, le nombre de cartes à acheter est
, qui suit une loi géométrique, de paramètre
, où
cartes, et qu’il m’en manque simplement une, alors la probabilité de la tirer est
), le coût moyen de ma stratégie est









is said to be an admissible solution for starting point x,
and
. If we reformulate the dynamic programming idea, we obtain that if
, then for all
, sequence
is a solution to problem
. It comes that function v is a solution of Bellman’s equation



,
and
.












. The dynamic is then















. Assuming that

has a form that can be explicited).




that maximizes the following function
(a state space) and
(i.e. some dynamic constraints). Assume further that the starting point is given, i.e.
. In economic application, note that frequently
i.e. we consider a discounted version of the value.




is a solution of the problem
then, for all
,
is a solution of problem
.




is solution of
if and only if for all
, 
denote the consumption at period 
. Let
be capital he got at time t. Assume that his initial capital is a given amount
, and suppose that this period’s capital and consumption determine next period’s capital as
is a positive constant and
. Assume further that capital cannot be negative. Then the consumer’s problem is simply
. Bellman’s equation is then
and
. And to solve that problem, we use backwards induction techniques.
is known, we can derive easily
, and so on until
. More precisely, given
, we can get
which is the maximum of function
. One can see that the following function is a possible solution
is a constant. Further, the optimal amount to consume at time
is



was the variable of interest,
was the control variable, and
. Thus, the programm




Pourtant l’explorateur n’a pas franchement donné d’information… et pourtant cette information a suffi à décimer l’île ! Étonnant non ?