Simulation avec contraintes

dedekind93

Bonjour à tous.
Je vous souhaite avant tout la meilleure santé qui soit à vous et vos proches.

Comment créeriez-vous algorithmiquement parlant $n+1$ réels $a_0,a_1,\dots,a_n\in[0,1]$ de somme $1$ et tels que $\displaystyle\sum_{k=0}^nka_k=m$ où $m$ est fixé à l'avance ?
Le tout le plus uniformément possible, même si ceci n'a peut-être pas grand sens.
Merci de vos lumières !

lourrran

$m$ est fixé à l'avance. Tu le dis noir sur blanc.
$n$ est-il fixé à l'avance ?

purple

Ça a tout l'air d'une loi de Dirichlet de paramètre $(1, 2, \ldots, n)$ à un coefficient multiplicatif près.
On peut la simuler avec des distributions Gamma je crois.
Ton $a_0$ ne sert à rien en l'état :-D.

YvesM

Bonjour,

On peut faire un calcul exact si l'on pose $a_k = a k+b$... $a = {6 (2 m - n-1)\over n (n^2-1)}$ et $b = {-6m+4n+2 \over n(n-1)}$ pour $n \geq 2.$

dedekind93

Merci pour vos réponses. J'ai oublié de préciser que $n$ était fixé, et que je souhaitais quelque chose d'aléatoire, mille excuses.

purple

Dans ce cas je pense que c'est bien une loi de Dirichlet que tu cherches à simuler, et que tu pourras construire avec des variables de loi Gamma indépendantes. Comme les paramètres sont entiers tu devrais pouvoir programmer ça en utilisant des variables exponentielles, qui se simulent très facilement à partir d'un générateur uniforme.

Saturne

Pour tout $\alpha > 0$ et tout $\beta > 0$, on a $$\sum_{k=0}^n {n \choose k}\frac{B(k+\alpha, n-k+\beta)}{B(\alpha,\beta)} = 1$$ et $$\sum_{k=0}^n k{n \choose k}\frac{B(k+\alpha, n-k+\beta)}{B(\alpha,\beta)} = \frac{n\alpha}{\alpha+\beta},$$ où $B$ est la fonction bêta.

Il me semble que la valeur maximale possible de $m$ est $m = n$ pour que ton problème ait une solution (correspondant alors à $a_0 = \cdots = a_{n-1} = 0$ et $a_n = 1$), et $m=0$ est la valeur minimale possible ($a_0=1$).

Supposons alors $0 < m < n$. Voici une méthode :

- prendre $\alpha > 0$ aléatoire, avec la loi que tu veux ;

- prendre $\beta = \frac{n-m}{m} \alpha$ (de sorte que $\frac{n\alpha}{\alpha+\beta} = m$) ;

- prendre $a_k = {n \choose k}\frac{B(k+\alpha, n-k+\beta)}{B(\alpha,\beta)}$.

Pour l'aspect "le plus uniformément possible", je ne sais pas. Ça va dépendre de la loi de $\alpha$.

dedekind93

Merci beaucoup c'est très clair...et finalement pas si évident !

Saturne

En fait, j'ai pris $a_k = \Pr(X=k)$ où $X$ suit la loi bêta-binomiale, ainsi $\sum a_k = 1$ et $\sum k a_k = \mathbb{E}(X)$.

Si tu fais le même truc avec une loi sur $\{0, \ldots, n\}$ qui a plus que deux paramètres, tu auras plus de possibilités pour la distribution de $(a_0, \ldots, a_n)$. Il y a une loi bêta-binomiale généralisée avec 3 paramètres, mais ses masses de probabilité s'expriment avec la fonction hypergéométrique de Gauss.

Saturne

J'ai trouvé une librairie R, hitandrun, qui permet de simuler la loi uniforme sur un polytope convexe défini par des contraintes linéaires. Voilà comment réaliser tes simulations avec cette librairie :

library(hitandrun)

n <- 2
m <- 1.5
constraint1 <- simplexConstraints(n+1)
constraint2 <- list(constr = t(0:n), rhs = m, dir = "=")
constraints <- mergeConstraints(constraint1, constraint2)
sims <- hitandrun(constraints, n.samples = 10000)

Saturne

Et polyapost, une librairie qui permet de simuler une loi de Dirichlet contrainte à des égalités ou inégalités linéaires.