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$ – AccidentalFourierTransform Aug 14 '19 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 Aug 14 '19 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$ – Thomas Ahle Aug 14 '19 at 23:20

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$ – Thomas Ahle Aug 14 '19 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$ – Thomas Ahle Aug 14 '19 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 Aug 14 '19 at 23:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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