SciPy 几何分布完全指南scipy.stats.geom 的数学定义、实现原理与实战用法【免费下载链接】scipySciPy library main repository项目地址: https://gitcode.com/gh_mirrors/sc/scipy几何分布Geometric Distribution是概率论中刻画首次成功所需试验次数的基础离散分布当每次独立试验的成功概率为p时随机变量k表示直到首次成功所进行的试验次数。本文以 几何分布官方教程 为骨架结合 SciPy 中scipy.stats.geom的源码实现scipy/stats/_discrete_distns.py完整推导其概率质量函数PMF、累积分布函数CDF、分位数、矩与矩生成函数并给出参数校验、数值稳定性设计、随机抽样与边界行为等源码级分析帮助读者在工程中正确、高效地使用该分布。一、几何分布的定义与背景几何分布属于 SciPy 离散统计分布家族描述伯努利试验序列中首次成功所需的试验次数。其核心设定为每次试验相互独立每次试验的成功概率恒为p失败概率为1 - p随机变量k表示获得一次成功所需的试验次数取值k ≥ 1。与首次成功前失败的次数即 Planck 分布scipy.stats.planck不同几何分布从 1 开始计数。例如掷硬币p0.5连续抛掷直到出现正面抛掷次数即为几何分布随机变量。参数p的合法区间为p ∈ (0, 1)在源码中由geom_gen._argcheck强制校验# scipy/stats/_discrete_distns.py def _argcheck(self, p): return (p 1) (p 0)即只接受0 p 1区间端点 0 被拒绝、1 被接受p1时必然首次即成功退化为退化分布。二、核心数学公式全解1. 概率质量函数PMF$$ p(k; p) (1-p)^{k-1} p, \quad k \geq 1 $$含义前k-1次全部失败概率(1-p)^{k-1}第k次成功概率p。源码实现为def _pmf(self, k, p): return np.power(1-p, k-1) * p def _logpmf(self, k, p): return special.xlog1py(k - 1, -p) log(p)其中对数概率使用了special.xlog1py即x·log(1y)的数值稳定实现避免1-p在p接近 1 时发生灾难性抵消保障logpmf的精度。2. 累积分布函数CDF$$ F(x; p) 1 - (1-p)^{\lfloor x \rfloor}, \quad x \geq 1 $$由于是离散分布CDF 对非整数参数x取下整。源码实现利用expm1与log1p保证小概率端数值稳定性def _cdf(self, x, p): k floor(x) return -expm1(log1p(-p)*k) def _logsf(self, x, p): k floor(x) return k*log1p(-p) def _sf(self, x, p): return np.exp(self._logsf(x, p))生存函数SF 1 − CDF同样采用log1pexp的稳定路线避免大k下直接计算1 - (1-p)^k的舍入误差。3. 分位数函数PPF / 逆 CDF$$ G(q; p) \left\lceil \frac{\log(1-q)}{\log(1-p)} \right\rceil $$对给定概率q ∈ [0, 1]返回满足F(k) ≥ q的最小整数k。源码在取整后额外做了一次 CDF 回查校正处理边界与浮点误差def _ppf(self, q, p): vals ceil(log1p(-q) / log1p(-p)) temp self._cdf(vals-1, p) return np.where((temp q) (vals 0), vals-1, vals)4. 矩与矩生成函数统计量公式备注均值 μ1/p成功概率越小期望试验次数越多方差 μ₂(1-p)/p²偏度 γ₁(2-p)/√(1-p)恒为正右偏峰度 γ₂(p²-6p6)/(1-p)超额峰度excess kurtosis矩生成函数 M(t)p / (e⁻ᵗ − (1-p))定义域需满足e⁻ᵗ 1-p矩在源码_stats中实现峰度系数用多项式求值np.polyval([1, -6, 6], p)计算def _stats(self, p): mu 1.0/p qr 1.0-p var qr / p / p g1 (2.0-p) / sqrt(qr) g2 np.polyval([1, -6, 6], p)/(1.0-p) return mu, var, g1, g2四个返回值依次对应mean、variance、skew、kurtosis超额峰度与上表完全一致。三、源码级实现剖析1. 类定义与分布实例化几何分布封装在 scipy/stats/_discrete_distns.py 的geom_gen类中继承自rv_discrete基类模块底部实例化geom geom_gen(a1, namegeom, longnameA geometric)参数a1声明了支撑集下界最小取值为 1即首次试验就成功的情形这与 CDF 公式中x ≥ 1的条件一一对应。支撑集上界默认取无穷因为理论上k可以任意大连续失败的次数没有上限。2. 形状参数声明_shape_info声明了形状参数p的定义域为开区间(0, 1)def _shape_info(self): return [_ShapeInfo(p, False, (0, 1), (True, True))]该信息用于通用参数校验、数组广播与自动生成文档保证调用geom.pmf(k, p)时p必须满足0 p ≤ 1。3. 随机数生成与 int64 边界处理_rvs直接调用 NumPy 的geometric采样并做了一处关键加固——将因溢出回绕成负数的采样值替换为int64最大值def _rvs(self, p, sizeNone, random_stateNone): res random_state.geometric(p, sizesize) # RandomState.geometric can wrap around to negative values; make behavior # consistent with Generator.geometric by replacing with maximum integer. max_int np.iinfo(res.dtype).max return np.where(res 0, max_int, res)这一点与geom_gen文档字符串中记录的边界行为呼应当p 10⁻¹⁷时观测值超过np.iinfo(np.int64).max的概率迅速上升当p 10⁻²⁰时几乎全部观测值都会超过该上限。由于输出 dtype 固定为int64这些值会被裁剪到最大值。实际使用时若非专门研究极端小概率场景建议保持p ≥ 10⁻¹⁷。4. 熵Entropy源码额外提供了闭式熵公式def _entropy(self, p): return -np.log(p) - np.log1p(-p) * (1.0-p) / p该式由几何分布熵的级数求和得出log1p的使用再次体现对小概率参数p → 0场景的数值稳健设计。四、实战用法示例1. 概率质量与累积概率import numpy as np from scipy import stats # 成功概率 0.3 的几何分布 p 0.3 geom stats.geom(p) # 首次成功恰好在第 4 次试验的概率 print(geom.pmf(4)) # 0.7**3 * 0.3 0.1029 # 在前 4 次试验内成功的概率CDF print(geom.cdf(4)) # 1 - 0.7**4 ≈ 0.7599 # 生存函数4 次试验仍未成功的概率 print(geom.sf(4)) # 0.7**4 0.2401 # 逆 CDF90% 概率必然成功所需的试验次数 print(geom.ppf(0.9)) # ceil(log(0.1)/log(0.7)) 7 # 描述性统计量 print(geom.mean()) # 1/0.3 ≈ 3.333 print(geom.var()) # 0.7/0.09 ≈ 7.778 print(geom.stats(momentsmvsk)) # (均值, 方差, 偏度, 峰度)2. 随机抽样# 生成 10000 个几何分布样本seed 固定以便复现 rvs geom.rvs(size10000, random_state42) print(rvs.min(), rvs.mean()) # 均值应接近 1/p ≈ 3.33 # 指定 p 的向量化调用不同 p 对应不同样本 ps np.array([0.1, 0.5, 0.9]) samples stats.geom.rvs(ps, size1000, random_state0)3. 与相关离散分布的对照stats.planck几何分布的从 0 计数变体geometric planck 1二者互为平移关系stats.nbinom负二项分布当失败次数参数n1时负二项分布退化为几何分布。几何分布可视为负二项分布在n1的特例stats.bernoulli单次试验的 0/1 结果是构成几何分布序列的基本试验单元。4. 对数域计算的推荐用法对于尾概率计算优先使用logpmf与logsf避免极小概率下溢# p 很小、k 很大时直接 pmf 可能下溢为 0 log_p stats.geom.logpmf(1000, 0.001) # ≈ -1.0005 print(np.exp(log_p)) # 精确恢复 ≈ 3.68e-2五、数值稳定性设计要点从源码可归纳出geom实现的三大数值稳健策略源码依据见 scipy/stats/_discrete_distns.pylog1p/expm1组合所有形如(1-p)^k的计算都改写为exp(k·log1p(-p))形式避免p接近 1 时1-p的有效数字损失对数域优先_logpmf、_logsf作为一等公民实现_sf由exp(_logsf)得出而非1 - _cdf规避大数相减的灾难性抵消分位数回查校正_ppf在向上取整后再次调用_cdf验证保证离散分位数单调性在浮点边界处不被破坏。历史版本曾修复geom.logpmf(1, 1)返回nan的问题应返回0.0相关记录见 0.15.0 版本发布说明本次源码中的_argcheck允许p1正是该修复的直接体现。六、应用场景小结几何分布广泛用于首次成功等待时间类问题可靠性工程设备连续无故障运行周期数、首次故障前的检查次数质量控制抽检流水线中首次发现次品所需的检验件数通信与随机接入ALOHA 类随机多址协议中节点首次成功发送所需的时隙数随机游走与排队论首次命中特定状态所需的步数近似。使用 scipy.stats.geom 时只需牢记三点参数p ∈ (0,1]是成功概率随机变量从 1 开始计数当p小于10⁻¹⁷时注意抽样结果的int64裁剪边界。其余统计推断、矩计算与向量化调用均可直接借助rv_discrete基类的完整方法族完成。【免费下载链接】scipySciPy library main repository项目地址: https://gitcode.com/gh_mirrors/sc/scipy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
