使用scipy、numpy等进行sigmoidal回归
我有两个变量（x和y），它们之间的关系有点像S形曲线。我需要找到一个预测方程，这样我就可以根据x的值来预测y的值。这个预测方程需要能体现这两个变量之间的S形关系。因此，我不能仅仅使用一个线性回归方程，因为那样只会画出一条直线。我需要看到在图的左右两边，斜率逐渐变化的曲线。
我开始使用numpy.polyfit，搜索了曲线回归和python，但结果很糟糕。如果你运行下面的代码，就会看到这些糟糕的结果。有没有人能告诉我如何修改下面的代码，以得到我想要的S形回归方程？
如果你运行下面的代码，你会发现它给出的结果是一个向下的抛物线，这并不是我变量之间应该有的关系。实际上，我的两个变量之间应该有更明显的S形关系，而且要与我在下面代码中使用的数据紧密贴合。下面代码中的数据是来自一个大样本研究的均值，所以它们的统计能力比五个数据点看起来要强。我没有这个大样本研究的实际数据，但我有下面的均值和它们的标准差（我没有显示标准差）。我更希望能用下面的均值数据画一个简单的函数，但如果复杂的代码能带来显著的改进，那也可以。
我该如何修改我的代码，以显示一个最佳拟合的S形函数，最好是使用scipy、numpy和python？ 这是我当前的代码版本，需要修正：
import numpy as np
import matplotlib.pyplot as plt
# Create numpy data arrays
x = np.array([821,576,473,377,326])
y = np.array([255,235,208,166,157])
# Use polyfit and poly1d to create the regression equation
z = np.polyfit(x, y, 3)
p = np.poly1d(z)
xp = np.linspace(100, 1600, 1500)
pxp=p(xp)
# Plot the results
plt.plot(x, y, '.', xp, pxp, '-')
plt.ylim(140,310)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
下面的编辑：（重新表述了问题）
你的回复和速度都让我很惊讶。谢谢你，unutbu。不过，为了得到更有效的结果，我需要重新整理我的数据值。这意味着要把x值转换为最大x值的百分比，同时把y值转换为原始数据中x值的百分比。我尝试用你的代码做到这一点，得到了以下结果：
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
# Create numpy data arrays
'''
# Comment out original data
#x = np.array([821,576,473,377,326])
#y = np.array([255,235,208,166,157])
'''
# Re-calculate x values as a percentage of the first (maximum)
# original x value above
x = np.array([1.000,0.702,0.576,0.459,0.397])
# Recalculate y values as a percentage of their respective x values
# from original data above
y = np.array([0.311,0.408,0.440,0.440,0.482])
def sigmoid(p,x):
x0,y0,c,k=p
y = c / (1 + np.exp(-k*(x-x0))) + y0
return y
def residuals(p,x,y):
return y - sigmoid(p,x)
p_guess=(600,200,100,0.01)
(p,
cov,
infodict,
mesg,
ier)=scipy.optimize.leastsq(residuals,p_guess,args=(x,y),full_output=1,warning=True)
'''
# comment out original xp to allow for better scaling of
# new values
#xp = np.linspace(100, 1600, 1500)
'''
xp = np.linspace(0, 1.1, 1100)
pxp=sigmoid(p,xp)
x0,y0,c,k=p
print('''\
x0 = {x0}
y0 = {y0}
c = {c}
k = {k}
'''.format(x0=x0,y0=y0,c=c,k=k))
# Plot the results
plt.plot(x, y, '.', xp, pxp, '-')
plt.ylim(0,1)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
你能告诉我如何修正这个修改后的代码吗？
注意：通过重新整理数据，我实际上是将二维（x,y）S形图绕z轴旋转了180度。此外，1.000并不是真正的x值的最大值，而是不同测试参与者在最大测试条件下的值范围的均值。
第二次编辑：
谢谢你，ubuntu。我仔细阅读了你的代码，并查阅了scipy文档中的相关内容。因为你的名字似乎是scipy文档的作者之一，我希望你能回答以下问题：
1.) leastsq()是否调用了residuals()，然后返回输入的y向量与sigmoid()函数返回的y向量之间的差异？如果是这样，它是如何处理输入y向量和sigmoid()函数返回的y向量长度不同的问题的？
2.) 看起来我可以对任何数学方程调用leastsq()，只要我通过一个residuals函数来访问那个数学方程，而这个residuals函数又调用了数学函数。这是真的吗？
3.) 我注意到p_guess的元素数量与p相同。这是否意味着p_guess的四个元素分别对应于x0、y0、c和k返回的值？
4.) 作为参数传递给residuals()和sigmoid()函数的p与leastsq()输出的p是同一个p吗？leastsq()函数在返回之前是否在内部使用了这个p？
5.) p和p_guess可以有任意数量的元素吗？只要p的元素数量等于p_guess的元素数量，就可以根据所用模型的复杂性来决定元素数量吗？
4 个回答
我觉得用多项式来拟合数据，不管是几次多项式，结果都不会太好。因为所有的多项式在X值非常大或非常小时，结果都会无限增大或减小。而S型曲线在两个方向上会逐渐接近某个有限的值，不会像多项式那样无限。
我不是Python程序员，所以不太清楚numpy有没有更通用的曲线拟合方法。如果你需要自己动手做，或许可以看看这篇关于逻辑回归的文章，里面可能会给你一些灵感。
正如@unutbu上面提到的，scipy现在提供了一个叫做 scipy.optimize.curve_fit 的功能，它的使用方法更简单。如果有人想快速了解这个过程是怎样的，我在下面提供了一个简单的例子：
from scipy.optimize import curve_fit
def sigmoid(x, k, x0):
return 1.0 / (1 + np.exp(-k * (x - x0)))
# Parameters of the true function
n_samples = 1000
true_x0 = 15
true_k = 1.5
sigma = 0.2
# Build the true function and add some noise
x = np.linspace(0, 30, num=n_samples)
y = sigmoid(x, k=true_k, x0=true_x0)
y_with_noise = y + sigma * np.random.randn(n_samples)
# Sample the data from the real function (this will be your data)
some_points = np.random.choice(1000, size=30) # take 30 data points
xdata = x[some_points]
ydata = y_with_noise[some_points]
# Fit the curve
popt, pcov = curve_fit(sigmoid, xdata, ydata)
estimated_k, estimated_x0 = popt
# Plot the fitted curve
y_fitted = sigmoid(x, k=estimated_k, x0=estimated_x0)
# Plot everything for illustration
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, y_fitted, '--', label='fitted')
ax.plot(x, y, '-', label='true')
ax.plot(xdata, ydata, 'o', label='samples')
ax.legend()
这个例子的结果在下图中展示：
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
def sigmoid(p,x):
x0,y0,c,k=p
y = c / (1 + np.exp(-k*(x-x0))) + y0
return y
def residuals(p,x,y):
return y - sigmoid(p,x)
def resize(arr,lower=0.0,upper=1.0):
arr=arr.copy()
if lower>upper: lower,upper=upper,lower
arr -= arr.min()
arr *= (upper-lower)/arr.max()
arr += lower
return arr
# raw data
x = np.array([821,576,473,377,326],dtype='float')
y = np.array([255,235,208,166,157],dtype='float')
x=resize(-x,lower=0.3)
y=resize(y,lower=0.3)
print(x)
print(y)
p_guess=(np.median(x),np.median(y),1.0,1.0)
p, cov, infodict, mesg, ier = scipy.optimize.leastsq(
residuals,p_guess,args=(x,y),full_output=1,warning=True)
x0,y0,c,k=p
print('''\
x0 = {x0}
y0 = {y0}
c = {c}
k = {k}
'''.format(x0=x0,y0=y0,c=c,k=k))
xp = np.linspace(0, 1.1, 1500)
pxp=sigmoid(p,xp)
# Plot the results
plt.plot(x, y, '.', xp, pxp, '-')
plt.xlabel('x')
plt.ylabel('y',rotation='horizontal')
plt.grid(True)
plt.show()
得到的结果是
这里用的是sigmoid（S形曲线）参数
x0 = 0.826964424481
y0 = 0.151506745435
c = 0.848564826467
k = -9.54442292022
需要注意的是，对于更新版本的scipy（比如0.9），还有一个更简单的函数 scipy.optimize.curve_fit，比 leastsq 更容易使用。关于使用 curve_fit 拟合sigmoid的相关讨论可以在 这里找到。
补充：添加了一个 resize 函数，这样原始数据可以重新缩放和移动，以适应任何想要的边界框。
"你的名字似乎经常出现在scipy文档的作者名单中"
声明：我并不是scipy文档的作者。我只是一个用户，而且还是个新手。我对 leastsq 的了解大部分来自阅读 这个教程，是Travis Oliphant写的。
1.) leastsq() 是不是调用了 residuals()，然后返回输入的y向量和sigmoid()函数返回的y向量之间的差值？
没错！就是这样。
如果是这样，那它是如何处理输入y向量和sigmoid()函数返回的y向量长度不同的问题的？
长度是相同的：
In [138]: x
Out[138]: array([821, 576, 473, 377, 326])
In [139]: y
Out[139]: array([255, 235, 208, 166, 157])
In [140]: p=(600,200,100,0.01)
In [141]: sigmoid(p,x)
Out[141]:
array([ 290.11439268, 244.02863507, 221.92572521, 209.7088641 ,
206.06539033])
Numpy的一个奇妙之处在于，它允许你编写“向量”方程，能够对整个数组进行操作。
y = c / (1 + np.exp(-k*(x-x0))) + y0
看起来像是对浮点数进行操作（确实是这样），但如果你把 x 设为一个numpy数组，而 c、k、x0、y0 设为浮点数，那么这个方程就定义了 y 为与 x 形状相同的numpy数组。所以 sigmoid(p,x) 返回的是一个numpy数组。关于这个是如何工作的，有更详细的解释在 numpybook 中（这是认真使用numpy的人必读的书）。
2.) 看起来我可以对任何数学方程调用leastsq()，只要我通过一个residuals函数来访问这个数学方程，而这个函数又调用了数学函数。这是真的吗？
是真的。leastsq 试图最小化残差（差值）的平方和。它在参数空间（所有可能的 p 值）中搜索，寻找能够最小化这个平方和的 p。发送给 residuals 的 x 和 y 是你的原始数据值。它们是固定的，不会改变。是 p（sigmoid函数中的参数）被 leastsq 尝试最小化。
3.) 我还注意到 p_guess 的元素数量和 p 是一样的。这是否意味着 p_guess 的四个元素分别对应于 x0、y0、c 和 k 的返回值？
完全正确！就像牛顿法一样，leastsq 需要一个 p 的初始猜测。你通过 p_guess 提供这个初始值。当你看到
scipy.optimize.leastsq(residuals,p_guess,args=(x,y))
你可以把它想象成leastsq算法的一部分（实际上是Levenburg-Marquardt算法），作为第一次尝试，leastsq调用了 residuals(p_guess,x,y)。
注意到
(residuals,p_guess,args=(x,y))
和
residuals(p_guess,x,y)
之间的视觉相似性，这可能会帮助你记住 leastsq 的参数顺序和含义。
residuals 和 sigmoid 一样，返回的是一个numpy数组。数组中的值会被平方，然后求和。这就是需要优化的数值。然后 p_guess 会被改变，leastsq 会寻找一组值，使得 residuals(p_guess,x,y) 最小。
4.) 作为参数传递给 residuals() 和 sigmoid() 函数的 p 和最后由 leastsq() 输出的 p 是同一个 p 吗？leastsq() 函数在返回之前内部使用的是这个 p 吗？
其实不完全是。你现在应该知道，p_guess 在 leastsq 搜索最小化 residuals(p,x,y) 的过程中会被改变。发送给 leastsq 的 p（也就是 p_guess）和 leastsq 返回的 p 形状是一样的。显然，除非你是个非常厉害的猜测者，否则它们的值应该是不同的 :)
5.) p 和 p_guess 可以有任意数量的元素吗？只要 p 的元素数量和 p_guess 的元素数量相等，取决于所用模型方程的复杂性？
可以的。我还没有对 leastsq 进行过大量参数的压力测试，但它是一个非常强大的工具。