← Domaines d'expertise/Choix discret · Estimation/2025

BHAMSLE : +17,2 % de gain de log-vraisemblance sur la mixture mode-choice de Londres, en moitié moins de temps que Multistart Biogeme

Un Breakpoint Heuristic Algorithm pour la Maximum Simulated Likelihood Estimation des modèles de choix discret avancés — exploite la même structure combinatoire qui a propulsé nos algorithmes de tarification OR Spectrum

L'estimation des modèles de choix modernes à classes latentes et à mélanges discrets-continus est un cauchemar de calcul silencieux. Des centaines de maxima locaux, du bruit de simulation, des surfaces de vraisemblance qui ressemblent aux Alpes. Les solveurs standards trouvent un sommet, généralement le mauvais.

Nous avons pris la structure de points de rupture découverte pour la tarification basée sur les choix et l'avons portée vers le problème d'estimation. Le résultat, BHAMSLE, trouve des log-vraisemblances plus élevées, des segments latents plus interprétables et une meilleure récupération des vrais paramètres que CMA-ES, l'estimateur MILP existant et même Biogeme multistart avec 50 redémarrages aléatoires.

+17,2 %

gain de log-vraisemblance sur la mixture discrète-continue London mode-choice

Moitié du temps

de Multistart Biogeme M=20, le nombre minimum de redémarrages qui atteint la solution BHAMSLE

7,7×

plus rapide que Multistart Biogeme M=100 à log-vraisemblance finale identique

Retrouve la vérité

split de classe latente (0,69, 0,31) là où Biogeme par défaut reste bloqué à uniforme (0,50, 0,50)

Le défi

L'estimation par maximum de vraisemblance simulée des DCMs avancés, modèles à classes latentes avec mélanges continus, ensembles de choix restreints, et spécifications flexibles d'utilité, est plombée par la multimodalité et le bruit induit par la simulation. La surface de vraisemblance d'un simple logit à deux classes peut montrer plusieurs bassins distincts à haute vraisemblance ; les estimateurs de type Newton convergent vers celui qui se trouve le plus proche du point de départ.

La seule approche structurée dans la littérature est la reformulation MILP de Fernández Antolín (2018), qui garantit l'optimalité globale en principe mais ne passe pas à l'échelle au-delà de minuscules instances. Les contournements standards des praticiens, multistart depuis des initialisations aléatoires, CMA-ES, PSO, paient la multimodalité par brute force et convergent régulièrement vers de l'insensé interprétable (une classe avec 100 % de probabilité, ‘distributions' à variance nulle, paramètres de mélange effondrés).

Notre démarche

Nous avons observé que la MSLE basée sur la simulation partage le même squelette combinatoire que la tarification basée sur les choix : des variables indicateurs par individu ω_inr qui basculent à un nombre fini de points de rupture dépendant des paramètres, avec l'objectif (log-vraisemblance ici, revenu là) constant par morceaux entre les bascules.

Sur cette observation nous avons construit BHAMSLE : un algorithme d'ascension par coordonnées qui parcourt l'espace des paramètres point de rupture par point de rupture, évaluant la log-vraisemblance simulée exactement à chaque bascule, exploitant une structure dont CMA-ES, PSO et les méthodes de gradient sont aveugles. Deux sources de points de rupture sont traitées séparément : points de rupture des paramètres d'utilité (où ω bascule quand β change, à assignation de classe fixée) et points de rupture de fonction de score (où l'appartenance à la classe elle-même bascule quand α ou γ change).

BHAMSLE n'est pas positionné comme un estimateur final mais comme une stratégie d'initialisation : placez-le en amont de la routine Newton standard de Biogeme et regardez la solution finale s'améliorer, la recherche exploitant la structure coûtant des minutes au lieu des heures du multistart Biogeme.

Le résultat

Livrable : BHAMSLE, un schéma d'initialisation par heuristique en points de rupture pour la Maximum Simulated Likelihood Estimation des modèles de choix discret avancés (classes latentes, mixtures continues, mixtures discrètes-continues). La méthode exploite la même structure combinatoire qui a propulsé nos algorithmes de tarification OR Spectrum, appliquée ici à la log-vraisemblance au lieu du revenu.

