从百分位数查找概率分布参数

John D. Cook

4.18/5 (8投票s)

2010年2月3日

BSD

4分钟阅读

74785

668

给定两个百分位数约束，如何确定概率分布的参数。

下载源代码 - 1.43 KB

引言

假设某件事有 10% 的概率小于 30，90% 的概率小于 60。一旦您选择了概率分布族（正态分布、伽马分布等），您就需要一种确定哪些参数将满足您的两个需求的方法。

更确切地说，假设随机变量 X 有一个两参数分布。问题是找到这些参数的值，使得 Pr(X < x₁) = p₁ 且 Pr(X < x₂) = p₂。

本文将展示如何使用 Python 的 SciPy 库计算正态分布、柯西分布、威布尔分布、伽马分布和逆伽马分布的这些参数。

查找满足两个百分位数方程的参数的问题是实际的。例如，在确定贝叶斯统计中的先验分布时，就会出现这个问题。本文也作为使用 SciPy 的简要介绍。

背景

本文提出的问题的数学解决方案在技术报告从分位数确定分布参数中进行了描述。在这里，我们解释如何在 Python 中实现该报告中的计算。这是一篇关于该问题的更多动机的文章。

Using the Code

此处描述的代码非常简单易于调用。假设您想找到正态分布的均值和标准差。您只需使用适当的参数调用 normal_parameters。它返回作为一对的均值和标准差。

(mean, stdev) = normal_parameters(x1, p1, x2, p2)

其他分布的函数也类似

cauchy_parameters
weibull_parameters
gamma_parameters
inverse_gamma_parameters
beta_parameters

所有函数都接受相同的四个参数，并且都返回两个参数。但是，这些参数的含义对于每个分布都是不同的。对于正态分布，问题是找到 μ 和 σ，使得 F(x₁) = p₁ 且 F(x₂) = p₂，其中 F(x) 定义为

$F(x) = \frac{1}{\sqrt{2\pi}\sigma} \int_{-\infty}^x \exp\left( -\frac{(t-\mu)^2}{2\sigma^2}\right) \, dt$

对于柯西分布，问题类似，除了现在

$F(x) = \frac{1}{\pi\sigma} \int_{-\infty}^x \frac{dt}{ 1 + \left(\frac{t-\mu}{\sigma} \right)^2 }}$

对于威布尔分布，问题是找到 γ 和 β，其中

$F(x) = \frac{\gamma}{\beta^\gamma} \int_0^x t^{\gamma-1} \exp\left(-(t/\beta)^\gamma \right)\, dt$

对于伽马分布，问题是找到 α 和 β，其中

$F(x) = \frac{1}{\Gamma(\alpha)\,\beta^\alpha} \int_0^x t^{\alpha-1} \exp(-x/\beta)\, dt$

对于逆伽马分布，问题是找到 α 和 β，其中

$F(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} \int_0^x t^{-\alpha-1}\exp(-\beta/t)\, dt$

对于贝塔分布，问题是找到 a 和 b，其中

$F(x) = \frac{\Gamma(a+b)}{\Gamma(a)\,\Gamma(b)} \int_0^x t^{a-1}(1-t)^{b-1}\, dt$

查看代码内部

此处提供的代码严重依赖于 SciPy，这是一个用于 Python 科学计算的大型库。有关 SciPy 的介绍，请参阅 CodeProject 文章开始使用 SciPy (科学 Python) 库。

对于正态分布和柯西分布，位置参数由下式给出

$\frac{x_1 F^{-1}(p_2) - x_2 F^{-1}(p_1)}{F^{-1}(p_2) - F^{-1}(p_1)}$

并且尺度参数由下式给出

$\frac{x_2 - x_1}{F^{-1}(p_2) - F^{-1}(p_1)}$

其中 F(x) 是正态分布或柯西分布的 CDF，如上一节所述。CDF 的逆由 SciPy 中的 ppt 方法给出。这里“ppt”代表“百分点函数”。其他库可能会将其称为分位数函数。概率分布类位于 scipy.stats 中。

威布尔分布的参数可以通过一个简单的公式给出，不需要任何 SciPy 功能。逆伽马参数也很容易找到，因为逆伽马问题可以简化为找到伽马分布参数的问题。

伽马分布参数不能如此简单地获得。用于查找这些参数的代码说明了 SciPy 和 functools 模块。

def gamma_parameters(x1, p1, x2, p2):
    
    # Standardize so that x1 < x2 and p1 < p2
    if p1 > p2:
        (p1, p2) = (p2, p1)
        (x1, x2) = (x2, x1)
    
    # function to find roots of for gamma distribution parameters
    def objective(alpha):
        return stats.gamma.ppf(p2, alpha) / stats.gamma.ppf(p1, alpha) - x2/x1
    
    # The objective function we're wanting to find a root of is decreasing.
    # We need to find an interval over which is goes from positive to negative.
    left = right = 1.0
    while objective(left) < 0.0:
        left /= 2
    while objective(right) > 0.0:
        right *= 2
    alpha = optimize.bisect(objective, left, right)
    beta = x1 / stats.gamma.ppf(p1, alpha)
    
    return (alpha, beta)

请注意，辅助函数 objective 定义在 gamma_parameters 函数内部。Python 的一个很好的特性是，您可以使用闭包轻松定义这样的小函数。只有 gamma_parameters 需要这个函数，因此我们将其定义在 gamma_parameters 内部，以避免弄乱外部作用域。

要找到函数 objective 为 0 的位置，我们使用 scipy.optimize 中的一种求根方法。此处的代码使用了 bisect 方法，与其他求根算法相比，该算法速度较慢，但非常可靠。

对于贝塔分布，我们有一种直接的方法来求解分布参数。相反，我们使用 SciPy 函数 fmin 求解一个优化问题。我们搜索最接近满足百分位数约束的参数。实际上，这些参数将完全满足约束。

def beta_parameters(x1, p1, x2, p2):
    "Find parameters for a beta random variable X 
	so that P(X > x1) = p1 and P(X > x2) = p2."

    def square(x): 
        return x*x

    def objective(v):
        (a, b) = v
        temp  = square( stats.beta.cdf(x1, a, b) - p1 )
        temp += square( stats.beta.cdf(x2, a, b) - p2 )
        return temp
    
    # arbitrary initial guess of (3, 3) for parameters
    xopt = optimize.fmin(objective, (3, 3))
    return (xopt[0], xopt[1])

用于计算贝塔分布参数的方法也可以用于其他分布。实际上，我们可以将这种方法用于所有分布。但是，用于其他分布的专用方法更准确、更高效。

历史

2010 年 2 月 3 日：初始版本
2010 年 2 月 4 日：文章更新
2010 年 7 月 26 日：添加贝塔分布