Python Numpy Poisson Distribution

Question

I am generating a Gaussian, for the sake of completeness, that's my implementation:

from numpy import *
x=linspace(0,1,1000)
y=exp(-(x-0.5)**2/(2.0*(0.1/(2*sqrt(2*log(2))))**2))

with peak at 0.5 and fwhm=0.1. So far so not interesting. In the next step I calculate the poisson distribution of my set of data using numpys random.poisson implementation.

poi = random.poisson(lam=y)

I'm having two major problems.

1. A specialty of poisson is that the variance equals the exp. value, comparing the output of mean() and var() does confuse me as the outputs are not equal.
2. When plotting this, the poisson dist. takes up integer values only and the max. value is around 7, sometimes 6, whilst my old function y has its max. at 1. Afai understand, the poisson-function should give me sort of a 'fit' of my actual function y. How come the max. values are not equal? Sorry for my mathematical incorrectness, actually I'm doing this to emulate poisson-distributed noise but I guess you understand 'fit' in this context.

---

EDIT: 3. question: What's the 'size' variable used for in this context? I've seen different types of usage but in the end they did not give me different results but failing when choosing it wrong...

EDIT2: OK, from the answer I got I think that I was not clear enough (although it already helped me correct some other stupid errors I did, thanks for that!). What I want to do is apply poisson (white) noise to the function y. As described by MSeifert in the post below, I now use the expectation value as lam. But this only gives me the noise. I guess I have some understanding problems on the level of how th{is,e} noise is applied (and maybe it's more physics related?!).

Answer 1

First of all, I'll write this answer assuming you import numpy as np because it clearly distinguishes numpy functions from the builtins or those of the math and random package of python.

I think it is not necessary to answer your specified questions because your basic assumption is wrong:

Yes, the poisson-statistics has a mean that equals the variance but that assumes you use a constant lam. But you don't. You input the y-values of your gaussian, so you cannot expect them to be constant (they are by your definition gaussian!).

Use np.random.poisson(lam=0.5) to get one random value from a poisson distribution. But be careful since this poisson distribution is not even approximately identical to your gaussian distribution because you are in the "low-mean" interval where both of these are significantly different, see for example the Wikipedia article about Poisson distribution.

Also you are creating random numbers, so you shouldn't really plot them but plot a np.histogram of them. Since statistical distributions are all about probabilitiy density functions (see Probability density function).

Before, I already mentioned that you create a poisson distribution with a constant lam so now it is time to talk about the size: You create random numbers, so to approximate the real poisson distribution you need to draw a lot of random numbers. There the size comes in: np.random.poisson(lam=0.5, size=10000) for example creates an array of 10000 elements each drawn from a poissonian probability density function for a mean value of 0.5.

And if you haven't read it in the Wikipedia article mentioned before the poisson distribution gives by definition only unsigned (>= 0) integer as result.

So I guess what you wanted to do is create a gaussian and poisson distribution containing 1000 values:

gaussian = np.random.normal(0.5, 2*np.sqrt(2*np.log(2)), 1000)
poisson = np.random.poisson(0.5, 1000)

and then to plot it, plot the histograms:

import matplotlib.pyplot as plt
plt.hist(gaussian)
plt.hist(poisson)
plt.show()

or use the np.histogram instead.

To get statistics from your random samples you can still use np.var and np.mean on the gaussian and poisson samples. And this time (at least on my sample run) they give good results:

print(np.mean(gaussian))
0.653517935138
print(np.var(gaussian))
5.4848398775
print(np.mean(poisson))
0.477
print(np.var(poisson))
0.463471

Notice how the gaussian values are almost exactly what we defined as parameters. On the other hand poisson mean and var are almost equal. You can increase the precision of the mean and var by increasing the size above.

Why the poisson distribution doesn't approximate your original signal

Your original signal contains only values between 0 and 1, so the poisson distribution only allows positive integer and the standard deviation is linked to the mean value. So far from the mean of the gaussian your signal is approximatly 0, so the poisson distribution will almost always draw 0. Where the gaussian has it's maximum the value is 1. The poisson distribution for 1 looks like this (left is the signal + poisson and on the right the poisson distribution around a value of 1)

so you'll get a lot of 0 and 1 and some 2 in that region. But also there is some probability that you draw values up to 7. This is exactly the antisymmetry that I mentioned. If you change the amplitude of your gaussian (multiply it by 1000 for example) the "fit" is much better since the poisson distribution is almost symmetric there: