Python: 两个函数之间的重叠（kde 和 normal 的 PDF）
Python: Overlap between two functions (PDF of kde and normal)
简短摘要：我正在尝试弄清楚如何计算两个函数之间的重叠。一个是高斯分布，另一个是基于数据的核密度。然后，我想做一个小算法来选择高斯的均值和方差，最大化重叠
首先，需要导入：
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats.kde import gaussian_kde
import scipy
我有一些大致正常的数据（右尾有点重）。我计算了这个数据的核密度、cdf 和 pdf （在这个例子中，数据是从制服中提取的，因为我无法提供真实数据） 像这样：
def survivalFunction():
data = np.random.normal(7,1,100) #Random data
p = sns.kdeplot(data, shade=False, lw = 3)
x,y = p.get_lines()[0].get_data()
cdf = scipy.integrate.cumtrapz(y, x, initial=0)
plt.hist(data,50,normed = 1,facecolor='b',alpha = 0.3)
然后我有另一个函数，它只是一个简单的高斯：
def surpriseFunction(mu,variance):
hStates = np.linspace(0,20,100)
sigma = math.sqrt(variance)
plt.plot(hStates,scipy.stats.norm.pdf(hStates, mu, sigma))
调用函数
surpriseFunction(5,1)
survivalFunction()
给出这个情节
您可能已经注意到，交换不同的 mu 值，围绕法线移动或多或少地与内核估计重叠。现在，我的问题是双重的：
1) 如何计算两个函数之间的重叠？
2) 我将如何制作一个小算法，为高斯选择均值和方差，以最大化重叠？
好的，所以我进行了相当大的改组，我认为它分离了主要部分，并且可以很容易地在各种功能中进行模块化。我之前给出的答案的原始代码是 here.
这是新内容，希望它是不言自明的。
# Setup our various global variables
population_mean = 7
population_std_dev = 1
samples = 100
histogram_bins = 50
# And setup our figure.
from matplotlib import pyplot
fig = pyplot.figure()
ax = fig.add_subplot(1,1,1)
from numpy.random import normal
hist_data = normal(population_mean, population_std_dev, samples)
ax.hist(hist_data, bins=histogram_bins, normed=True, color="blue", alpha=0.3)
from statsmodels.nonparametric.kde import KDEUnivariate
kde = KDEUnivariate(hist_data)
kde.fit()
#kde.supprt and kde.density hold the x and y values of the KDE fit.
ax.plot(kde.support, kde.density, color="red", lw=4)
#Gaussian function - though you can replace this with something of your choosing later.
from numpy import sqrt, exp, pi
r2pi = sqrt(2*pi)
def gaussian(x, mu, sigma):
return exp(-0.5 * ( (x-mu) / sigma)**2) / (sigma * r2pi)
#interpolation of KDE to produce a function.
from scipy.interpolate import interp1d
kde_func = interp1d(kde.support, kde.density, kind="cubic", fill_value=0)
你想做的只是标准曲线拟合——有很多方法可以做到，你说你想通过最大化两个函数的重叠来拟合曲线（为什么？）。 curve_fir scipy 例程是最小二乘法拟合，它试图最小化两个函数之间的差异 - 差异是微妙的：最大化重叠不会惩罚拟合函数大于数据，而 curve_fit 会。
我已经包含了使用这两种技术的解决方案，并对它们进行了分析：
#We need to *maximise* the overlap integral
from scipy.integrate import quad as integrate
def overlap(func1, func2, limits, func1_args=[], func2_args=[]):
def product_func(x):
return min(func1(x, *func1_args),func2(x, *func2_args))
return integrate(product_func, *limits)[0] # we only care about the absolute result for now.
limits = hist_data.min(), hist_data.max()
def gaussian_overlap(args):
mu, sigma = args
return -overlap(kde_func, gaussian, limits, func2_args=[mu, sigma])
现在两种不同的方法，重叠度量：
import cProfile, pstats, StringIO
pr1 = cProfile.Profile()
pr1.enable()
from scipy.optimize import fmin_powell as minimize
mu_overlap_fit, sigma_overlap_fit = minimize(gaussian_overlap, (population_mean, population_std_dev))
pr1.disable()
s = StringIO.StringIO()
sortby = 'cumulative'
ps = pstats.Stats(pr1, stream=s).sort_stats(sortby)
ps.print_stats()
print s.getvalue()
3122462 function calls in 6.298 seconds
Ordered by: cumulative time
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 6.298 6.298 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:2120(fmin_powell)
1 0.000 0.000 6.298 6.298 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:2237(_minimize_powell)
57 0.000 0.000 6.296 0.110 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:279(function_wrapper)
57 0.000 0.000 6.296 0.110 C:\Users\Will\Documents\Python_scripts\hist_fit.py:47(gaussian_overlap)
57 0.000 0.000 6.296 0.110 C:\Users\Will\Documents\Python_scripts\hist_fit.py:39(overlap)
57 0.000 0.000 6.296 0.110 C:\Python27\lib\site-packages\scipy\integrate\quadpack.py:42(quad)
57 0.000 0.000 6.295 0.110 C:\Python27\lib\site-packages\scipy\integrate\quadpack.py:327(_quad)
57 0.069 0.001 6.295 0.110 {scipy.integrate._quadpack._qagse}
66423 0.154 0.000 6.226 0.000 C:\Users\Will\Documents\Python_scripts\hist_fit.py:41(product_func)
4 0.000 0.000 6.167 1.542 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:2107(_linesearch_powell)
4 0.000 0.000 6.166 1.542 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:1830(brent)
4 0.000 0.000 6.166 1.542 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:1887(_minimize_scalar_brent)
4 0.001 0.000 6.166 1.542 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:1717(optimize)
和scipy方法curve_fit:
pr2 = cProfile.Profile()
pr2.enable()
from scipy.optimize import curve_fit
(mu_curve_fit, sigma_curve_fit), _ = curve_fit(gaussian, kde.support, kde.density, p0=(population_mean, population_std_dev))
pr2.disable()
s = StringIO.StringIO()
sortby = 'cumulative'
ps = pstats.Stats(pr2, stream=s).sort_stats(sortby)
ps.print_stats()
print s.getvalue()
122 function calls in 0.001 seconds
Ordered by: cumulative time
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.001 0.001 C:\Python27\lib\site-packages\scipy\optimize\minpack.py:452(curve_fit)
1 0.000 0.000 0.001 0.001 C:\Python27\lib\site-packages\scipy\optimize\minpack.py:256(leastsq)
1 0.000 0.000 0.001 0.001 {scipy.optimize._minpack._lmdif}
19 0.000 0.000 0.001 0.000 C:\Python27\lib\site-packages\scipy\optimize\minpack.py:444(_general_function)
19 0.000 0.000 0.000 0.000 C:\Users\Will\Documents\Python_scripts\hist_fit.py:29(gaussian)
1 0.000 0.000 0.000 0.000 C:\Python27\lib\site-packages\scipy\linalg\basic.py:314(inv)
1 0.000 0.000 0.000 0.000 C:\Python27\lib\site-packages\scipy\optimize\minpack.py:18(_check_func)
您可以看到 curve_fit 方法快很多，结果：
from numpy import linspace
xs = linspace(-1, 1, num=1000) * sigma_overlap_fit * 6 + mu_overlap_fit
ax.plot(xs, gaussian(xs, mu_overlap_fit, sigma_overlap_fit), color="orange", lw=2)
xs = linspace(-1, 1, num=1000) * sigma_curve_fit * 6 + mu_curve_fit
ax.plot(xs, gaussian(xs, mu_curve_fit, sigma_curve_fit), color="purple", lw=2)
pyplot.show()
非常相似。我会推荐 curve_fit。在这种情况下，速度提高了 6000 倍。当底层数据更复杂时，区别是 一点点 更多，但差别不大，你仍然可以获得巨大的加速。这是一个适合 6 个均匀分布的正态分布的示例：
选择 curve_fit!
简短摘要：我正在尝试弄清楚如何计算两个函数之间的重叠。一个是高斯分布，另一个是基于数据的核密度。然后，我想做一个小算法来选择高斯的均值和方差，最大化重叠
首先，需要导入：
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats.kde import gaussian_kde
import scipy
我有一些大致正常的数据（右尾有点重）。我计算了这个数据的核密度、cdf 和 pdf （在这个例子中，数据是从制服中提取的，因为我无法提供真实数据） 像这样：
def survivalFunction():
data = np.random.normal(7,1,100) #Random data
p = sns.kdeplot(data, shade=False, lw = 3)
x,y = p.get_lines()[0].get_data()
cdf = scipy.integrate.cumtrapz(y, x, initial=0)
plt.hist(data,50,normed = 1,facecolor='b',alpha = 0.3)
然后我有另一个函数，它只是一个简单的高斯：
def surpriseFunction(mu,variance):
hStates = np.linspace(0,20,100)
sigma = math.sqrt(variance)
plt.plot(hStates,scipy.stats.norm.pdf(hStates, mu, sigma))
调用函数
surpriseFunction(5,1)
survivalFunction()
给出这个情节
您可能已经注意到，交换不同的 mu 值，围绕法线移动或多或少地与内核估计重叠。现在，我的问题是双重的：
1) 如何计算两个函数之间的重叠？
2) 我将如何制作一个小算法，为高斯选择均值和方差，以最大化重叠？
好的，所以我进行了相当大的改组，我认为它分离了主要部分，并且可以很容易地在各种功能中进行模块化。我之前给出的答案的原始代码是 here.
这是新内容，希望它是不言自明的。
# Setup our various global variables
population_mean = 7
population_std_dev = 1
samples = 100
histogram_bins = 50
# And setup our figure.
from matplotlib import pyplot
fig = pyplot.figure()
ax = fig.add_subplot(1,1,1)
from numpy.random import normal
hist_data = normal(population_mean, population_std_dev, samples)
ax.hist(hist_data, bins=histogram_bins, normed=True, color="blue", alpha=0.3)
from statsmodels.nonparametric.kde import KDEUnivariate
kde = KDEUnivariate(hist_data)
kde.fit()
#kde.supprt and kde.density hold the x and y values of the KDE fit.
ax.plot(kde.support, kde.density, color="red", lw=4)
#Gaussian function - though you can replace this with something of your choosing later.
from numpy import sqrt, exp, pi
r2pi = sqrt(2*pi)
def gaussian(x, mu, sigma):
return exp(-0.5 * ( (x-mu) / sigma)**2) / (sigma * r2pi)
#interpolation of KDE to produce a function.
from scipy.interpolate import interp1d
kde_func = interp1d(kde.support, kde.density, kind="cubic", fill_value=0)
你想做的只是标准曲线拟合——有很多方法可以做到，你说你想通过最大化两个函数的重叠来拟合曲线（为什么？）。 curve_fir scipy 例程是最小二乘法拟合，它试图最小化两个函数之间的差异 - 差异是微妙的：最大化重叠不会惩罚拟合函数大于数据，而 curve_fit 会。
我已经包含了使用这两种技术的解决方案，并对它们进行了分析：
#We need to *maximise* the overlap integral
from scipy.integrate import quad as integrate
def overlap(func1, func2, limits, func1_args=[], func2_args=[]):
def product_func(x):
return min(func1(x, *func1_args),func2(x, *func2_args))
return integrate(product_func, *limits)[0] # we only care about the absolute result for now.
limits = hist_data.min(), hist_data.max()
def gaussian_overlap(args):
mu, sigma = args
return -overlap(kde_func, gaussian, limits, func2_args=[mu, sigma])
现在两种不同的方法，重叠度量：
import cProfile, pstats, StringIO
pr1 = cProfile.Profile()
pr1.enable()
from scipy.optimize import fmin_powell as minimize
mu_overlap_fit, sigma_overlap_fit = minimize(gaussian_overlap, (population_mean, population_std_dev))
pr1.disable()
s = StringIO.StringIO()
sortby = 'cumulative'
ps = pstats.Stats(pr1, stream=s).sort_stats(sortby)
ps.print_stats()
print s.getvalue()
3122462 function calls in 6.298 seconds
Ordered by: cumulative time
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 6.298 6.298 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:2120(fmin_powell)
1 0.000 0.000 6.298 6.298 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:2237(_minimize_powell)
57 0.000 0.000 6.296 0.110 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:279(function_wrapper)
57 0.000 0.000 6.296 0.110 C:\Users\Will\Documents\Python_scripts\hist_fit.py:47(gaussian_overlap)
57 0.000 0.000 6.296 0.110 C:\Users\Will\Documents\Python_scripts\hist_fit.py:39(overlap)
57 0.000 0.000 6.296 0.110 C:\Python27\lib\site-packages\scipy\integrate\quadpack.py:42(quad)
57 0.000 0.000 6.295 0.110 C:\Python27\lib\site-packages\scipy\integrate\quadpack.py:327(_quad)
57 0.069 0.001 6.295 0.110 {scipy.integrate._quadpack._qagse}
66423 0.154 0.000 6.226 0.000 C:\Users\Will\Documents\Python_scripts\hist_fit.py:41(product_func)
4 0.000 0.000 6.167 1.542 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:2107(_linesearch_powell)
4 0.000 0.000 6.166 1.542 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:1830(brent)
4 0.000 0.000 6.166 1.542 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:1887(_minimize_scalar_brent)
4 0.001 0.000 6.166 1.542 C:\Python27\lib\site-packages\scipy\optimize\optimize.py:1717(optimize)
和scipy方法curve_fit:
pr2 = cProfile.Profile()
pr2.enable()
from scipy.optimize import curve_fit
(mu_curve_fit, sigma_curve_fit), _ = curve_fit(gaussian, kde.support, kde.density, p0=(population_mean, population_std_dev))
pr2.disable()
s = StringIO.StringIO()
sortby = 'cumulative'
ps = pstats.Stats(pr2, stream=s).sort_stats(sortby)
ps.print_stats()
print s.getvalue()
122 function calls in 0.001 seconds
Ordered by: cumulative time
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.001 0.001 C:\Python27\lib\site-packages\scipy\optimize\minpack.py:452(curve_fit)
1 0.000 0.000 0.001 0.001 C:\Python27\lib\site-packages\scipy\optimize\minpack.py:256(leastsq)
1 0.000 0.000 0.001 0.001 {scipy.optimize._minpack._lmdif}
19 0.000 0.000 0.001 0.000 C:\Python27\lib\site-packages\scipy\optimize\minpack.py:444(_general_function)
19 0.000 0.000 0.000 0.000 C:\Users\Will\Documents\Python_scripts\hist_fit.py:29(gaussian)
1 0.000 0.000 0.000 0.000 C:\Python27\lib\site-packages\scipy\linalg\basic.py:314(inv)
1 0.000 0.000 0.000 0.000 C:\Python27\lib\site-packages\scipy\optimize\minpack.py:18(_check_func)
您可以看到 curve_fit 方法快很多，结果：
from numpy import linspace
xs = linspace(-1, 1, num=1000) * sigma_overlap_fit * 6 + mu_overlap_fit
ax.plot(xs, gaussian(xs, mu_overlap_fit, sigma_overlap_fit), color="orange", lw=2)
xs = linspace(-1, 1, num=1000) * sigma_curve_fit * 6 + mu_curve_fit
ax.plot(xs, gaussian(xs, mu_curve_fit, sigma_curve_fit), color="purple", lw=2)
pyplot.show()
非常相似。我会推荐 curve_fit。在这种情况下，速度提高了 6000 倍。当底层数据更复杂时，区别是 一点点 更多，但差别不大，你仍然可以获得巨大的加速。这是一个适合 6 个均匀分布的正态分布的示例：
选择 curve_fit!