I want to compute the sampling distribution of the mean for $k$ samples from an arbitrary, known probability density $f(x)$, with $x \in \mathbb{R}$. What is the most efficient way to do so numerically?

(For the record, this is part of a MATLAB toolbox I am writing, but the question is language-independent.)

This is what I thought so far on top of my head.

Let $f_k$ be the sampling distribution of the mean for $k$ samples ($f_1 = f$). Then we can write: $$ f_k(x) = \int \int \delta\left(x - \frac{(k-1) \cdot y + z}{k}\right) \; f_{k-1}(y) \; f(z) \, dy \, dz \\ = \int f_{k-1}(y) \; f\left(k \cdot x - (k-1) \cdot y\right) dy \\ \propto \int f_{k-1}\left(\frac{y}{k-1}\right) \; f\left(k \cdot x - y\right) dy. \\ $$ By appropriately defining two auxiliary densities: $$g_k (x) \propto f_k(k \cdot x) \\ h_k(y) \propto f_k(y/(k-1)),$$ the equation for $f_k$ can be rewritten as: $$ g_k(\tilde{x}) \propto \int h_{k-1}(y) \; f\left(\tilde{x} - y\right) dy, $$ which is a convolution, and therefore can be efficiently solved by FFT.

This way, I can write an algorithm that iterates its way up to any $k$. The disadvantage of this method is that, in addition to the FFT, at every step I need to recompute the $g_k$ and $h_k$ (e.g., via numerical interpolation).

The naive algorithm would compute every $k$ iteratively, with total computation time $O(k)$, but it just occurred to me that with a divide et impera approach I need only to compute $k$ for powers of 2 and then fill in the rest, with computation time $O(\log_2 k)$. I haven't fledged this out but it should work pretty much the same.

My questions are:

  1. Is there already a "standard" algorithm or a package (in any language - e.g., R) to solve this problem? I'd be curious to see how they approach it.
  2. Can you think of a better way than the one I am proposing? (I am not particularly happy with the $g_k$ and $h_k$, but I don't see a way around it.)

PS: This is my first post on any stackexchange, so I am happy to take suggestions on how to improve the format.


1 Answer 1


I don't know about "most efficient", but consider that the kth power of the Fourier transform of $f$ is the FT of the sum of the $k$ random variates. So once you have the FT of $f$, the remaining calculations should involve only a power and the inverse transform.

See also Wikipedia on the Characteristic function (which up to sign of the argument of $\exp$) is just a Fourier transform.

  • $\begingroup$ Yes, that's a good point.This is essentially what I was doing in a less efficient way. $\endgroup$
    – lacerbi
    Commented Jun 23, 2015 at 15:03
  • $\begingroup$ (pressed enter by mistake, this is the full comment) Yes, that's a very good point, thanks! This is essentially what I was doing, but my method was less efficient and less clear. The only downside is that if I do a naive direct computation of the $k$-th power via FT I would need to allocate $O(k)$ space (otherwise the tails of the distribution would get 'clipped'). $\endgroup$
    – lacerbi
    Commented Jun 23, 2015 at 15:09
  • $\begingroup$ Yes, you need to have enough room for the solution, which is "larger" than the thing you start with. If your original $f$ is done on a grid of $m$ bins, you may need as many as $km$ bins, depending on how rapidly the tail goes down and how large $k$ is. $\endgroup$
    – Glen_b
    Commented Jun 24, 2015 at 1:18

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.