4

我一直在尝试通过最小化卡方来将线性模型拟合到一组应力/应变数据。不幸的是，使用下面的代码不能正确地最小化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

感谢您阅读我的问题，任何帮助将不胜感激。

干杯，威尔

编辑：我目前使用的数据集：链接到数据

4

1 回答 1

5

问题是您最初的猜测与实际解决方案相去甚远。如果你在chisqfunc()like中添加一个 print 语句print (a,b)，然后重新运行你的代码，你会得到类似的东西：

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

这意味着minimize仅在这些点评估函数。

如果您现在尝试评估chisqfunc()这 3 对值，您会发现它们完全匹配，例如

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

发生这种情况是因为舍入浮点算术。换句话说，在求值时stress - model，varstress比 大太多数量级model，结果被截断。

然后可以尝试对其进行暴力破解，提高浮点精度，并data=data.astype(np.float128)在加载数据后立即写入loadtxt. minimize失败，带有result.success=False，但带有有用的消息

由于精度损失，不一定能达到所需的误差。

然后一种可能性是提供更好的初始猜测，以便在减法中stress - model该model部分具有相同的数量级，另一种可能性是重新调整数据，以便解决方案更接近您的初始猜测(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)

您的线性模型非常好，即您的材料在此变形范围内具有非常线性的行为（无论如何它是什么材料？）： 在此处输入图像描述

于 2014-03-05T16:17:33.227 回答