python - 使用scipy，numpy，python等进行S型回归

我有两个变量(x和y)，它们之间存在某种S型关系，并且我需要找到某种预测方程，使我能够在给定任何x值的情况下预测y的值。我的预测方程式需要显示两个变量之间的某种S形关系。因此，我无法解决产生一条线的线性回归方程。我需要看到在两个变量的曲线图的左右两侧都发生了斜率的逐渐曲线变化。

在谷歌搜索曲线回归和python之后，我开始使用numpy.polyfit，但这给了我可怕的结果，如果您运行下面的代码，您可以看到。 谁能告诉我如何重新编写以下代码以获得所需的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()

编辑如下:(重新构建问题)

您的回应及其速度令人印象深刻。谢谢你
但是，为了产生更有效的结果，我需要重新构造数据值。这意味着将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)的S型曲线沿z轴旋转了180度。同样，1.000并不是x值的最大值。取而代之的是，1.000是最大测试条件下不同测试参与者的值范围的平均值。

下面的第二个编辑:

谢谢你，ubuntu。我仔细阅读了您的代码，并在scipy文档中查找了代码的各个方面。由于您的名字似乎是作为scipy文档的作者而冒出来的，所以我希望您可以回答以下问题:

1.)Minimumsq()是否会调用残差()，然后该残差返回输入y vector 和sigmoid()函数返回的y vector 之间的差？如果是这样，它如何解决输入y vector 和sigmoid()函数返回的y vector 的长度差异？

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

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

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

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

最佳答案

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

产量

具有S型参数

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

请注意，对于较新版本的scipy(例如0.9)，还有scipy.optimize.curve_fit函数，它比leastsq易于使用。有关使用curve_fit拟合Sigmoids的相关讨论，可以找到here。

编辑:添加了resize函数，以便可以重新缩放原始数据并转移到适合任何所需边界框的位置。



免责声明:我不是scipy文档的作者。我只是一个用户，而且还是一个新手。我对leastsq的了解大部分来自阅读Travis Oliphant撰写的this tutorial。



是!究竟。



长度是相同的:

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的妙处之一是，它允许您编写在整个数组上运行的“ vector ”方程。

y = c / (1 + np.exp(-k*(x-x0))) + y0

看起来好像可以在浮点数上工作(实际上是这样)，但是如果将x设置为numpy数组，并将c，k，x0和y0设置为浮点数，则该公式将y定义为与x形状相同的numpy数组。因此，sigmoid(p,x)返回一个numpy数组。在numpybook中有一个更完整的解释(对于numpy的严重用户必读)。



真正。 leastsq尝试最小化残差(差)的平方和。它搜索参数空间(p的所有可能值)，以寻找使该平方和最小的p。发送到x的y和residuals是您的原始数据值。它们是固定的。他们没有改变。 p尝试最小化的是leastsq s(Sigmoid函数中的参数)。



正是如此!像牛顿的方法一样，leastsq需要对p进行初始猜测。您将其作为p_guess提供。当你看到

scipy.optimize.leastsq(residuals,p_guess,args=(x,y))

您可以认为，作为minimumsq算法(实际上是Levenburg-Marquardt算法)的一部分，第一次通过，minimumsq调用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)随之变化。



好吧，不完全是。您现在知道，随着p_guess搜索最小化leastsq的p值，residuals(p,x,y)有所不同。发送到p的p_guess(er，leastsq)具有与p返回的leastsq相同的形状。显然，除非您是个猜测家，否则值应该是不同的:)



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