I'm trying to compute some p-values for samples from a distribution of sums of ~1000 random variables. The exact distribution of these random variables isn't known, but I have empirical estimates that I think are pretty accurate.

So far I've been using the central limit theorem to produce a normal approximation for this sum, which does OK since n is relatively large, but not great. Computing the exact distribution of the sum via convolutions is too slow, but I only really care about having high accuracy in the tails; everywhere else a rough approximation is fine.

Are there any methods that will allow me to improve my estimates in the tails (or just one tail) without having to compute the entire convolutions? I'm not sure if this would be a variant on the CLT implementation or something completely different. My gut feeling is that there isn't any way to do this, so I'd be open to any kind of solution at all!

  • $\begingroup$ Where do the empirical estimates come from and in what form are they? (In a brief test Mathematica computed the convolution of 1000 discrete RVs, each specified by ten parameters, in just a few seconds, so I don't understand why computing the exact distribution is so slow.) $\endgroup$ – whuber Nov 30 '10 at 16:42
  • $\begingroup$ @whuber They're stored as histograms, i.e. lists of frequencies. My problem is that I have to ~20000 of these computations in a reasonable time, so a few seconds is unfortunately not fast enough. $\endgroup$ – bnaul Nov 30 '10 at 18:40
  • 1
    $\begingroup$ Summarize each histogram with a low-order truncation of its cumulant generating function and add those functions to obtain an Edgeworth-type expansion of the sum (en.wikipedia.org/wiki/Edgeworth_series ). The computational effort is proportional to the total number of bars in all those histograms. $\endgroup$ – whuber Nov 30 '10 at 19:21
  • $\begingroup$ Have you considered calculating the convolution by using (fast) Fourier transform (and its inverse) ? $\endgroup$ – Andre Holzner Dec 1 '10 at 20:36

$n=1000$ gets you actually extremely close to normal if you do not go to the extreme tails ($n=10$ is often close enough to normal in the central region), but if you need to go there, estimating a PDF might well be impossible.

To get an event very far out in the tails, extreme fluctuations in your sumands need to happen. Since you say your distributions were "empirical estimates" you should ask yourself if you have enough input data to understand the tails of your summand distributions. In many applications empirical distributions are being truncated in the tails, since they are not well understood and might be the only contributors in some parameter regions.

No matter if I truncate my distributions or not, I often use Monte Carlo simulations to compute convolutions of (not parametrized) empirical distributions since computing power is cheap.

  • $\begingroup$ Thanks for the answer. These are significances of large sets of genome regions; typically a hypergeometric or binomial test is used in cases like these, and p-values on the order of 10^-50 to 100 are not at all uncommon (if you pass in a set of regions that turn out to all affect limb development, you get a pretty small p-value for a "limb" test). I initially dismissed Monte Carlo as too slow, but I don't care too much about speed so I will definitely experiment some with that. $\endgroup$ – bnaul Nov 30 '10 at 6:02
  • $\begingroup$ @bnaul: No idea how people work in your field, but contemplating differences on a 10^-50 or 10^-100 level?? Even with n=1000 I image that would be hard to back up with actual data. Hope you people can sleep well ;) $\endgroup$ – Benjamin Bannier Nov 30 '10 at 6:41

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.