'Return predicted probability of each count level for each observation
Parameters
n : array-like or int
The counts for which you want the probabilities. If n is None
then the probabilities for each count from 0 to max(y) are
given.
Returns
ndarray
A nobs x n array where len(`n`) columns are indexed by the count
n. If n is None, then column 0 is the probability that each
observation is 0, column 1 is the probability that each
observation is 1, etc.'

def predict_prob(self, n=None, exog=None, exposure=None, offset=None, transform=True):

    if (n is not None):
        counts = np.atleast_2d(n)
    else:
        counts = np.atleast_2d(np.arange(0, (np.max(self.model.endog) + 1)))
    mu = self.predict(exog=exog, exposure=exposure, offset=offset, transform=transform, linear=False)[:, None]
    return stats.poisson.pmf(counts, mu)
