31

我有两个变量(x 和 y),它们之间的关系有点像 sigmoid,我需要找到某种预测方程,让我能够在给定 x 的任何值的情况下预测 y 的值。我的预测方程需要显示两个变量之间的某种 S 型关系。因此,我不能满足于产生一条直线的线性回归方程。我需要看到在两个变量的图表的左右两侧出现的斜率逐渐的曲线变化。

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

如果你运行下面的代码,你可以看到它给出了一个向下的抛物线,这不是我的变量之间的关系应该看起来的样子。相反,我的两个变量之间应该有更多的 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 值重新转换为最大 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 是在最大测试条件下来自不同测试参与者的值范围的平均值。


下面的第二次编辑:

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

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

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

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 中的元素数?

4

4 回答 4

42

使用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. 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. 所以sigmoid(p,x)返回一个numpy数组。在numpybook中有更完整的解释它是如何工作的(认真的 numpy 用户需要阅读)。

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

真的。leastsq试图最小化残差(差)的平方和。它搜索参数空间(所有可能的值p)寻找p最小化平方和的值。和发送到x,是您的原始数据值。它们是固定的。他们不会改变。试图最小化的是s(sigmoid 函数中的参数) 。yresidualspleastsq

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

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

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

您可以认为作为第一次通过的最小平方算法(实际上是列文伯格-马夸特算法)的一部分,最小平方调用residuals(p_guess,x,y). 注意两者之间的视觉相似性

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

residuals(p_guess,x,y)

它可以帮助您记住参数的顺序和含义leastsq

residuals, likesigmoid返回一个 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_guessleastsqpleastsq

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

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

于 2010-11-29T22:06:12.367 回答
4

正如上面@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()

结果如下图所示:

在此处输入图像描述

于 2018-09-04T13:49:55.607 回答
2

对于 Python 中的逻辑回归,scikits-learn公开了高性能拟合代码:

http://scikit-learn.sourceforge.net/modules/linear_model.html#logistic-regression

于 2011-04-28T10:02:26.697 回答
2

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

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

于 2010-11-29T21:34:55.907 回答