# How can I calculate variance of a very large random variable?

I'm implementing an algorithm which recieves as input samples from a random variable with an unknown distribution.

The random variable is extremely large so my input is logarithmic, and still large (in the tens of thousands).

I want to calculate its mean and variance, or even log[mean] & log[var], but am having trouble because I can't use an exponent on the input.

When faced with a similiar problem, but with very small numbers, to perserve accuracy I used the following trick to calculate log[mean]:

def log_expectation(log_samples):

a_max = max(log_samples)
b = log_samples - a_max
n = len(log_samples)

log_mean = a_max + np.log(sum(np.exp(b))) - np.log(n)

return log_mean


Can you suggest a similar trick for large numbers? What else can I do?

• What sort of data are you dealing with? There are no numbers in $e^{1000}$ range that we deal with in business or even science Commented Jul 9, 2016 at 0:28
• @Aksakal, I'm trying to calculate $E[\frac{1}{P(X|G)}]$, where $X$ is The aligned full-sequence dna of individuals from different populations (of the same species), and $G$ are genealogies over independent loci across the whole genome. I have tens-of-thousands of genealogies, and the Data-Likelihood in each of them is small. Thats why the original expression ends up being very large. What will eventually interest me is the ratio between two such calculation, so hopefully my results will remain sane. Thanks! Commented Jul 9, 2016 at 21:10

Not sure how this really helps but you can use the same trick to get something proportional to the mean, with known proportionality constant.

E.g. say the log values are $1000$ and $1001$. Subtract the max to get the log of the scaled values: $-1, 0$. So, the scaled average is $(1 + e^{-1})/2$. Scaling up by the proportionality constant, the average is $e^{1001}(1 + e^{-1})/2$.

Equivalently, the log of the mean is $1001 + \log(1 + e^{-1}) - \log(2)$

• And how would you calculate Variance? Commented Jul 10, 2016 at 0:29
• Not sure. I assume you could factor out a large term and use a similar trick to calculate the log expectation of $(X-E(X)) ^2$. Don't have time to think through it right now. If you do figure it out, please post something here- I'd be interested to see it. Commented Jul 10, 2016 at 2:41

I crunched it out and am satisfied with my results.

Given the large samples of $ln(X)$, I calculate $ln(Var(X))$ using the following expansion-

$ln(Var(X)) = ln(Var(\frac{X}{C})*C^2) = ln(Var(\frac{X}{C})) + 2ln(C) =$ $ln(Var(e^{ln(X)-ln(C)})) + 2ln(C)$

When I set $ln(C):= Max(ln(X))$, the calculation is stable, since most of the calculations are in log-scale, and the exponent is applied on small, rather than large numbers.

My Code:

def log_variance(log_samples):

ln_C = max(log_samples)

a = log_samples - ln_C
b = np.exp(a)
c = b.var()
d = np.log(c)
e = d + 2*mx

return e