Python-最小化卡方

我一直在尝试通过最小化卡方使线性模型适合一组应力/应变数据.不幸的是,使用下面的代码不能正确最小化chisqfunc函数.它在初始条件x0处找到最小值,这是不正确的.我浏览了scipy.optimize文档,并进行了测试,以尽量减少其他已正常运行的功能.您能否建议下面的代码修复方法,或者建议我可以使用另一种方法通过最小化卡方来使线性模型适合数据？

import numpy
import scipy.optimize as opt

filename = 'data.csv'

data = numpy.loadtxt(open(filename,"r"),delimiter=",")

stress = data[:,0]
strain = data[:,1]
err_stress = data[:,2]

def chisqfunc((a, b)):
    model = a + b*strain
    chisq = numpy.sum(((stress - model)/err_stress)**2)
    return chisq

x0 = numpy.array([0,0])

result =  opt.minimize(chisqfunc, x0)
print result

感谢您阅读我的问题,我们将不胜感激.

干杯,威尔

编辑：我当前正在使用的数据集：Link to data

解决方法:

问题是您的最初猜测与实际解决方案相距甚远.如果在chisqfunc()中添加一条print语句(例如print(a,b)),然后重新运行代码,则会得到类似以下内容的内容：

(0, 0)
(1.4901161193847656e-08, 0.0)
(0.0, 1.4901161193847656e-08)

这意味着最小化仅在这些点上评估功能.

如果现在尝试使用这3对值评估chisqfunc(),则会看到它们完全匹配,例如

print chisqfunc((0,0))==chisqfunc((1.4901161193847656e-08,0))
True

发生这种情况是因为舍入浮点算法.换句话说,当评估应力模型时,var应力比模型大太多数量级,结果被截断.

然后,可以在通过loadtxt加载数据后立即尝试执行暴力破解,以提高浮点精度,方法是写入data = data.astype(np.float128).最小化失败,结果为result.success = False,但有一条有用的消息

Desired error not necessarily achieved due to precision loss.

然后,一种可能是提供更好的初始猜测,以便在减法应力中建模-模型零件具有相同的数量级,另一种可能性是重新缩放数据,以便解决方案更接近您的初始猜测(0 ,0).

如果只是重新缩放数据,例如相对于某个应力值进行无量纲化(例如,这种材料的屈服/破裂),则效果会更好

这是拟合的示例,使用最大测量应力作为应力标度.您的代码几乎没有更改：

import numpy
import scipy.optimize as opt

filename = 'data.csv'

data = numpy.loadtxt(open(filename,"r"),delimiter=",")

stress = data[:,0]
strain = data[:,1]
err_stress = data[:,2]


smax = stress.max()
stress = stress/smax
#I am assuming the errors err_stress are in the same units of stress.
err_stress = err_stress/smax

def chisqfunc((a, b)):
    model = a + b*strain
    chisq = numpy.sum(((stress - model)/err_stress)**2)
    return chisq

x0 = numpy.array([0,0])

result =  opt.minimize(chisqfunc, x0)
print result
assert result.success==True
a,b=result.x*smax
plot(strain,stress*smax)
plot(strain,a+b*strain)

您的线性模型非常好,也就是说,在此变形范围内,您的材料具有非常线性的行为(反正是什么材料？)：

标签：minimization,scipy,statistics,chi-squared,python
来源： https://codeday.me/bug/20191122/2056442.html