# Unexpected under-dispersion in Perl simulations of Poisson RV

The following piece of Perl code randomly maps a set of ranges onto a circumference of a circle. In the example, the circumference is of length 1000 and legal ranges are e.g. (0,8)=0,1,2,...,8 and (995,2)=995,996,...,999,0,1,2 (i.e. zero-based coordinates; both start and end are inclusive).

I take some arbitrary position on the circumference (e.g. 36) and count how many ranges cover it in each simulation.

finally, I calculate the mean and variance of this statistic.

use strict;
use warnings;

use Statistics::Descriptive;

my $n_simulations = 1000; my$circumference   = 1000;
my @lengths_distrib = (100) x 100;    # distibution of range lengths
my $some_pos = 36; # arbitrary position my$stat = Statistics::Descriptive::Full->new();

foreach my $sim ( 1 ..$n_simulations ) {

# randomly map ranges onto circumference

my @random_ranges =
map { my $start = int( rand($circumference) ); [ $start, ($start + $_ -1 ) %$circumference ] }
@lengths_distrib;

# count how many range contain $some_pos my$num_covering_ranges = scalar(
grep { ( $_->[0] <=$some_pos and $_->[1] >=$some_pos ) or ( $_->[1] <$_->[0] and $_->[1] >$some_pos ) }
@random_ranges
);
$stat->add_data($num_covering_ranges);
}
print $stat->mean, ' ',$stat->variance, "\n";


To the best of my knowledge, this kind of random variable should follow Poisson distribution (law of rare events and so on). Hence, the mean and variance should be equal. However, the variance seems to systematically be a bit lower than the mean.

What am I missing?

UPDATE Following whuber response, a couple of notes and additions:

1. the Poisson distribution is a good approximation of the binomial distribution if n is at least 20 and p is smaller than or equal to 0.05, and an excellent approximation if n ≥ 100 and np ≤ 10. source: wikipedia. The example I gave here is a toy one. I usually use much larger circumferences (length ~ 3M), many more ranges (around 25k) an relatively small ranges (size ~ 2.5k), so I think in general Poisson distribution should be appropriate. I will give it a look again on some real data.

2. In my real data, I have e.g. some 25k ranges of different sizes. Currently, I'm simulating random mapping of these clones on the circumference, count the coverage of some point of interest in each simulation (this could be an arbitrary fixed point since all point are equal), then get the mean of this coverage over some 100 or 1000 simulations.

Now, I use a CDF of a Poisson distribution with this mean to get a P-value for the actual observation (number of ranges covering a specific point of interest). I take lower tail since I'm interested in under-covered positions.

How would you suggest going around this? Should I use a Binomial distribution? How can I evaluate its parameters? Note this is not exactly a Binomial RV, since while I do have n independent trials, the chances for each trial to succeed is not equal (it depends on the range size). Perhaps I should simply take the average success probability (average length) and not even have to use simulations, right? In other words, I guess what I'm asking is how does the sum of independent Bernoulli RV with different p's look like?

I should also note that while this example only deals with the number of ranges covering a single position, I'm also interested in the number of ranges covering some range. For now I use the simulations scheme exactly the same way described before.

Thanks!

• Whether the Poisson distribution is a "good" approximation depends on your purpose. In this situation, a general rule of thumb sheds no light. You have noted a discrepancy; one obvious place to look for an explanation lies in using the Poisson as an approximation to the correct distribution, which is Binomial (that's not in question, I hope). Another place, indicated in your update, lies in using ranges of "different sizes": now all bets are off. You have a mixture of Binomials, which should create some overdispersion: en.wikipedia.org/wiki/Overdispersion
– whuber
Nov 30, 2010 at 17:26
• @whuber Indeed, I have a mixture of Binomials, and it seems sometimes the Poisson approximation is not very good. My question is what are my options for getting a p-value for an actual observation. One way would be to run many simulations and just get an empirical p-value. This is the simplest way I know, but it uses nothing of the info I have -- I do now it's a mixture of Binomials and I know all their p's -- isn't there a way to use this valuable information? Nov 30, 2010 at 19:13

If I understood correctly, a "range" of length $k$ (such as $k=100$) within a circumference of length $n$ (such as $n = 1000$) has a chance of $(2k+1)/n$ of covering a given point on the circumference, and all the chances in a simulation are independent. Therefore in a simulation of $N$ trials (such as $N=1000$) the count of ranges that cover a given point should follow a $\text{Binomial}(\frac{2k+1}{n}, N)$ distribution. This distribution has mean $(2k+1)\frac{N}{n}$ and variance $(2k+1)(1 - \frac{2k+1}{n})\frac{N}{n}$. The variance of your data (the simulation trial outcomes) will deviate from this latter value due to chance.

These considerations identify two reasons for your observations:

1. The variance ought to be a fraction of the mean; specifically, it should only be about $1 - \frac{2k+1}{n}$ of the mean. In the case $k=100$ and $n=1000$, the variance should only be around 79.9% of the mean.

2. When you conduct a small number of trials (and $N=1000$ is quite small for simulations of this type), expect noticeable chance deviations between the results and the expectations.

• Why (2k+1)/n? I think it's k/n. For example, if my range is of length k=10 and my given point is 50, the range must fall on [41,50] or [42,51] or [43,52] or ... [50,60] in order to cover this point (10 options). The simplest case - the range is of length k=1 - it has only one option to cover the point (fall exactly on it). Nov 30, 2010 at 9:07
• And in the general case, when we look at the probability a randomly mapped range of length k covers a given range of length m on a circumference of length g, we get p=(k-m+1)/n. When m=1, as in the original post, this is indeed k/n. Nov 30, 2010 at 9:12
• @David You're right, thank you. I was (confusedly) equating "length" with "radius." But this does not render the explanation any less valid. The question is whether the discrepancies you are seeing are as predicted (by the corrected formulas).
– whuber
Nov 30, 2010 at 15:58