如何在 python 中 return 通过 sklearn 的函数 KernelDensity 估计的分布的平均值（或期望值）？
How to return mean value (or expectation value) of a distribution estimated via function KernelDensity of sklearn, in python?
我的问题是，如何return估计"kde"的均值和方差？或者有没有你知道的任何其他包可以轻松输出平均值或方差值，如 print kde.mean() 或 print kde.get_parameter(mean)?
import numpy as np
from scipy.stats import norm
from sklearn.neighbors import KernelDensity
N = 100
np.random.seed(1)
X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
kde = KernelDensity(kernel='gaussian', bandwidth=0.5).fit(X)
一般来说，你需要用数字来做这个。我建议两种不同的方法：
- 整合
- Monte Carlo 模拟
这些方法适用于任何 内核 和任何 带宽 。
整合
利用一旦我们知道概率密度函数这一事实，我们就可以通过积分轻松计算均值和方差。
请注意，在 scikit-learn 中，方法 score_samples returns 记录 pdf，因此需要 "exp" 它。
Monte Carlo 模拟
这里的想法是简单地从 KDE 中抽样并通过样本均值和方差估计总体均值和方差。
代码
import numpy as np
from scipy.integrate import quad
from sklearn.neighbors import KernelDensity
N = 100
np.random.seed(1)
X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
kde = KernelDensity(kernel='gaussian', bandwidth=0.5).fit(X)
# Mean and Variance - Integration
pdf = lambda x : np.exp(kde.score_samples([[x]]))[0]
mean_integration = quad(lambda x: x * pdf(x), a=-np.inf, b=np.inf)[0]
variance_integration = quad(lambda x: (x ** 2) * pdf(x), a=-np.inf, b=np.inf)[0] - mean_integration ** 2
# Mean and Variance - Monte Carlo
n_samples = 10000000
samples = kde.sample(n_samples)
mean_mc = samples.mean()
variance_mc = samples.var()
print('Mean:\nIntegration: {}\nMonte Carlo: {}\n'.format(mean_integration, mean_mc))
print('Variance\nIntegration: {}\nMonte Carlo: {}\n'.format(variance_integration, variance_mc))
输出：
Mean:
Integration: 3.560582852075697
Monte Carlo: 3.5595633705830934
Variance:
Integration: 6.645066811078639
Monte Carlo: 6.646732489654485
原则
我 运行 遇到了同样的情况，你可以在 statsmodels 库中找到期望的工作实现，但是有趣的是它正在修补 scipy 而不是做同样的功能在它自己的图书馆里。
他们似乎有 expect()、expect_v2() 和 expect_discrete()，它们在另一个使用 scipy.integrate.quad 进行数值积分的答案中使用了类似的技术。
TL;DR 版本
此外，如果您像我一样懒惰并使用 statsmodels 推断的默认采样，您可以通过将数值积分替换为 density 的点积来快速获得期望值，并且support。大多数情况下可能 'close enough'。
import numpy as np
import statsmodels.api as sm
import scipy
N = 100
np.random.seed(1)
X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
kde = sm.nonparametric.KDEUnivariate(X)
kde.fit(kernel='gau', bw=0.5)
mean1 = np.dot(kde.density, kde.support) / kde.density.sum()
mean2 = scipy.integrate.quad(lambda x: kde.evaluate(x) * x, kde.support[0], kde.support[-1])
print('TL;DR version - Mean:', mean1)
print('Integration version - Mean:', mean2)
print('TL;DR version - Variance:', np.dot(kde.density, kde.support**2) / kde.density.sum() - mean1**2)
print('Integration version - Variance:', scipy.integrate.quad(lambda x: kde.evaluate(x) * x**2, kde.support[0], kde.support[-1])[0] - mean2[0]**2)
输出：
TL;DR version - Mean: 3.5605148164179368
Integration version - Mean: (3.5604536291684905, 1.9311947816995413e-08)
TL;DR version - Variance: 6.646077637181225
Integration version - Variance: 6.644042199345121
四边形和点的区别
请注意，由于缺少 evaluate()，quad 方法不适用于 Gaussian 以外的内核。比如你fit的时候传kernel='epa', fft=False，目前只能用numpy.dot。这可能会在未来的版本中改变。
作为一个有趣的效果，seaborn 库依靠 statsmodels 来显示它是 kde 的，但是由于上述原因只适用于 Gaussian。对于 cdf，即使指定了不同的 cdf，它也会显示高斯 cdf； pdf 似乎仍然遵守 kerenel 规范。
参考
如果存储库更改 [1 月 22 日的最新提交 bfa3e69]，也会在此处转载，但是如果将此算法与 statsmodels 一起使用，您可能必须将 pdf 调用替换为 evaluate 来电。
我的问题是，如何return估计"kde"的均值和方差？或者有没有你知道的任何其他包可以轻松输出平均值或方差值，如 print kde.mean() 或 print kde.get_parameter(mean)?
import numpy as np
from scipy.stats import norm
from sklearn.neighbors import KernelDensity
N = 100
np.random.seed(1)
X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
kde = KernelDensity(kernel='gaussian', bandwidth=0.5).fit(X)
一般来说，你需要用数字来做这个。我建议两种不同的方法：
- 整合
- Monte Carlo 模拟
这些方法适用于任何 内核 和任何 带宽 。
整合
利用一旦我们知道概率密度函数这一事实，我们就可以通过积分轻松计算均值和方差。
请注意，在 scikit-learn 中，方法 score_samples returns 记录 pdf，因此需要 "exp" 它。
Monte Carlo 模拟
这里的想法是简单地从 KDE 中抽样并通过样本均值和方差估计总体均值和方差。
代码
import numpy as np
from scipy.integrate import quad
from sklearn.neighbors import KernelDensity
N = 100
np.random.seed(1)
X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
kde = KernelDensity(kernel='gaussian', bandwidth=0.5).fit(X)
# Mean and Variance - Integration
pdf = lambda x : np.exp(kde.score_samples([[x]]))[0]
mean_integration = quad(lambda x: x * pdf(x), a=-np.inf, b=np.inf)[0]
variance_integration = quad(lambda x: (x ** 2) * pdf(x), a=-np.inf, b=np.inf)[0] - mean_integration ** 2
# Mean and Variance - Monte Carlo
n_samples = 10000000
samples = kde.sample(n_samples)
mean_mc = samples.mean()
variance_mc = samples.var()
print('Mean:\nIntegration: {}\nMonte Carlo: {}\n'.format(mean_integration, mean_mc))
print('Variance\nIntegration: {}\nMonte Carlo: {}\n'.format(variance_integration, variance_mc))
输出：
Mean: Integration: 3.560582852075697 Monte Carlo: 3.5595633705830934
Variance: Integration: 6.645066811078639 Monte Carlo: 6.646732489654485
原则
我 运行 遇到了同样的情况，你可以在 statsmodels 库中找到期望的工作实现，但是有趣的是它正在修补 scipy 而不是做同样的功能在它自己的图书馆里。
他们似乎有 expect()、expect_v2() 和 expect_discrete()，它们在另一个使用 scipy.integrate.quad 进行数值积分的答案中使用了类似的技术。
TL;DR 版本
此外，如果您像我一样懒惰并使用 statsmodels 推断的默认采样，您可以通过将数值积分替换为 density 的点积来快速获得期望值，并且support。大多数情况下可能 'close enough'。
import numpy as np
import statsmodels.api as sm
import scipy
N = 100
np.random.seed(1)
X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]
X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]
kde = sm.nonparametric.KDEUnivariate(X)
kde.fit(kernel='gau', bw=0.5)
mean1 = np.dot(kde.density, kde.support) / kde.density.sum()
mean2 = scipy.integrate.quad(lambda x: kde.evaluate(x) * x, kde.support[0], kde.support[-1])
print('TL;DR version - Mean:', mean1)
print('Integration version - Mean:', mean2)
print('TL;DR version - Variance:', np.dot(kde.density, kde.support**2) / kde.density.sum() - mean1**2)
print('Integration version - Variance:', scipy.integrate.quad(lambda x: kde.evaluate(x) * x**2, kde.support[0], kde.support[-1])[0] - mean2[0]**2)
输出：
TL;DR version - Mean: 3.5605148164179368
Integration version - Mean: (3.5604536291684905, 1.9311947816995413e-08)
TL;DR version - Variance: 6.646077637181225
Integration version - Variance: 6.644042199345121
四边形和点的区别
请注意，由于缺少 evaluate()，quad 方法不适用于 Gaussian 以外的内核。比如你fit的时候传kernel='epa', fft=False，目前只能用numpy.dot。这可能会在未来的版本中改变。
作为一个有趣的效果，seaborn 库依靠 statsmodels 来显示它是 kde 的，但是由于上述原因只适用于 Gaussian。对于 cdf，即使指定了不同的 cdf，它也会显示高斯 cdf； pdf 似乎仍然遵守 kerenel 规范。
参考
如果存储库更改 [1 月 22 日的最新提交 bfa3e69]，也会在此处转载，但是如果将此算法与 statsmodels 一起使用，您可能必须将 pdf 调用替换为 evaluate 来电。