How do I propagate correlated errors numerically?

Question

I'm facing an error propagation problem in fitting some experimental data. I have measured several quantities, $m_i$, and I know from theory that $\sum_{i=0}^{n} m_i = 1$. Each of the $m_i$ has its own measurement error, $\delta_i$ reported as a standard deviation of $m_i$. The covariance $\sigma_{ij}$ of any two $m_i, m_j$ can be calculated from experimental data; we'll assume they follow natural correlation (see comment).

I'm using the $m_i$ to derive another property, $R_i$. The transformation is complicated enough that I can't analytically propagate the errors in $m_i$. Instead, I hope to do so numerically: for $N$ cycles, I want to perturb the $m_i$ within measurement error and recalculate the $R_i$; I will save the calculated $R_i$ values and use these to calculate the errors in my products.

However, I'm not sure how to properly perturb the $m_i$ measurements, as they are correlated with each other. I've considered re-normalization (i.e. perturb all within error and then divide such that $\sum_{i=0}^{n} m_i = 1$). I've also considered using the covariances, in the following manner: suppose $m_1 = 0.4 \pm 0.1 $, $m_2 = 0.3 \pm 0.1$, and $m_3 = 0.3 \pm 0.1$. I first perturb $m_1$ within error; maybe I get $m_1 = 0.45$. I then propagate this perturbation into $m_2$ and $m_3$; assuming equal covariances, I get $m_1 = 0.45, m_2 = 0.275, m_3=0.275$. I then proceed to perturb $m_2$, propagate based on covariance to satisfy the closure condition, then perturb $m_3$, and do the same.

But I don't have any theoretical justification for these ideas (or indeed formal training in statistics). What is a theoretically sound way to do this perturbation (or at least a reference that discusses this?)

Answer 1

Random proportion vectors $\boldsymbol{p}=(p_1, \ldots, p_n)$ are usually negatively correlated. The usual way to sample proportions is to assume an underlying Dirichlet distribution. In your case, you may consider sampling from $Dir(N\boldsymbol{m})$ where $\boldsymbol{m}$ are the observed proportions and $N$ is a constant you choose such that the standard deviation/variance match what you expect. The approach of sampling from $m_i \pm \delta_i$ and normalizing might work, except for the difficulty in choosing an appropriate distribution for sampling. For example, if we assume $N(m_i, \delta_i^2)$, then there is the possibility of obtaining a negative result.