2010-06-22 86 views
0

我需要用直线拟合不同数据集中的一些点。从每个数据集中,我想要适合一条线。所以我得到了描述i线的参数ai和bi:ai + bi * x。问题是我想强制每个AI都是平等的,因为我需要相同的intercepta。我在这里找到了一个教程:http://www.scipy.org/Cookbook/FittingData#head-a44b49d57cf0165300f765e8f1b011876776502f。不同的是,我不知道我有多少数据集。我的代码是这样的:scipy智能优化

from numpy import * 
from scipy import optimize 

# here I have 3 dataset, but in general I don't know how many dataset are they 
ypoints = [array([0, 2.1, 2.4]), # first dataset, 3 points 
      array([0.1, 2.1, 2.9]), # second dataset 
      array([-0.1, 1.4])]  # only 2 points 

xpoints = [array([0, 2, 2.5]),  # first dataset 
      array([0, 2, 3]),  # second, also x coordinates are different 
      array([0, 1.5])]   # the first coordinate is always 0 

fitfunc = lambda a, b, x: a + b * x 
errfunc = lambda p, xs, ys: array([ yi - fitfunc(p[0], p[i+1], xi) 
            for i, (xi,yi) in enumerate(zip(xs, ys)) ]) 


p_arrays = [r_[0.]] * len(xpoints) 
pinit = r_[[ypoints[0][0]] + p_arrays] 
fit_parameters, success = optimize.leastsq(errfunc, pinit, args = (xpoints, ypoints)) 

Traceback (most recent call last): 
    File "prova.py", line 19, in <module> 
    fit_parameters, success = optimize.leastsq(errfunc, pinit, args = (xpoints, ypoints)) 
    File "/usr/lib64/python2.6/site-packages/scipy/optimize/minpack.py", line 266, in leastsq 
    m = check_func(func,x0,args,n)[0] 
    File "/usr/lib64/python2.6/site-packages/scipy/optimize/minpack.py", line 12, in check_func 
    res = atleast_1d(thefunc(*((x0[:numinputs],)+args))) 
    File "prova.py", line 14, in <lambda> 
    for i, (xi,yi) in enumerate(zip(xs, ys)) ]) 
ValueError: setting an array element with a sequence. 

回答

1

如果你只是需要一个线性拟合,那么最好用线性回归而不是非线性优化器来估计它。 更适合的统计数据可以通过使用scikits.statsmodels来获得。

import numpy as np 
from numpy import array 

ypoints = np.r_[array([0, 2.1, 2.4]), # first dataset, 3 points 
      array([0.1, 2.1, 2.9]), # second dataset 
      array([-0.1, 1.4])]  # only 2 points 

xpoints = [array([0, 2, 2.5]),  # first dataset 
      array([0, 2, 3]),  # second, also x coordinates are different 
      array([0, 1.5])]   # the first coordinate is always 0 

xp = np.hstack(xpoints) 
indicator = [] 
for i,a in enumerate(xpoints): 
    indicator.extend([i]*len(a)) 

indicator = np.array(indicator) 


x = xp[:,None]*(indicator[:,None]==np.arange(3)).astype(int) 
x = np.hstack((np.ones((xp.shape[0],1)),x)) 

print np.dot(np.linalg.pinv(x), ypoints) 
# [ 0.01947973 0.98656987 0.98481549 0.92034684] 

回归量的基体有一个共同的截距,但不同列的每个数据集:

>>> x 
array([[ 1. , 0. , 0. , 0. ], 
     [ 1. , 2. , 0. , 0. ], 
     [ 1. , 2.5, 0. , 0. ], 
     [ 1. , 0. , 0. , 0. ], 
     [ 1. , 0. , 2. , 0. ], 
     [ 1. , 0. , 3. , 0. ], 
     [ 1. , 0. , 0. , 0. ], 
     [ 1. , 0. , 0. , 1.5]]) 
+0

谢谢,我喜欢它。问题是现在我需要使用错误,我的意思是:y点有错误,我需要用1 /错误^ 2加权。我怎样才能用你的代码做到这一点? – 2010-06-28 18:31:38

+0

最好的方法是使用scikits.statsmodels,因为在这种情况下,所有的预测,预测和结果统计都是预制的。 http://pypi.python.org/pypi/scikits.statsmodels/0.2.0和链接 获得预测的y PARAMS = np.dot(np.linalg.pinv(x)中,ypoints中) ypred = np.dot(x,params) errors = ypoints - ypred ... 如果您是指称误差,使用加权最小二乘法,那么x和y点都需要除以误差标准差,或者使用scikits.statsmodels中的WLS类。 – user333700 2010-06-28 22:15:58

1

(旁注:使用def,不lambda分配给一个名字 - 这是完全愚蠢和有什么,但不足之处,lambda唯一使用正在制作匿名函数!)。

errfunc应该返回浮点数字序列(数组或其他方式),但它不是,因为你试图把作为阵列的项目,它们分别y点(记住不同的阵列,ypoints又名ys是数组列表!)和拟合函数的结果。因此,您需要将表达式yi - fitfunc(p[0], p[i+1], xi)“折叠”为单个浮点数,例如, norm(yi - fitfunc(p[0], p[i+1], xi))

+0

我无法理解在链路 – 2010-06-22 20:37:40

+0

@Wiso的例子之间的差异,该代码是相当晦涩难懂(带有未记录的'r_'等),但我确信它最终会从errfunc返回一个浮点数序列,因为这就是文档需要的 - 使用调试器或使用'print'对其进行检查验证一下,如果你可以自己成功运行这个代码(我不能 - 什么是'r_'?!)。 – 2010-06-22 20:47:31

+0

好的,我解决了。通常scipy是有据可查的,'r_'函数在这里:http://docs.scipy.org/doc/numpy/reference/generated/numpy.r_.html 解决的办法是将'array'更改为'concatenate'像你说的那样获得一个唯一的数组。特别是我得到的系数:'[0.01947973 0.98656987 0.98481549 0.92034684]',他们似乎正在纠正 – 2010-06-22 21:16:16