Python Numpy Poisson Distribution
Question: So basically 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)
    

So 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.

Thanks in advance.

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: 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 the poisson distribution is not even
approximatly identical to you 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](https://en.wikipedia.org/wiki/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](https://en.wikipedia.org/wiki/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`](http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.histogram.html)
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)

[![enter image description
here](http://i.stack.imgur.com/hDxf4.png)](http://i.stack.imgur.com/hDxf4.png)

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:

[![enter image description
here](http://i.stack.imgur.com/LLSza.png)](http://i.stack.imgur.com/LLSza.png)

