Poisson分布绘图异常求助：大数计算溢出致r≥8后图形错误

你的问题核心是直接计算mu1^r和factorial(r)时出现大数溢出，导致后续概率计算失真。下面提供几种可靠的解决办法：

方法1：使用scipy.stats的泊松PMF函数（推荐） #

scipy内置的泊松分布概率质量函数（PMF）已经做了数值稳定性处理，能避免大数溢出问题，直接调用即可：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import poisson, norm

# 参数设置
mu1 = 20
sigma = np.sqrt(20)
x1 = np.arange(1, 21, 1)

# 计算理论概率值
# 高斯分布：用norm.pdf，参数loc=均值，scale=标准差
gauss_vals = norm.pdf(x1, loc=mu1, scale=sigma)
# 泊松分布：用poisson.pmf，参数mu=均值
poisson_vals = poisson.pmf(x1, mu=mu1)

# 绘图
plt.bar(x1, poisson_vals, color='b', alpha=1, label='Poisson-Verteilung')
plt.bar(x1, gauss_vals, color='r', alpha=0.2, label='Gauss-Verteilung')
plt.xlim(-0.1, 20.1)
plt.xticks(x1)
plt.xlabel('r')
plt.ylabel('Wahrscheinlichkeit')
plt.title(r'Poisson- und Gauss-Verteilung für $(\mu = 20; r = 1,2,...,20)$')
plt.legend()
plt.show()

方法2：手动用对数计算避免溢出 #

如果不想依赖scipy的stats模块，可以通过对数转换，把乘法转为加法，避免直接计算超大数：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial

sigma = np.sqrt(20)
mu1 = 20

def fGauss(r):
    return 1/(np.sqrt(2*np.pi)*sigma) * np.exp(-(r - mu1)**2 / (2*sigma**2))

def fPoisson_stable(r):
    # 取对数计算：ln(mu^r / r!) = r*ln(mu) - ln(r!)
    log_prob = r * np.log(mu1) - np.log(factorial(r))
    # 再加上ln(exp(-mu)) = -mu，最后指数还原
    return np.exp(log_prob - mu1)

x1 = np.arange(1, 21, 1)
Gauss = [fGauss(r) for r in x1]
Poisson2 = [fPoisson_stable(r) for r in x1]

plt.bar(x1, Poisson2, color='b', alpha=1, label='Poisson-Verteilung')
plt.bar(x1, Gauss, color='r', alpha=0.2, label='Gauss-Verteilung')
plt.xlim(-0.1, 20.1)
plt.xticks(x1)
plt.xlabel('r')
plt.ylabel('Wahrscheinlichkeit')
plt.title(r'Poisson- und Gauss-Verteilung für $(\mu = 20; r = 1,2,...,20)$')
plt.legend()
plt.show()

关于np.random.poisson的说明 #

你朋友提到的np.random.poisson是用来生成泊松分布的随机样本，而不是计算理论概率值。如果需要用样本直方图来近似泊松分布，可以这样用：

# 生成10000个泊松分布样本
poisson_samples = np.random.poisson(mu=mu1, size=10000)
# 绘制直方图
plt.hist(poisson_samples, bins=np.arange(0.5, 21.5, 1), color='b', alpha=0.7, label='Poisson样本直方图', density=True)
# 叠加高斯分布曲线
x = np.linspace(1,20,100)
plt.plot(x, norm.pdf(x, loc=mu1, scale=sigma), 'r-', label='Gauss-Verteilung')
plt.xlim(-0.1,20.1)
plt.xticks(x1)
plt.xlabel('r')
plt.ylabel('相对频率/概率密度')
plt.title('Poisson样本分布与Gauss-Verteilung')
plt.legend()
plt.show()

这种方式得到的是样本的近似分布，适合展示统计特性，但如果你需要的是精确的理论概率值，还是用前两种方法更合适。

内容的提问来源于stack exchange，提问作者NoobOfAllTrades -