我有两个变量(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中的元素数量?
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
拟合sigmoids的相关讨论可以找到here。
编辑:添加了resize
函数,以便可以重新缩放原始数据并移动以适合任何所需的边界框。
“你的名字似乎弹出作为scipy文档的作者”
免责声明:我不是scipy文档的作者。我只是一个用户,也是一个新手。我对leastsq
的了解大部分来自阅读由Travis Oliphant撰写的this tutorial。
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
定义为与x
相同形状的numpy数组。所以sigmoid(p,x)
返回一个numpy数组。有一个更完整的解释,这是如何在numpybook(严重的numpy用户必读)。
2.)看起来我可以为任何数学方程式调用leastsq(),只要我通过残差函数访问该数学方程式,而残差函数又调用数学函数。这是真的?
真正。 leastsq
试图最小化残差的平方和(差异)。它搜索参数空间(p
的所有可能值),寻找最小化该平方和的p
。发送给x
的y
和residuals
是您的原始数据值。他们是固定的。他们不会改变。它是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))
你可以认为作为最小化算法(实际上是Levenberg-Marquardt算法)的一部分作为第一遍,最小的调用qazxsw i。注意之间的视觉相似性
residuals(p_guess,x,y)
和
(residuals,p_guess,args=(x,y))
它可以帮助你记住residuals(p_guess,x,y)
的参数的顺序和含义。
leastsq
,像residuals
一样返回一个numpy数组。数组中的值是平方的,然后求和。这是要击败的数字。然后sigmoid
变化,因为p_guess
寻找一组最小化leastsq
的值。
4.)作为参数发送到residuals()和sigmoid()的p是否与将由leastsq()输出的p相同,而leastsq()函数在返回之前在内部使用该p?
好吧,不完全是。正如你现在所知,residuals(p_guess,x,y)
是多种多样的,因为p_guess
搜索最小化leastsq
的p
值。发往residuals(p,x,y)
的p
(呃,p_guess
)与leastsq
返回的p
具有相同的形状。显然,值应该是不同的,除非你是一个猜测者的地狱:)
5.)p和p_guess可以有任意数量的元素,这取决于用作模型的方程的复杂性,只要p中的元素数等于p_guess中的元素数量?
是。我没有对非常大量的参数进行压力测试leastsq
,但它是一个非常强大的工具。
我不认为你会得到任何程度的多项式拟合得到好的结果 - 因为所有多项式对于足够大和小的X都会变为无穷大,但是S形曲线将渐近地在每个方向上接近某个有限值。
我不是Python程序员,所以我不知道numpy是否有更通用的曲线拟合程序。如果你必须自己动手,也许这篇关于leastsq
的文章会给你一些想法。
对于Python中的逻辑回归,Logistic regression公开了高性能拟合代码:
正如@unutbu所指出的那样,http://scikit-learn.sourceforge.net/modules/linear_model.html#logistic-regression现在提供的scipy
拥有一个不那么复杂的呼叫。如果有人想要快速了解相同过程在这些术语中的样子,我将在下面给出一个最小的例子:
scipy.optimize.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(return_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()