algorithm [@robert2013monte]. Note that to draw $  (\gamma_1, \ldots,\gamma_d)
\sim \mathcal{D}ir(\alpha,\dots,\alpha) $ is equivalent to draw $ \gamma_i = Y_i/(Y_1+\ldots+Y_d) $ with $ Y_i \overset{i.i.d}{\sim} Gamma(\alpha,1), \forall i = 1,\ldots,d $. Thus, instead of $\gamma_i's $, we conduct a MH updating for $ Y_i's $. So the objective is to produce a Markov chain $(Y^{(t)}_1,\ldots,Y^{(t)}_d, V^{(t)}_1,\ldots, V^{(t)}_d) $. From this, we deduce obviously the sequence $(\gamma^{(t)}_1,\ldots,\gamma^{(t)}_d, V^{(t)}_1,\ldots, V^{(t)}_d) $ and use the following empirical mean as the Monte-Carlo approximation of our estimator: $$\hat{\rho}^{{\rm MH}} := \frac{1}{T}\sum_{t=1}^T\left(  \sum_{i=1}^d \gamma_i^{(t)} V_{i}^{(t)} (V_{i}^{(t)})^\dagger\right).$$

For $ t$ from $ 1$ to $ T$, we iteratively update through the following steps:

updating for $Y_i's $:

:   for $i $ from $1 $ to $d $,\
    Sample $\tilde{Y}_i  \sim h(y|Y^{(t-1)}_i) $ where $h$ is a proposal distribution given explicitely below.\
    Calculate $\tilde{\gamma_i} = \tilde{Y}_i / (\sum_{i=1}^d\tilde{Y}_i) $.\
    Set $$Y^{(t)}_i  =
    \begin{cases}
    \tilde{Y}_i    &\text{with probability }\min \left\{ 1,R(\tilde{Y}, Y^{(t-1)})\right\},
    \\
     Y^{(t-1)}_i        & \text{otherwise} 
    \end{cases}$$ where $R(\tilde{Y}, Y^{(t-1)})$ is the acceptance ratio given below.\
    Put $\gamma_i^{(t)} =Y^{(t)}_i /(\sum_{j=1}^dY^{(t)}_j) , i = 1,\ldots,d$.

updating for $V_i's $:

:   for $i $ from $1 $ to $d $,\
    Sample $ \tilde{V}_i $ from the uniform distribution on the unit sphere.\
    