使用 scipy.stats 拟合非标准化分布
Fitting an un-normalised distribution with scipy.stats
我正在尝试拟合直方图,但拟合仅适用于归一化数据,即在直方图中使用选项 normed=True
。有没有办法用 scipy 统计数据(或其他方法)做到这一点?这是一个使用均匀分布的 MWE:
import matplotlib.pyplot as plt
import numpy as np
import random
from scipy.stats import uniform
data = []
for i in range(1000):
data.append(random.uniform(-1,1))
loc, scale = uniform.fit(data)
x = np.linspace(-1,1, 1000)
y = uniform.pdf(x, loc, scale)
plt.hist(data, bins=100, normed=False)
plt.plot(x, y, 'r-')
plt.show()
我也尝试过定义自己的函数(如下),但我觉得不太合适。
import matplotlib.pyplot as plt
import numpy as np
import random
from scipy import optimize
data = []
for i in range(1000):
data.append(random.uniform(-1,1))
def unif(x,avg,sig):
return avg*x + sig
y, base = np.histogram(data,bins=100)
x = [0.5 * (base[i] + base[i+1]) for i in xrange(len(base)-1)]
popt, pcov = optimize.curve_fit(unif, x, y)
x_fit = np.linspace(x[0], x[-1], 100)
y_fit = unif(x_fit, *popt)
plt.hist(data, bins=100, normed=False)
plt.plot(x_fit, y_fit, 'r-')
plt.show()
请注意,将分布拟合到直方图通常不是一个好主意。与原始数据相比,直方图包含的信息较少,因此拟合很可能更差。因此,问题中的第一个 MWE 实际上包含最佳方法。简单地标准化直方图,它将匹配数据的分布:plt.hist(data, bins=100, normed=True)
.
但是,您似乎真的想使用非标准化直方图。在这种情况下,采用直方图通常使用的归一化并将其 inverted 应用于拟合分布。 documentation 将规范化描述为
n/(len(x)`dbin)
这是冗长的说法 除以观察次数乘以 bin 宽度。
将分布乘以该值得到每个 bin 的预期计数:
loc, scale = uniform.fit(data)
x = np.linspace(-1,1, 1000)
y = uniform.pdf(x, loc, scale)
n_bins = 100
bin_width = np.ptp(data) / n_bins
plt.hist(data, bins=n_bins, normed=False)
plt.plot(x, y * len(data) * bin_width, 'r-')
第二个 MWE 很有趣,因为您描述了 不合适 行,但实际上它是 非常合适 :)。您只是 过度拟合 直方图,因为尽管您期望一条水平线(一个自由度),但您拟合了一条任意线(两个自由度)。
所以如果你想要一条水平线适合水平线并且如果你适合其他东西也不要惊讶得到其他东西...
def unif(x, sig):
return 0 * x + sig # slope is zero -> horizontal line
但是,有一种更简单的方法来获取非标准化均匀分布的高度。只需对所有 bin 的直方图进行平均:
y, base = np.histogram(data,bins=100)
y_hat = np.mean(y)
print(y_hat)
# 10.0
或者,更简单的使用len(data) / n_bins == 10
的理论值。
我正在尝试拟合直方图,但拟合仅适用于归一化数据,即在直方图中使用选项 normed=True
。有没有办法用 scipy 统计数据(或其他方法)做到这一点?这是一个使用均匀分布的 MWE:
import matplotlib.pyplot as plt
import numpy as np
import random
from scipy.stats import uniform
data = []
for i in range(1000):
data.append(random.uniform(-1,1))
loc, scale = uniform.fit(data)
x = np.linspace(-1,1, 1000)
y = uniform.pdf(x, loc, scale)
plt.hist(data, bins=100, normed=False)
plt.plot(x, y, 'r-')
plt.show()
我也尝试过定义自己的函数(如下),但我觉得不太合适。
import matplotlib.pyplot as plt
import numpy as np
import random
from scipy import optimize
data = []
for i in range(1000):
data.append(random.uniform(-1,1))
def unif(x,avg,sig):
return avg*x + sig
y, base = np.histogram(data,bins=100)
x = [0.5 * (base[i] + base[i+1]) for i in xrange(len(base)-1)]
popt, pcov = optimize.curve_fit(unif, x, y)
x_fit = np.linspace(x[0], x[-1], 100)
y_fit = unif(x_fit, *popt)
plt.hist(data, bins=100, normed=False)
plt.plot(x_fit, y_fit, 'r-')
plt.show()
请注意,将分布拟合到直方图通常不是一个好主意。与原始数据相比,直方图包含的信息较少,因此拟合很可能更差。因此,问题中的第一个 MWE 实际上包含最佳方法。简单地标准化直方图,它将匹配数据的分布:plt.hist(data, bins=100, normed=True)
.
但是,您似乎真的想使用非标准化直方图。在这种情况下,采用直方图通常使用的归一化并将其 inverted 应用于拟合分布。 documentation 将规范化描述为
n/(len(x)`dbin)
这是冗长的说法 除以观察次数乘以 bin 宽度。
将分布乘以该值得到每个 bin 的预期计数:
loc, scale = uniform.fit(data)
x = np.linspace(-1,1, 1000)
y = uniform.pdf(x, loc, scale)
n_bins = 100
bin_width = np.ptp(data) / n_bins
plt.hist(data, bins=n_bins, normed=False)
plt.plot(x, y * len(data) * bin_width, 'r-')
第二个 MWE 很有趣,因为您描述了 不合适 行,但实际上它是 非常合适 :)。您只是 过度拟合 直方图,因为尽管您期望一条水平线(一个自由度),但您拟合了一条任意线(两个自由度)。
所以如果你想要一条水平线适合水平线并且如果你适合其他东西也不要惊讶得到其他东西...
def unif(x, sig):
return 0 * x + sig # slope is zero -> horizontal line
但是,有一种更简单的方法来获取非标准化均匀分布的高度。只需对所有 bin 的直方图进行平均:
y, base = np.histogram(data,bins=100)
y_hat = np.mean(y)
print(y_hat)
# 10.0
或者,更简单的使用len(data) / n_bins == 10
的理论值。