<p dir="auto">When you feed scipy.stats.poisson a mu &lt;= 0 (invalid according to <a href="http://en.wikipedia.org/wiki/Poisson_distribution" rel="nofollow">http://en.wikipedia.org/wiki/Poisson_distribution</a>) the outcome is that pmf or logpmf evaluate to nan for all data.  Using the Wikipedia definition, nonpositive mu should just raise an invalid value error.  This is what I'd probably expect.  Using a more general definition, poisson with zero variance could be a dirac delta distribution or something, but I don't think this generality is expected in scipy stats.</p>
<p dir="auto">Edit: because it's a discrete not continuous distribution the dirac delta wouldn't come into play for a generalized pmf, it would just have probability 1 for count=0 and probability 0 for all other counts.</p>
<p dir="auto">The Poisson PMF is well defined when mu is 0, it should be 0 for all k&gt;0 and 1 for k==0.</p>
<p dir="auto">But poisson.pmf fails (returns NaN) because it computes logpmf and then exponentiates:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="def _logpmf(self, k, mu):
    Pk = k*log(mu)-gamln(k+1) - mu
    return Pk

def _pmf(self, k, mu):
    return exp(self._logpmf(k, mu))

&gt;&gt;&gt; scipy.stats.poisson.pmf(0, 0)
nan
&gt;&gt;&gt; scipy.stats.poisson.pmf(1, 0)
nan 
&gt;&gt;&gt; mu = 0
&gt;&gt;&gt; k = 1
&gt;&gt;&gt; scipy.stats.poisson.pmf(k, mu)
nan
&gt;&gt;&gt; mu ** k * math.exp(-mu) / math.factorial(k)
0.0
&gt;&gt;&gt; k = 0
&gt;&gt;&gt; mu ** k * math.exp(-mu) / math.factorial(k)
1.0"><pre class="notranslate"><code class="notranslate">def _logpmf(self, k, mu):
    Pk = k*log(mu)-gamln(k+1) - mu
    return Pk

def _pmf(self, k, mu):
    return exp(self._logpmf(k, mu))

&gt;&gt;&gt; scipy.stats.poisson.pmf(0, 0)
nan
&gt;&gt;&gt; scipy.stats.poisson.pmf(1, 0)
nan 
&gt;&gt;&gt; mu = 0
&gt;&gt;&gt; k = 1
&gt;&gt;&gt; scipy.stats.poisson.pmf(k, mu)
nan
&gt;&gt;&gt; mu ** k * math.exp(-mu) / math.factorial(k)
0.0
&gt;&gt;&gt; k = 0
&gt;&gt;&gt; mu ** k * math.exp(-mu) / math.factorial(k)
1.0
</code></pre></div>