Sur la mixture discrète-continue London mode-choice (N = 500, R = 500), Biogeme initialisé par BHAMSLE atteint LL = −243,68 vs −294,38 avec l'initialisation par défaut, une amélioration de +17,2 %, en 3,8 heures. Multistart Biogeme a besoin de 20 redémarrages (7,4 heures, environ 2× le temps BHAMSLE) pour égaler ; avec 100 redémarrages, il prend 7,7× le temps BHAMSLE pour la même qualité finale. Sur le test de classes latentes discrètes, Biogeme initialisé par BHAMSLE retrouve le vrai split de classes (0,69, 0,31) là où Biogeme par défaut reste bloqué à uniforme (0,50, 0,50).

Plongée technique

Le modèle derrière le résultat.

Le modèle

L'estimation par maximum de vraisemblance simulée d'un DCM avancé a une structure fondamentalement différente du logit de manuel. Les probabilités de choix n'ont pas de forme close (les modèles de mélange intègrent sur des distributions continues de paramètres de goût), donc nous les approchons par Monte Carlo. Cela introduit R scénarios simulés par individu, et à chaque valeur de paramètre, chaque paire (n, r) a un indicateur déterministe ‘l'alternative observée gagne' qui bascule on/off selon que le paramètre est au-dessus ou en dessous d'un seuil spécifique au client et au scénario. La structure de points de rupture que BHAMSLE exploite est exactement l'ensemble des valeurs de paramètre auxquelles l'une de ces bascules se produit.

β, π max n \in N \sum ln (s \in S \sum π_{n s} P_{y_{n} n}^{s} (β)), P_{in}^{s} = P (U_{in}^{s} \geq U_{j n}^{s} \forall j \in C_{n}^{s}) .

Le problème MSLE sous forme stochastique. Somme sur les individus du log d'une moyenne pondérée par probabilité des probabilités de choix par classe. Les variables de décision sont les coefficients d'utilité β (souvent dépendant de la classe) et les probabilités d'appartenance de classe π. Pour les mélanges continus, P^s_in est une intégrale sur une densité normale de paramètres de goût distribués, d'où la nécessité de la simulation.

U_{in r}^{s} = k \in K_{s} \sum x_{ink} β_{k} + ε_{in r}, ω_{in r}^{s} = 1 [U_{in r}^{s} = j \in C_{n}^{s} max U_{j n r}^{s}] .

Cadre de simulation. Pour chaque scénario r et individu n, l'utilité de l'alternative i sous la classe s devient déterministe une fois le terme d'erreur ε_inr tiré. L'indicateur de choix ω^s_inr est 1 ssi i gagne l'arg-max pour (n, s, r). La moyenne de ω sur les R tirages donne un estimateur sans biais de P^s_in (Train 2003).

s LL (β, π) = - N ln R + n \in N \sum ln (Σ_{n} (β, π)), Σ_{n} := r \in R \sum s \sum 1_{class assignment} \cdot ω_{y_{n} n r}^{s} .

