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“的枕函数，它是PDF的逆函数，所以它会返回特定概率密度的值，但是我还没有找到一个函数来返回比例参数。

实现迭代需要花费太多的时间(包括创建和计算)。我的脚本将是巨大的，所以我应该节省计算时间。在这种情况下，lambda函数有帮助吗？我大概知道它在做什么，但到目前为止我还没有使用它。对此有什么想法吗？

谢谢!

Answer 1

给出了正态概率密度函数，f

​

​

考虑到我们希望解决的f和x问题。让我们问一下同情，它是否能解出这个方程：

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函数，W的角度看，似乎有一个封闭的解决方案，它满足

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

所有复值z。

我们也可以使用渐近方法来找到给定的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的实值解决方案感兴趣。按照医生的说法，special.lambertw有两个部分真实的分支，当k=0和k=1.因此，上面的代码检查返回的值(对于这两个分支)是否为真，并返回任何实际解决方案的列表(如果它们存在)。如果不存在真正的解决方案，则返回一个空列表。如果对于任何西格玛的实际值(给定的y值)，都无法获得pdf值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小于y_max，有两种解决方案。