Generate poisson numbers which sum a certain value?

Hi I want to create a vector with numbers following a poisson distribution but at the same time the sum of these numbers should be the same that a certain value.

r= 750;
while true
dis= poissrnd(lambda,1,1000);   
if(sum(dis)== r)
 break;
end
end
sum(dis)

This code works sometimes but usually it takes a lot of time to run it so, there is another way to do the same but in a less demanding iteration?

Answer 1

[John D'Errico]

I frequently see this happen, where someone wants to generate random numbers, subject to some constraint, that essentially invalidates one or more of the requirements for those random numbers.

In this case, the set will arguably not be truly random Poisson distributed, if they sum to some fixed value, because they are not independent. The argument is a subtle one, but not one I will bother to take on. The real problem is the events will not be independently Poisson, once that constraint comes into play.

At least, you need to be careful in how you generate them. Rejection methods, as you currently are doing, are a classically viable scheme, yielding random sets subject to a constraint. The problem is that rejection tends to be inefficient, generating many sets of numbers that are usually wasted, since they will fail to meet the requirement.

Instead, we might want to look for a constructive approach to such a problem. So what does a Poisson random number mean? It is often used to model the number of random events that fall in a region of some space, usually the real line.

https://en.wikipedia.org/wiki/Poisson_point_process

So, think of a process that generates points on the real line. If you want a total of N Poisson samples, subject to the constraint that they sum to S?

In your example, the sum would be 750, but you wish to find a set of 1000 Poisson random numbers, such that they sum to 750.

I recall the mean of a poisson distribution with parameter lambda is lambda. You can check that here:

https://en.wikipedia.org/wiki/Poisson_distribution

But if we try a random set, it misses.

lambda = 750/1000;
P = poissrnd(lambda,1,1000);
sum(P)
ans =
   726

The Poisson distribution for small parameter is pretty skewed, so that did pretty poorly. If I take a sufficient number of samples, eventually I might get lucky and hit 750, but slow it will be.

Again, the point is we can use a different stochastic process. So, instead, generate 750 points, uniformly on the unit interval.

Now, just count the number of events that lie in each sub-interval, of length 1/1000. We will expect to see on average 750/1000 such events to lie in each sub-interval. The total number of events will be 750 of course. So the effective Poisson parameter will be 750/1000. The solution is as simple as:

S = 750;
nevents= 1000;
u = rand(S,1);
P = accumarray(ceil(u*nevents),1,[nevents,1]);

Is the sum correct? Of course. It must be so by construction, as well as the number of Poisson random numbers generated.

sum(P)
ans =
   750

There are 1000 events, also as desired.

numel(P)
ans =
   1000

And the rate parameter is roughly as expected.

[LAMBDAHAT, LAMBDACI] = poissfit(P)
LAMBDAHAT =
         0.751503006012024
LAMBDACI =
         0.697719614918322
         0.805286397105726

P
ans =
     0
     0
     3
     0
     1
     3
     0
     1
     1
     1
     0
     2
     1
     1
     2
     0
     1
     1
     0
     0
...

The resulting vector P should follow effectively a Poisson distribution, but subject to the desired sum constraint. No rejection was needed at all. It has been a while since I worked with Poisson processes, but this should be a viable solution to your problem. And it will be fast.

[John D'Errico]

Hmm. I need to think about whether exponential inter-arrivals are required. It would be trivially fixed though. It has been about 40 years since I worked with these things, so I'd need to do some reading.

[pauldjn]

Thanks sir this answer was more than I expected I guess I completely forgot the essentials of the distribution