sigmoidal回归与scipy,numpy,python等

Question

我有两个变量(x和y)彼此之间有一些S形关系,我需要找到某种预测方程,这将使我能够在给定任何x值的情况下预测y的值.我的预测方程需要显示两个变量之间的某种S形关系.因此,我不能满足于产生线的线性回归方程.我需要看到两个变量图的右侧和左侧出现的斜率的逐渐曲线变化.

我在googling曲线回归和python之后开始使用numpy.polyfit,但是如果你运行下面的代码,这给了我可怕的结果. 任何人都可以告诉我如何重新编写下面的代码,以获得我想要的S形回归方程式吗？

如果你运行下面的代码,你可以看到它给出了一个向下的抛物线,这不是我的变量之间的关系应该是什么样子.相反,我的两个变量之间应该存在更多的S形关系,但是与我在下面的代码中使用的数据紧密相符.下面代码中的数据来自大样本研究的手段,因此它们包含的统计功效比五个数据点所暗示的要多.我没有大样本研究的实际数据,但我确实有下面的方法和他们的标准偏差(我没有显示).我更愿意用下面列出的平均数据绘制一个简单的函数,但如果复杂性会带来实质性的改进,代码可能会变得更加复杂.

如何更改我的代码以显示最适合的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()

---

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

您的反应及其速度令人印象深刻.谢谢你,unutbu.但是,为了产生更有效的结果,我需要重新构建我的数据值.这意味着将x值重新转换为max 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()

你能告诉我如何解决这个修改后的代码吗？
注意:通过重新投射数据,我基本上围绕z轴旋转了2d(x,y)sigmoid 180度.此外,1.000实际上不是x值的最大值.相反,1.000是最大测试条件下来自不同测试参与者的值范围的平均值.

---

下面的第二个编辑:

谢谢你,ubuntu.我仔细阅读了你的代码,并在scipy文档中查看了它的各个方面.由于您的名字似乎弹出作为scipy文档的作者,我希望您可以回答以下问题:

1.)leastsq()是否调用residuals(),然后返回输入y-vector和sigmoid()函数返回的y-vector之间的差异？如果是这样,它如何解释输入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中的元素数量？

Answer 1

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

产量

用sigmoid参数

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

请注意,对于较新版本的scipy(例如0.9),还有scipy.optimize.curve_fit函数,它比使用起来更容易leastsq.有关拟合sigmoids使用的相关讨论curve_fit可以在这里找到.

编辑:resize添加了一个函数,以便可以重新缩放原始数据并移动以适合任何所需的边界框.

"你的名字似乎弹出作为scipy文档的作者"

免责声明:我不是scipy文档的作者.我只是一个用户,也是一个新手.我所了解的大部分内容leastsq来自阅读本教程,由Travis Oliphant撰写.

1.)leastsq()是否调用residuals(),然后返回输入y-vector和sigmoid()函数返回的y-vector之间的差异？

是! 究竟.

如果是这样,它如何解释输入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是相同形状的numpy的阵列x.所以sigmoid(p,x)返回一个numpy数组.有关如何在numpybook中工作的更完整的解释(严格的numpy用户需要阅读).

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

真正.leastsq试图最小化残差平方和(差异).它搜索参数空间(所有可能的值p),寻找p最小化该平方和的方法.在x与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))

你可以认为作为最小化算法(实际上是Levenburg-Marquardt算法)的一部分,作为第一遍,最小调用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相同,而leastsq()函数在返回之前在内部使用该p？

好吧,不完全是.正如您现在所知,p_guess随着leastsq对p最小化值的搜索而变化residuals(p,x,y).发送到的p(er,p_guess)leastsq具有与p返回的形状相同的形状leastsq.显然,值应该是不同的,除非你是一个猜测者的地狱:)

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

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

Answer 2

我认为任何阶次的多项式拟合都不会得到好的结果——因为对于足够大和足够小的 X，所有多项式都会趋于无穷大，但 sigmoid 曲线将在每个方向上渐近地接近某个有限值。

我不是Python程序员，所以我不知道numpy是否有更通用的曲线拟合例程。如果您必须自己动手，也许这篇关于逻辑回归的文章会给您一些想法。