Log-vraisemblance simulée après Monte Carlo. Σ_n compte combien des R scénarios choisissent l'alternative observée y_n pour l'individu n (sous l'assignation de classe du scénario). L'objectif est la somme des log-Σ_n sur les individus, une fonction constante par morceaux en β et π, sautant à chaque point de rupture où un ω bascule.

[s_{1}, s_{2}] : U_{y_{n} n r}^{s} (β_{j}) \geq U_{in r}^{s} (β_{j}) \forall i \in C_{n}^{s} . (entry breakpoint s_{1}, exit breakpoint s_{2})

Construction de point de rupture pour un paramètre d'utilité β_j. En figeant tous les autres paramètres, le segment de β_j pour lequel l'alternative observée y_n reste dominante pour l'individu n en scénario r sous la classe s est borné par deux points d'indifférence. La borne inférieure s_1 est un point de rupture ‘entrée' (ω bascule on quand β_j croît au-dessus) ; la borne supérieure s_2 est une ‘sortie' (ω bascule off). Sur tous les (n, r, s) on obtient au plus 2NRS points de rupture par paramètre.

At an entry breakpoint: Δ s LL = ln (Σ_{n} + 1) - ln (Σ_{n}); at an exit: Δ s LL = ln (Σ_{n} - 1) - ln (Σ_{n}) .

Mises à jour incrémentales de log-vraisemblance. En parcourant les points de rupture dans l'ordre trié sur la j-ième coordonnée, le changement de log-vraisemblance simulée à chacun est juste la différence de deux log-comptages, travail O(1) par point de rupture, vs. le coût O(NR) qu'un optimiseur boîte-noire paierait pour ré-évaluer l'objectif complet. C'est ce qui permet à BHAMSLE de monter en charge : un passage complet de coordonnée coûte O(NRS) au lieu de O((NR)²) pour une recherche en ligne naïve.

β_{j}^{(t + 1)} = ar g β_{j} \in B_{j} (β_{- j}^{(t)}) max s LL (β_{- j}^{(t)}, β_{j}), j (t + 1) = (j (t) mod (K + S)) + 1.

Boucle externe BHAMSLE. Ascension par coordonnées sur les K paramètres d'utilité et les S − 1 paramètres de probabilité de classe. À chaque pas nous figeons K + S − 2 paramètres, énumérons les points de rupture pour la coordonnée active, et choisissons le point qui maximise la log-vraisemblance simulée. Terminer quand un passage complet de K + S coordonnées ne produit aucune amélioration au-dessus de τ = 10⁻⁶.

T_{BHAMSLE} = O (C \cdot (K + S) \cdot N R S) ≪ T_{Multistart Biogeme} = O (M \cdot T_{Newton}) .

Comparaison de complexité. C est le nombre de cycles externes d'ascension par coordonnées (empiriquement < 20 sur les benchmarks) ; le terme interne NRS est l'énumération des points de rupture et la mise à jour incrémentale. Multistart Biogeme doit payer une convergence Newton complète par redémarrage M, et M doit croître avec le nombre d'optima locaux, qui pour les mélanges discrets-continus se compte en centaines. Sur le mélange Londres du Test 4, BHAMSLE atteint multistart Biogeme à M = 20 en environ la moitié du temps wall-clock et surpasse M = 10.

Ce que nous avons fait ici : prendre une observation structurelle d'une classe de problème totalement différente, la tarification, et remarquer qu'elle se transfère à l'estimation, car les deux problèmes calculent un objectif externe sur des indicateurs de choix simulés. La même énumération de points de rupture qui pilotait BEA et BHA dans les articles de tarification pilote ici BHAMSLE, avec deux nouvelles subtilités : (1) l'objectif a un log au lieu d'une somme, donc nous utilisons des mises à jour incrémentales en ln au lieu d'additions, et (2) les paramètres de classes latentes vivent à l'intérieur d'indicateurs pour l'assignation de classe, pas pour le choix d'alternative, donc ils ont besoin de leur propre logique de points de rupture. Le gain est cohérent sur chaque test de l'article : Biogeme initialisé par BHAMSLE trouve des log-vraisemblances plus élevées que l'initialisation par défaut, les trouve plus vite que CMA-ES ou PSO ne le peuvent, et les trouve avec moins de budget de calcul que le contournement multistart de manuel.

Benchmark

Test 4 de l'article : mélange discret-continu de logit sur le dataset London mode-choice (4 alternatives : marche, vélo, transport public, voiture ; N = 500, R = 500 ; vraie répartition de classes 0,70/0,30 ; sensibilité au coût normalement distribuée). Chaque méthode est utilisée comme initialisation pour l'estimateur de type Newton de Biogeme ; nous rapportons la log-vraisemblance finale atteinte (plus haute = meilleure, donc moins négative) et le temps wall-clock total. MSB = Multistart Biogeme avec M redémarrages aléatoires dans [-12, 12].

Méthode d'initialisation	Log-vraisemblance finale	ΔLL vs default	Temps total (s)	ΔTemps vs default
Biogeme (default)	−294,38	-	583	-
Biogeme + BHAMSLE (notre)	−243,68	+50,69 (+17,2 %)	13 796	+13 212 s (+2 265 %)
MSB (M = 10)	−243,78	+50,60 (+17,2 %)	15 003	+14 419 s (+2 472 %)
MSB (M = 20)	−243,69	+50,69 (+17,2 %)	26 706	+26 122 s (+4 478 %)
MSB (M = 50)	−243,68	+50,69 (+17,2 %)	68 062	+67 479 s (+11 568 %)
MSB (M = 100)	−243,68	+50,69 (+17,2 %)	106 438	+105 854 s (+18 151 %)

Deux choses se détachent. Premièrement : l'initialisation Biogeme par défaut reste coincée dans un optimum local 17 % pire que le global, et les méthodes Newton ne peuvent pas le détecter. La recherche structurée de BHAMSLE corrige cela en routant l'estimateur au-delà du mauvais bassin avant qu'il ne démarre. Deuxièmement : BHAMSLE retrouve la solution finale que multistart trouve avec M = 20 redémarrages, en environ moitié du temps wall-clock. Avec M = 10 redémarrages, multistart est même légèrement en deçà de BHAMSLE en qualité et le dépasse en temps. L'implication : les redémarrages brute force paient une lourde taxe parce qu'ils ne savent pas où chercher, alors que BHAMSLE sait où chercher car la structure du problème le lui dit. Le même schéma se répète sur les quatre montages expérimentaux de l'article : mélanges discrets avec choix observés et synthétiques, mélanges discrets-continus avec les deux, sur les données Swiss Metro et London mode-choice.

Tiré du dossier

Section 2.1 de l'article : la mise en place du modèle de choix à classes latentes. Chaque individu n est supposé appartenir à l'une des S classes latentes ; chaque classe a son propre vecteur de paramètres β^(s) et son ensemble de choix C_n^s. Les probabilités d'appartenance π_ns sont elles-mêmes estimées sur les données, soit directement, soit via des fonctions de score f_ns. La surface de vraisemblance, même de ce simple exemple à deux classes et une variable explicative, présente plusieurs bassins distincts à haute vraisemblance, motivation de tout ce qui suit.

Formulations (1)–(6) de l'article : le problème MSLE comme maximisation sous contraintes sur les coefficients d'utilité β, les coefficients de scoring de classe α, les probabilités de choix P^s et les probabilités d'appartenance π. Le ℙ(U_in^s ≥ U_jn^s) sans forme close dans la troisième contrainte est exactement où la simulation entre, pour les DCMs avancés cette intégrale n'a pas de solution close, donc nous la remplaçons par une moyenne d'échantillon Monte Carlo et le problème devient un MILP.

Étapes 1–4 de l'algorithme BHAMSLE : partant d'un vecteur de paramètres initial, fixer toutes les coordonnées sauf β_j, puis itérer sur toutes les paires (individu, scénario) et calculer l'intervalle [s_1, s_2] des valeurs β_j pour lesquelles l'alternative observée reste dominante en utilité. Les extrémités de cet intervalle sont un point de rupture d'entrée et un de sortie pour la vraisemblance simulée. L'ensemble des points de rupture sur (n, r, s) donne les mouvements candidats pour cette coordonnée.

Étape 4 de l'algorithme suite : quand la coordonnée active est un paramètre de probabilité de classe γ_g (plutôt qu'un coefficient d'utilité), la logique de points de rupture change. Ici l'algorithme demande quelle assignation de classe gagne actuellement pour chaque (n, r), puis calcule le segment des valeurs γ_g pour lequel cette assignation continue de gagner, et enregistre les points de rupture d'entrée/sortie en conséquence. Cette construction double de points de rupture (utilité vs appartenance de classe) est l'une des deux généralisations non triviales que BHAMSLE apporte sur le BHA originel des articles de tarification.

Techniques

Énumération de points de rupture pour vraisemblance simulée
Ascension par coordonnées sur paramètres MSLE
Estimation de modèles à classes latentes / mélange
Mises à jour incrémentales ln(Σ ± 1) − ln(Σ)
Comparaison vs CMA-ES, multistart Biogeme

Stack

Python
Biogeme 3.2.14
Julia (CMA-ES via Evolutionary.jl)
Gurobi (baseline MILP)

Un problème de ce type ?

Nous serions ravis d'en entendre parler.

contact@optimally.ch