# Expectations of sums

I would like to use Mathematica to derive some bounds on empirical estimators, such as $$E[Y]$$ where $$Y = \tfrac1n\sum_{i=1}^n (X_i - X)^2$$ and $$X = \tfrac1n\sum_{i=1}^n X_i$$.

For a moment I thought this worked:

Expectation[Sum[X[i], {i, 1, n}], {X[i] \[Distributed] BernoulliDistribution[p]}]

Out[*] = n p


However, it quickly turned out that Mathematica just considers this a single variable named X[i].

I wonder if Expectation can be made to handle a variable number of random variables?

• For fixed n, you can always use Expectation[Sum[X[i], {i, 1, n}], Table[X[i] \[Distributed] BernoulliDistribution[p], {i, 1, n}]], naturally. I'm not sure there is a built-in way to make this work for symbolic n, though. – AccidentalFourierTransform Aug 14 '19 at 12:42
• You are seeking the expectation of the $2^\text{nd}$ sample central moment, which for any distribution whose first 2 moments exist, is equal to $\frac{ (n-1)}{n} \mu _2$, where $\mu_2$ denotes the $2^\text{nd}$ population central moment (i.e. the population variance). Your problem would be more interesting/challenging if it extended outside of textbook base cases. – wolfies Aug 14 '19 at 14:09
• @wolfies I was trying to keep it simple to focus on the Mathematica code rather than the mathematics. But indeed I'm actually trying to find the variance of the sample variance. (and a few other quantities). It seems simple enough to do by hand, but tedious and easy to get wrong. (My problem is somewhat related to stats.stackexchange.com/questions/196689 .) – Thomas Ahle Aug 14 '19 at 23:20

Clear["Global*"]


Use TransformedDistribution to find the distribution for the sum. Assuming that the x[i] are independent and identically distributed

dist[p_, n_Integer?Positive] :=
TransformedDistribution[Sum[x[i], {i, 1, n}],
Table[x[i] \[Distributed] BernoulliDistribution[p], {i, 1, n}]]

dist[p, #] & /@ Range

(* {BernoulliDistribution[p], BinomialDistribution[2, p],
BinomialDistribution[3, p], BinomialDistribution[4, p],
BinomialDistribution[5, p]} *)


The Mean is then

seq = Mean[dist[p, #]] & /@ Range

(* {p, 2 p, 3 p, 4 p, 5 p} *)


Generalizing,

mean[p_, n_] = FindSequenceFunction[seq, n]

(* n p *)


Similarly for the Variance

variance[p_, n_] = FindSequenceFunction[
Variance[dist[p, #]] & /@ Range, n]

(* n p - n p^2 *)

• So the approach is to calculate a number of values, and then use FindSequenceFunction to guess a closed form? In either case, that function is a catch! I had no idea Mathematica had something like that, even if the solutions don't prove anything. – Thomas Ahle Aug 14 '19 at 22:50
• What would you do if the x[i] weren't independent? Say I want to find E[Y] as a function of the variances and covariances of the x[i]? (Note in this case E[Y] = (1-1/n)\sigma^2, so the covariances don't matter, but for E[(Y-EY)^2] it does.) – Thomas Ahle Aug 14 '19 at 23:15
• Look at the documentation for TransformedDistribution. TransformedDistribution[expr, {x1, x2, ...} \[Distributed] dist] where dist` is the multivariate distribution. – Bob Hanlon Aug 14 '19 at 23:22