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?

  • 1
    $\begingroup$ 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. $\endgroup$ Commented Aug 14, 2019 at 12:42
  • 1
    $\begingroup$ 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. $\endgroup$
    – wolfies
    Commented Aug 14, 2019 at 14:09
  • $\begingroup$ @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 .) $\endgroup$ Commented Aug 14, 2019 at 23:20

1 Answer 1


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[5]

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

The Mean is then

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

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


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

(* n p *)

Similarly for the Variance

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

(* n p - n p^2 *)
  • $\begingroup$ 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. $\endgroup$ Commented Aug 14, 2019 at 22:50
  • $\begingroup$ 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.) $\endgroup$ Commented Aug 14, 2019 at 23:15
  • $\begingroup$ Look at the documentation for TransformedDistribution. TransformedDistribution[expr, {x1, x2, ...} \[Distributed] dist] where dist is the multivariate distribution. $\endgroup$
    – Bob Hanlon
    Commented Aug 14, 2019 at 23:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.