Python（Scipy）：查找高斯分布的尺度参数（标准差）

Question

在概率密度函数 (PDF) 中计算一个值的概率密度是很常见的。想象一下，我们有一个均值 = 40、标准差为 5 的高斯分布，现在想要得到值 32 的概率密度。我们会这样：

In [1]: import scipy.stats as stats
In [2]: print stats.norm.pdf(32, loc=40, scale=5)
Out [2]: 0.022

--> 概率密度为2.2%。

但是现在，让我们考虑一下逆问题。我有平均值，概率密度为 0.05，我想得到标准偏差（即比例参数）。

我可以实现的是一种数值方法：多次创建 stats.norm.pdf 并逐步增加比例参数，然后使用该参数使结果尽可能接近。

在我的例子中，我将值 30 指定为 5% 标记。所以我需要解决这个“方程”：

stats.norm.pdf(30, loc=40, scale=X) = 0.05

有一个叫做“ppf”的scipy函数，它是PDF的倒数，所以它会返回特定概率密度的值，但我还没有找到返回scale参数的函数 em>。

实施迭代会花费太多时间（包括创建和计算）。我的脚本会很大，所以我应该节省计算时间。在这种情况下，lambda 函数可以提供帮助吗？我大致知道它在做什么，但到目前为止我还没有使用它。对此有何想法？

谢谢！

Answer 1

normal probability density 函数 f 由下式给出

给定f 和x，我们希望求解?。请问sympy能不能解方程：

import sympy as sy
from sympy.abc import x, y, sigma

expr = (1/(sy.sqrt(2*sy.pi)*sigma) * sy.exp(-x**2/(2*sigma**2))) - y
ans = sy.solve(expr, sigma)[0]
print(ans)
# sqrt(2)*exp(LambertW(-2*pi*x**2*y**2)/2)/(2*sqrt(pi)*y)

因此，就LambertW function、W 而言，似乎存在一个封闭形式的解决方案，它满足

z = W(z) * exp(W(z))

对于所有复值z。

我们也可以使用 sympy 来找到给定 x 和 y 的数值结果，但是 也许做数值工作会更快 scipy.special.lambertw:

import numpy as np
import scipy.special as special

def sigma_func(x, y):
    results = set([np.real_if_close(
        np.sqrt(2)*np.exp(special.lambertw(-2*np.pi*x**2*y**2, k=k)/2)
        /(2*np.sqrt(np.pi)*y)).item() for k in (0, -1)])
    results = [s for s in results if np.isreal(s)]
    return results

一般来说，LambertW 函数返回复数，但我们只 对sigma 的实值解决方案感兴趣。 Per the docs, special.lambertw 有两个部分真实的分支，分别是 k=0 和 k=1。所以 上面的代码检查返回的值（对于这两个分支）是否是真实的，并且 如果存在，则返回任何实际解决方案的列表。如果没有真正的解决方案， 然后返回一个空列表。如果 pdf 值 y 不是 对于 sigma 的任何实际值（对于给定的 x 值）达到。

---

你可以这样使用它：

x = 30.0
loc = 40.0
y = 0.02
s = sigma_func(loc-x, y)
print(s)
# [16.65817044316178, 6.830458938511113]

import scipy.stats as stats
for si in s:
    assert np.allclose(stats.norm.pdf(x, loc=loc, scale=si), y)

---

在您给出的示例中，使用y = 0.025，sigma 没有解决方案：

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

x = 30.0
loc = 40.0
y = 0.025
s = np.linspace(5, 20, 100)
plt.plot(s, stats.norm.pdf(x, loc=loc, scale=s))
plt.hlines(y, 4, 20, color='red')  # the horizontal line y = 0.025
plt.ylabel('pdf')
plt.xlabel('sigma')
plt.show()

所以sigma_func(40-30, 0.025) 返回一个空列表：

In [93]: sigma_func(40-30, 0.025)
Out [93]: []

---

上面的情节是典型的，当y太大时，零 解决方案，在曲线的最大值处（我们称之为y_max）有一个 解决方案

In [199]: y_max = np.nextafter(np.sqrt(1/(np.exp(1)*2*np.pi*(10)**2)), -np.inf)

In [200]: y_max
Out[200]: 0.024197072451914336

In [201]: sigma_func(40-30, y_max)
Out[201]: [9.9999999776424]

对于小于 y_max 的 y，有两种解决方案。

Answer 2

这将是两个解决方案，因为普通 PDF 围绕均值对称。 就目前而言，您需要求解一个单变量方程。 它没有封闭形式的解决方案，因此您可以使用例如scipy.optimize.fsolve来解决。

编辑：请参阅 @unutbu 对 Lambert W 函数的封闭形式解决方案的回答。