Numerical Precision Weak Law of Large numbers

I want to use the sample average (X_1 + .... X_n)/n as a substitute for the expectation, i.e. E(X). As claimed by the weak law of large numbers, as n increases the sample average should converge to 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 - 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.

Answer 1

[Jan]

You do not loose precision as long as you do not cast the result to single. You loose accuracy, when you sum elements of different magnitudes, e.g. 1e17 + 1e0 - 1e17 is 0 when calculated using doubles. You can use a compensated sum, which stores the intermediate values using a higher precision, e.g. FEX: XSum.

[Lal]

I will see if XSum helps in anyway

Answer 2

[Peter Perkins]

Slow convergence is the reason why variance reduction methods such as importance sampling, antithetic sampling or quasi-random sequences are sometimes used.

[Lal]

Thanks. I will look into these.

Answer 3

[Tom Lane]

While I don't know exactly what you are encountering, I do know it is sometimes possible to repeat the mean calculation on residuals to recover more precision. Example:

>> x = exprnd(1000,1e6,1);
>> m = mean(x)
m =
     9.990841964835327e+02
>> mean(x-m)
ans =
     7.003778591752052e-12
>> m = m + mean(x-m)
m =
     9.990841964835397e+02
>> mean(x-m)
ans =
    -3.860623110085726e-14

[Lal]

I want the error between 1000 and m to be small. Your idea suggests one can improvise on estimate of m

[Tom Lane]

When you wrote "precision is lost" I took that to mean you are concerned about the calculation itself. If you are concerned about the difference between the sample mean m and the theoretical mean 1000, that's different. The standard deviation of an exponential random variable with mean M is sqrt(M). The standard deviation of the average of a sample of size N is therefore sqrt(M/N). That means to gain an extra digit of accuracy, you need to increase N by a factor of 100.