我有遵循高斯分布的数据。然而，数据仅对于一系列值 [xa,xb] 来说才是真正的高斯分布，所以我想使用以下方法拟合截断正态分布scipy.stats.truncnorm同时利用我知道范围 [xa,xb] 的事实。我的目标是找到地点和规模。
我不明白如何在拟合中修复 xa 和 xb 。形状参数是“a”和“b”，但这些参数取决于位置和比例，这是我的未知数。而且，似乎无法对'a'和'b'进行初步猜测（它们只能用fa和fb冻结？）。当我做：
par = truncnorm.fit(r, a=a_guess, b=b_guess, scale= scale_guess, loc = loc_guess)
I get
未知参数：{'a'：0.0，'b'：2.4444444444444446}。
而且，我得到的配合非常不稳定。这是一个例子：
from scipy.stats import truncnorm
import matplotlib.pyplot as plt
xa, xb = 30,250
loc, loc_guess = 50, 30
scale, scale_guess = 75, 90
a,b = (xa-loc)/scale, (xb-loc)/scale
fig, ax = plt.subplots(1, 1)
x = np.linspace(xa,xb,10000)
ax.plot(x, truncnorm.pdf(x, a, b, loc=loc, scale=scale),
'r-', lw=5, alpha=0.6, label='truncnorm pdf')
r = truncnorm.rvs(a, b, loc=loc, scale=scale, size=10000)
par = truncnorm.fit(r, scale= scale_guess, loc = loc_guess)
ax.plot(x, truncnorm.pdf(x, *par),
'b-', lw=1, alpha=0.6, label='truncnorm fit')
ax.hist(r, density=True, histtype='stepfilled', alpha=0.3)
plt.legend()
plt.show()
第一个例子 第二个例子
我也经常收到这样的警告：
/home/elie/anaconda2/envs/py36/lib/python3.6/site-packages/scipy/stats/_continuous_distns.py:5823：RuntimeWarning：日志中遇到除以零 self._logdelta = np.log(self._delta)
正如您所发现的，问题在于您想要保持固定的参数，xa and xb，不是本机参数truncnorm. truncnorm具有形状参数a and b，通过设置 x 间隔来确定形状standard正态分布。然后这个形状被移动和缩放loc and scale参数。关系是
xa = a*scale + loc
xb = b*scale + loc
To fix xa and xb，您可以使用接受等式约束的 SciPy 最小化器之一。这里我将使用scipy.optimize.fmin_slsqp。 （您可以使用“综合”功能scipy.optmize.minimize，其中包括 SLSQP 求解器作为其选项之一。）
这是一个演示如何使用的脚本fmin_slsqp对于这个问题。功能func是要最小化的目标函数。它只是一个包装truncnorm.nnlf，负对数似然函数。功能constraint返回包含两个值的数组。当满足约束时这些值为 0。
import numpy as np
from scipy.stats import truncnorm
from scipy.optimize import fmin_slsqp
import matplotlib.pyplot as plt
def func(p, r, xa, xb):
return truncnorm.nnlf(p, r)
def constraint(p, r, xa, xb):
a, b, loc, scale = p
return np.array([a*scale + loc - xa, b*scale + loc - xb])
xa, xb = 30, 250
loc = 50
scale = 75
a = (xa - loc)/scale
b = (xb - loc)/scale
# Generate some data to work with.
r = truncnorm.rvs(a, b, loc=loc, scale=scale, size=10000)
loc_guess = 30
scale_guess = 90
a_guess = (xa - loc_guess)/scale_guess
b_guess = (xb - loc_guess)/scale_guess
p0 = [a_guess, b_guess, loc_guess, scale_guess]
par = fmin_slsqp(func, p0, f_eqcons=constraint, args=(r, xa, xb),
iprint=False, iter=1000)
xmin = 0
xmax = 300
x = np.linspace(xmin, xmax, 1000)
fig, ax = plt.subplots(1, 1)
ax.plot(x, truncnorm.pdf(x, a, b, loc=loc, scale=scale),
'r-', lw=3, alpha=0.4, label='truncnorm pdf')
ax.plot(x, truncnorm.pdf(x, *par),
'k--', lw=1, alpha=1.0, label='truncnorm fit')
ax.hist(r, bins=15, density=True, histtype='stepfilled', alpha=0.3)
ax.legend(shadow=True)
plt.xlim(xmin, xmax)
plt.grid(True)
plt.show()
这是它生成的图。样本数据是随机的，因此每次运行的图都会不同。
注意：有时会生成随机数据集fmin_slsqp计算期间因“遇到无效值”而失败。我没有进一步研究这个问题，但您的数据可能会遇到这个问题。