我试图弄清楚使用lmfit包的截距估计如何根据起始值产生不同的结果。我正在使用statsmodels将估算结果与标准OLS估算的结果进行比较。
有人可以协助吗?

码:

from lmfit import minimize, Parameters
import numpy as np
import statsmodels.api as sm

x = np.linspace(0, 15, 10)
x_ols = sm.add_constant(x)
y = range(0,10)

model = sm.OLS(y,x_ols)
results = model.fit()
print "OLS: ", format(results.params[0], '.10f'), format(results.params[1], '.10f')


# define objective function: returns the array to be minimized
def fcn2min(params, x, data):
    a = params['a'].value
    b = params['b'].value

    model = a + b * x
    return model - data

for i in range(-2,3):
    # create a set of Parameters
    params = Parameters()
    params.add('a', value= i)
    params.add('b', value= 20)

    # do fit, here with leastsq model
    result = minimize(fcn2min, params, args=(x, y))
    # print "lmfit: ",result.values # older usage
    print "lmfit: ",result.params.values # newer syntax

最佳答案

如果查看结果,则会发现该常数接近零,近似于scipy.optimize.leastsq的默认容差。

OLS:  -0.0000000000 0.6000000000
lmfit:  {'a': 1.1967327242889913e-08, 'b': 0.59999999932429748}
lmfit:  {'a': 3.2643846243929039e-08, 'b': 0.59999999717671448}
lmfit:  {'a': 3.2644232427278463e-08, 'b': 0.59999999717679642}
lmfit:  {'a': 1.636450187584131e-08, 'b': 0.59999999900662315}
lmfit:  {'a': 3.3578611617157454e-08, 'b': 0.59999999693286854}


如果我使所需的公差更严格,则常数的估计值将接近于零,并且两个估计的系数都将更接近OLS解。

xtol添加到您的代码中之后

result = minimize(fcn2min, params, args=(x, y), xtol=1e-12)


我得到:

OLS:  -0.0000000000 0.6000000000
lmfit:  {'a': -3.5437341988915725e-32, 'b': 0.59999999999999998}
lmfit:  {'a': -1.8490806236697864e-32, 'b': 0.59999999999999998}
lmfit:  {'a': -9.8614500814838325e-32, 'b': 0.59999999999999998}
lmfit:  {'a': 8.328833474508222e-09, 'b': 0.59999999921839664}
lmfit:  {'a': -5.8547746020524404e-32, 'b': 0.59999999999999998}


OLS使用线性代数获得显式解(基于广义逆,pinv)。因此,数值精度很高。

非线性优化使用迭代解决方案,并在达到所需的公差范围内时停止。

在示例中仍然有一种情况仅在8e-9时才是正确的,但这也可能是由于对梯度的数值近似以及其他影响收敛性检查的因素。

10-07 19:11
查看更多