使用 scipy、numpy、python 等进行 sigmoidal 回归

我有两个变量（x 和 y），它们彼此之间存在某种 S 型关系，并且我需要找到某种预测方程，使我能够在给定 x 的任何值的情况下预测 y 的值。我的预测方程需要显示两个变量之间的 S 形关系。因此，我不能满足于产生一条直线的线性回归方程。我需要看到两个变量图表右侧和左侧发生的斜率逐渐曲线变化。

在谷歌搜索曲线回归和 python 之后，我开始使用 numpy.polyfit，但这给了我可怕的结果，如果你运行下面的代码，你可以看到。谁能告诉我如何重写下面的代码以获得我想要的 sigmoidal 回归方程类型？

如果运行下面的代码，您可以看到它给出了一条向下的抛物线，这不是我的变量之间的关系应有的样子。相反，我的两个变量之间应该有更多的 sigmoidal 关系，但与我在下面的代码中使用的数据紧密配合。下面代码中的数据来自大样本研究，因此它们比五个数据点可能暗示的统计能力更强。我没有大样本研究的实际数据，但我有下面的方法及其标准差（我没有显示）。我更愿意只用下面列出的平均数据绘制一个简单的函数，但如果复杂性能够带来实质性的改进，代码可能会变得更复杂。

如何更改代码以显示 sigmoidal 函数的最佳拟合，最好使用 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()

编辑如下：（重新提出问题）

您的回复及其速度令人印象深刻。谢谢你，乌努特布。 但是，为了产生更有效的结果，我需要重新构建我的数据值。这意味着将 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()

你能告诉我如何修复这个修改后的代码吗？
注意：通过重新转换数据，我基本上将 2d (x,y) sigmoid 绕 z 轴旋转了 180 度。此外，1.000 并不是真正的 x 值的最大值。相反，1.000 是最大测试条件下不同测试参与者的值范围的平均值。

下面第二次编辑：

谢谢你，乌班图。我仔细阅读了你的代码，并在 scipy 文档中查找了它的各个方面。由于您的名字似乎是 scipy 文档的作者，我希望您能回答以下问题：

1.）leastsq()是否调用residuals()，然后返回输入y向量与sigmoid()函数返回的y向量之间的差？如果是这样，它如何解释输入 y 向量和 sigmoid() 函数返回的 y 向量的长度差异？

2.) 看起来我可以为任何数学方程调用leastsq()，只要我通过残差函数访问该数学方程，该函数又调用数学函数。这是真的？

3.) 另外，我注意到 p_guess 与 p 具有相同数量的元素。这是否意味着 p_guess 的四个元素分别与 x0、y0、c 和 k 返回的值按顺序对应？

4.) 作为参数发送到residuals() 和sigmoid() 函数的p 与leastsq() 输出的p 是否相同，并且leastsq() 函数在返回之前在内部使用该p？

5.) p 和 p_guess 是否可以有任意数量的元素，具体取决于用作模型的方程的复杂性，只要 p 中的元素数量等于 p_guess 中的元素数量？

Using scipy.optimize.leastsq http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html#scipy.optimize.leastsq:

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()

yields

带 sigmoid 参数

x0 = 0.826964424481
y0 = 0.151506745435
c = 0.848564826467
k = -9.54442292022

请注意，对于较新版本的 scipy（例如 0.9），还有scipy.optimize.curve_fit http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy-optimize-curve-fit比这个更容易使用的功能leastsq。使用拟合 sigmoid 的相关讨论curve_fit可以被找寻到here http://comments.gmane.org/gmane.comp.python.scientific.user/26237.

Edit: A resize添加了函数，以便可以重新缩放和移动原始数据以适应任何所需的边界框。

“你的名字似乎以作家的身份出现 scipy 文档”

免责声明：我不是 scipy 文档的作者。我只是一个用户，而且还是个新手。我所知道的很多事情leastsq来自阅读本教程 http://www.tau.ac.il/~kineret/amit/scipy_tutorial/，特拉维斯·奥利芬特撰写。

1.）leastsq()是否调用residuals()，然后返回差值 输入 y 向量和 sigmoid() 返回的 y 向量 功能？

是的！确切地。

如果是这样，它是如何解释的 输入长度的差异 y 向量和返回的 y 向量 sigmoid() 函数？

长度相同：

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是一个形状相同的 numpy 数组x. So sigmoid(p,x)返回一个 numpy 数组。关于它如何工作的更完整的解释在numpybook http://web.mit.edu/dvp/Public/numpybook.pdf（numpy 认真用户必读）。

2.）看起来我可以为任何数学方程调用leastsq()，只要我 通过访问该数学方程 残差函数，进而 调用数学函数。这是真的？

True. leastsq尝试最小化残差（差值）的平方和。它搜索参数空间（所有可能的值p）寻找p最小化平方和。这x and y发给residuals，是您的原始数据值。它们是固定的。他们没有改变。这是ps（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))

您可以认为，作为第一遍的 lesssq 算法（实际上是 Levenburg-Marquardt 算法）的一部分，leastsq 调用residuals(p_guess,x,y)。 注意之间的视觉相似性

(residuals,p_guess,args=(x,y))

and

residuals(p_guess,x,y)

它可以帮助您记住论点的顺序和含义leastsq.

residuals, like sigmoid返回一个 numpy 数组。对数组中的值进行平方，然后求和。这是要击败的数字。p_guess然后变化​​为leastsq寻找一组最小化的值residuals(p_guess,x,y).

4.) 是作为参数发送给residuals() 的p 以及 sigmoid() 的功能与 p 相同 将由 lesssq() 输出，并且 lesssq() 函数正在使用 p 返回之前内部？

嗯，不完全是。正如你现在所知，p_guess变化为leastsq搜索p最小化的值residuals(p,x,y). The p (er, p_guess) 发送到leastsq具有相同的形状p由返回leastsq。显然，这些值应该是不同的，除非你是一个猜测者:)

5.) p 和 p_guess 可以有任意数量的元素，具体取决于 所使用方程的复杂性 作为一个模型，只要数量 p 中的元素等于数量 p_guess 中有多少个元素？

是的。我没有进行压力测试leastsq对于非常大量的参数，但它是一个非常强大的工具。