用scipy曲线拟合的多变量回归：总是以一个系统值为准

Question

我一直在用python做多元线性回归（方程的形式是：y=b0+b1*x1+b2*x2+…+bnxn）。我可以成功地解决以下函数：

def MultipleRegressionFunc(X,B):
    y=B[0]
    for i in range(0,len(X)): y += X[i]*B[i+1]
    return y

我现在跳过这个函数的细节。我只想说，使用curve_fit中的scipy包装器和这个函数，已经成功地使我能够求解具有许多变量的系统。你知道吗

现在我想考虑变量之间可能的相互作用，所以我修改了函数如下：

def MultipleRegressionFuncBIS(X,B):
    #Define terms of the equation
    #The first term is 1*b0 (intercept)
    terms=[1]
    #Adding terms for the "non-interaction" part of the equation
    for x in X: terms.append(x)
    #Adding terms for the 'interaction' part of the equations
    for x in list(combinations(X, 2)): terms.append(x[0]*x[1])
    #I'm proceeding in this way because I found that some operations on iterables are not well handled when curve_fit passes numpy arrays to the function
    #Setting a float object to hold the result of the calculation
    y = 0.0
    #Iterating through each term in the equation, and adding the value to y
    for i in range(0, len(terms)):  y += B[i]*terms[i]
    return y

我为上面的代码创建了一个包装函数，以便能够通过curve_fit将多个线性系数传递给它。你知道吗

def wrapper_func(X,*B):
    return MultipleRegressionFuncBIS(X,B)

下面是一些模拟输入，通过应用以下公式生成：1+2*x1+3*x2+4*x3+5*x1*x2+6*x1*x3+7*x2*x3

x1=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53]
x2=[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54]
x3=[3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55]
y=[91, 192, 329, 502, 711, 956, 1237, 1554, 1907, 2296, 2721, 3182, 3679, 4212, 4781, 5386, 6027, 6704, 7417, 8166, 8951, 9772, 10629, 11522, 12451, 13416, 14417, 15454, 16527, 17636, 18781, 19962, 21179, 22432, 23721, 25046, 26407, 27804, 29237, 30706, 32211, 33752, 35329, 36942, 38591, 40276, 41997, 43754, 45547, 47376, 49241, 51142, 53079]

然后通过调用以下公式得到线性系数：

linear_coeffs=list(curve_fit(wrapper_func,[x1,x2,x3],y,p0=[1.1,2.2,3.1,4.1,5.1,6.1,7.1],bounds=(0.0,10.0))[0])
print linear_coeffs

注意，这里p0估计值被手动设置为非常接近实际值的值，以排除curve_fit难以收敛的可能性。你知道吗

然而，这种特殊情况下的输出偏离实际值的程度超出了我的预期（预期值：[1.0,2.0,3.0,4.0,5.0,6.0,7.0]）：

[1.1020684140370627, 2.1149407566785214, 2.9872182044259676, 3.9734017072175436, 5.0575156518729969, 5.9605293645760549, 6.9819549835509491]

现在，我的问题来了。虽然系数与输入模型并不完全匹配，但这是第二个问题。我确实期望在实际例子中出现一些错误，尽管这在这个没有噪音的模拟例子中是令人费解的。我的主要问题是错误是系统性的。在上例中，使用由curve_fit估计的系数，对于x1,x2,x3的所有值，残差系统地等于0.10206841。其他模拟数据集产生不同的，但仍然是系统的残差。你知道吗

你能想出什么解释来解释这个系统性错误吗？你知道吗

我在这里发帖是因为我怀疑这是一个编码问题，而不是一个统计问题。我非常愿意把这个问题转移到交叉验证，如果它证明我犯了统计错误。你知道吗

非常感谢！你知道吗