Python - 最小化卡方值

Question

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

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

谢谢，Will

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

Answer 1

问题在于，你最开始的猜测离实际答案太远了。如果你在 chisqfunc() 里面加一句打印语句，比如 print (a,b)，然后重新运行你的代码，你会看到类似下面的内容：

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

这意味着 minimize 只在这些点上计算函数。

如果你现在尝试在这三对值上计算 chisqfunc()，你会发现它们完全一致，比如：

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

这种情况发生是因为浮点数运算的四舍五入。换句话说，当计算 stress - model 时，变量 stress 的数值远远大于 model，导致结果被截断。

你可以尝试一种简单的方法，增加浮点数的精度，在用 loadtxt 加载数据后，写上 data=data.astype(np.float128)。虽然 minimize 失败了，result.success=False，但它给出了一个有用的提示：

由于精度损失，期望的误差未必能达到。

一种解决办法是提供一个更好的初始猜测，这样在计算 stress - model 时，model 的数值就能和 stress 在同一个数量级。另一种办法是对数据进行重新缩放，这样得到的解会更接近你最初的猜测 (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)

你的线性模型相当不错，也就是说，在这个变形范围内，你的材料表现得非常线性（顺便问一下，这是什么材料？）：