# sampling distribution of the mean for arbitrary 1-D pdf

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.

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.
• (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'). – lacerbi Jun 23 '15 at 15:09
• 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. – Glen_b -Reinstate Monica Jun 24 '15 at 1:18