5
$\begingroup$

I want to use the sample average $(X_1 + .... X_n)/n$ as a substitute for the expectation $\mathbb{E}(X)$. As claimed by the weak law of large numbers, as n increases the sample average should converge to $\mathbb{E}(X)$.

I wish to use this logic in my project, where the $X_i$ are iid exponential random variables. A simple code, however, does not demonstrate this well. This is because when we add the $n$ numbers, $X_1, ..., X_n$, it usually results in a large number and most of the precision is lost. So when I divide by $n$, the difference $|(X_1 + ... X_n)/n - \mathbb{E}(X)|$ is never very small even when I increase $n$ by a large amount.

I have tried some simple manipulations such as taking sub sums and then taking the total average. Even here I seem to be suffering from the same problem.

Is there a neat way of achieving errors less than, say, $10^{-8}$.

$\endgroup$

migrated from codereview.stackexchange.com Mar 13 '12 at 15:02

This question came from our site for peer programmer code reviews.

  • $\begingroup$ a) When you say you "add the n numbers, X_1, ..., X_n" are you referring to summing n realizations from the X_i distribution, or are you forming the sum of the RVs and drawing realizations from that? b) I don't recall the formulation, but I believe the sum of n iid variables that follow an exponential distribution have a known distribution, so E(X) is known analytically. Unless I'm wrong (always a real possibility), why estimate it? $\endgroup$ – Barron Mar 13 '12 at 17:47
  • 1
    $\begingroup$ For distributions such as the exponential distribution, where the first moment is well-defined and the law of large numbers applies, the expected relative error in $\frac{\sum_i X_i}{n} - EX$ is on the order of magnitude of $1/\sqrt{n}$, so if you want a relative error in the order of $10^{-8}$, you'll need about $10^{16}$ realizations as a ballpark figure. Are you sure you want to invest that much effort into it? $\endgroup$ – Erik P. Mar 13 '12 at 20:13
5
$\begingroup$

Knuth (TAOCP Volume 2, 3rd ed., pg. 232) suggests using the formula $M_k = M_{k-1} + (x_k - M_{k-1})/k$ to calculate the mean, where $x$ is a vector of your samples.

See also: a stackoverflow question, an article on accurately computing running variance, and a paper on computing covariances and arbitrary-order statistical moments

$\endgroup$
2
$\begingroup$

Since you are summing up exponentials you may want to look at the answer provided to this SC question.

$\endgroup$
2
$\begingroup$

If you would compute in exact arithmetic, there is no way of rewriting your formula to get better accuracy; you cannot do better than $O(\sigma/\sqrt{N})$.

But if you get low accuracy just because of roundoff errors, you can use the formula $E(X)=x+ N^{-1}\sum (X_k-x)$, using for $x$ an arbitrary approximation of $E(X)$, e.g., the mean of the first few terms.

You can reduced the rounding errors further by summing the entries of this sum in the ordering by their absolute values, starting with the smallest.

$\endgroup$
1
$\begingroup$

I think you can get around the precision problem by computing it as ((X_1-E(X))+...+(X_n-E(X)))/n.

$\endgroup$
  • $\begingroup$ I'm not sure what you are saying here, since $E(X)$ is unknown until it is computed, unless you are using $E(X)$ here to refer to a running mean, in which case your answer would be very similar to JamesCuster's. $\endgroup$ – Geoff Oxberry Mar 14 '12 at 0:49
  • $\begingroup$ @Geoff: Yes, I assumed that E(X) is known, and that you just want to estimate the difference of the average and the expected value. In the case with the exponential distribution the expected value is surely known. In case it's not known, my method is of course not applicable, and then James's solution is probably the best. $\endgroup$ – Pontus von Brömssen Mar 14 '12 at 9:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy