# Calculate variance for a lottery?

I'm trying to calculate the variance for a lottery using the mathematica function "variance". In the following a simplified example with just a few values, and how I tried to solve it:

payouts = {7, 3, 1, 0};
probabilities = {1/6, 1/6, 1/3, 1/3};

Variance[WeightedData[payouts, probabilities]]
Mean[WeightedData[payouts, probabilities]]


The Mean gets calculated correctly (= 2), but the output for variance is 108/13 instead of 6. What am I doing wrong? Should I use another function?

• Looks like a bug to me. Send a report to: [email protected] Commented Jan 4, 2021 at 9:44
• @DanielHuber It's not a bug . Variance works properly. Commented Jan 4, 2021 at 9:46
• I need to stress again, in this case we have the population variance not the samples variance. The given data is the population, not a sample! Therefore, we need not guess the mean, it is given accurately by the population mean. Commented Jan 10, 2021 at 14:40

@ Mariusz Iwaniuk. If you are of the opinion that MMA is right, then please explain what is wrong with the following calculation:

Note that the given values are the whole population, not a sample!

payouts = {7, 3, 1, 0};
probabilities = {1/6, 1/6, 1/3, 1/3};
mean=payouts . probabilities
var= ((payouts - mean)^2) . probabilities
Variance[WeightedData[payouts, probabilities]]
(*6*)
(*108/13*)


Our calculation gives 6 for the value of the variance, but MMA claims: 108/13.

• Is this the good-old chestnut of maximum-likelhihood vs unbiased variance estimate? Commented Jan 4, 2021 at 11:11
• @Sjoerd, you are correct. In particular, Variance[] is using the unbiased estimate featured here: wd = WeightedData[payouts, probabilities]; probabilities . (payouts - Mean[wd])^2/(Total[probabilities] - probabilities . probabilities/Total[probabilities]). Always a good idea to check your definitions before calling things bugs. Commented Jan 4, 2021 at 11:14
• Variance[WeightedData[payouts, probabilities]] == var /(1-Total[probabilities^2) , that is.
– kglr
Commented Jan 4, 2021 at 11:24
• @DanielHuber You're calculating the mean with mean = payouts . probabilities, so it's estimated from the data. Just because the probabilities are known doesn't mean that the mean is known. That's why the correction is needed. Commented Jan 4, 2021 at 11:33
• @DanielHuber In that case you need to use EmpiricalDistribution, not WeightedData. Variance @ EmpiricalDistribution[probabilities -> payouts] gives you the result you seek. Commented Jan 4, 2021 at 11:41

In response to Daniel Huber's answer: this is a case of which variance estimator you want to use. The easiest way to see this, is to make the probabilities all the same, so that the weights don't actually do anything:

payouts = {7, 3, 1, 0};
probabilities = {1/4, 1/4, 1/4, 1/4};
Total[probabilities]


1

These two are now the same (as they should be):

Variance[payouts]
Variance[WeightedData[payouts, probabilities]]


115/12

115/12

The other calculation gives something different:

mean = payouts . probabilities
((payouts - mean)^2) . probabilities


115/16

That's because that formula is for the maximum-likelihood estimate of the variance, i.e. the 1/n estimate instead of the 1/(n - 1) one:

Total[((payouts - mean)^2)]/Length[payouts]


115/16

I'm not going into the details of the pros and cons of each estimator, but it suffices to say that Mathematica uses the unbiased one.

# Edit

It seems there is some confusion about the intended use case of WeightedData. This function is for representing data obtained from a sample, as the name suggests. If you want to calculate the variance of a theoretical distribution where the payouts and weights are exactly known, WeightedData is not the correct way to represent this distribution. You need to use EmpiricalDistribution instead:

payouts = {7, 3, 1, 0};
probabilities = {1/6, 1/6, 1/3, 1/3};
Variance @ EmpiricalDistribution[probabilities -> payouts]


6

• (Perhaps we should make the suggestion that this should be mentioned in the docs.) Commented Jan 4, 2021 at 11:26
• I thoroughly agree with this. Commented Jan 4, 2021 at 11:53
• @DanielHuber The documentation for WeightedData does make a short remark about this in the "Properties & Relations" section, but yeah: it's not super obvious. Commented Jan 4, 2021 at 12:07

You can look at one dataset as consisting of distinct elements occuring with certain probabilities or as containing the same elements as before with multiplicities according to their probability like:

data = {7, 3, 1, 0};
probabilities = {1/6, 1/6, 1/3, 1/3};
data1 = {7, 3, 1, 1, 0, 0};
probabilities1 = {1/6, 1/6, 1/6, 1/6, 1/6, 1/6};


Nevertheless Mathematica's variance of the same data gives different values :

{Variance[WeightedData[data, probabilities]],Variance[WeightedData[data1, probabilities1]]}
(* {108/13, 36/5}  *)


It seems to me that this is not explained by n vs n - 1. Exactly what correction is applied to the variance of weighted data?

Edit:

Ok, that means the estimation of the Variance depends to some extend on the grouping of the data.

• A manual implementation like wvar[dat_, prob_] := prob . (dat - prob . dat)^2/(1 - prob . prob) might clarify things for you. Then, compare e.g. wvar[data1, probabilities1] with Variance[WeightedData[data1, probabilities1]]. Commented Jan 4, 2021 at 14:25