Vtyjkyuk

Given a discrete probability distribution, $boldsymbollambda = (lambda_1,lambda_2, ..., lambda_n)$, how can we construct a random distribution,
$boldsymbolLambda = (Lambda_1, Lambda_2, ..., Lambda_n)$ such that $mathbbE(boldsymbolLambda) = boldsymbollambda$?

The goal is to generate random variations of the initial distribution, $boldsymbollambda$, with the above property. If it is possible to do this and have some degree of control over the variance/spread of each $lambda_i$, then even better!

I suspect that there are many ways to achieve this, so any approaches/answers are appreciated.

What I've tried:

My initial thought was to let each $Lambda_i sim Beta(alpha_i, beta_i)$, with $alpha_i$, $beta_i$ chosen so that $mathbbE(Lambda_i) = lambda_i$. By choosing larger $alpha_i$, and $beta_i$ values, we can also make the distributions tighter around the mean.

However, this has the issue that, in general,

$$
sum_i=1^n Lambda_i ne 1
$$

We may work around this by normalizing $boldsymbolLambda$ after evaluating each $Lambda_i$ (by dividing through by sum). Unfortunately I think I am right in saying that, in general, this means that $mathbbE(boldsymbolLambda) = boldsymbollambda$ no longer holds.

Similarly, we could 'normalize' $boldsymbolLambda$ by setting

$$
Lambda_n = 1 - sum_i=1^n-1 Lambda_i
$$

This preserves $mathbbE(boldsymbolLambda) = boldsymbollambda$ (at the expense of the distribution of $Lambda_n$ now being difficult to determine exactly, though we can approximate it). The problem with this method is that $Lambda_n$ could now be negative!

You could generate a random vector that sums to $1$ and subtract $frac1n$ from each component, then add some multiple of the resulting zero-sum vector to $boldsymbollambda$. A straightforward way to generate a random vector that sums to $1$ while keeping all components on an equal footing is to use the $n$ intervals formed by $n-1$ points independently uniformly chosen in the unit interval. Then $mathbbE(boldsymbolLambda) = boldsymbollambda$ follows by symmetry. You could ensure that the $Lambda_i$ lie in $[0,1]$ by adding a sufficiently small multiple of the random vector. (The factor would need to be determined beforehand, independent of the random vector, in order to preserve $mathbbE(boldsymbolLambda) = boldsymbollambda$.)

Rather embarrassingly, I seem to have overlooked the Dirichlet distribution!

Nevertheless, any other approaches, would be much appreciated! In particular, methods of generating $boldsymbolLambda$ that give easier control over the spread/variance of each $Lambda_i$ would be useful.