Never use scipy.stats!!

benchmark.ipy

import scipy.stats as ss
import math
import numpy as np
from timeit import timeit

np_log_poisson = lambda k, l: -l + k * np.log(l) - math.lgamma(k+1)
log_poisson = lambda k, l: -l + k * math.log(l) - math.lgamma(k+1)

%timeit -n 10000 ss.poisson.logpmf(5, 4)  # 164.0 us
%timeit -n 10000 ss.poisson._logpmf(5, 4) #  23.1 us
%timeit -n 10000 np_log_poisson(5, 4)     #  11.2 us
%timeit -n 10000 log_poisson(5, 4)        #   1.2 us

(to be pendantic about it)