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 

感謝您閱讀我的問題和任何幫助將不勝感激。

乾杯，威爾

編輯：Link to data

Answer 1

的問題是，你最初的猜測是從實際的解決方案很遠：數據我目前使用的設置。如果添加打印語句中chisqfunc()像print (a,b)，並重新運行你的代碼，你會得到這樣的：

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

這意味着minimize只在這些點計算函數。

，如果你現在嘗試在這3對值來評估chisqfunc()，你會發現它們完全匹配，例如

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

這是因爲四捨五入的浮動點算術的。換句話說，在評估stress - model時，變量stress的數量級比model大得多，並且結果被截斷。

然後可以嘗試bruteforcing它，提高浮點精度，在用loadtxt加載數據之後編寫。 minimize失敗，與result.success=False，但與一個有用的消息

由於精度損失不一定達到所需的錯誤。

一種可能性是再提供一個更好的初始猜測，所以，在減法stress - model的model部分是相同的數量級，另重新調整數據，因此該解決方案將是更接近你初步猜測(0,0)。

它是MUCH更好，如果你只是重新縮放數據，使得例如無量綱相對於一定的應力值（如yelding /該材料的開裂）

這是裝配件的一個例子，將最大測量的應力用作應力等級。有從您的代碼非常少的變化：

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) 

你的線性模型相當不錯的，即你的材料具有這個範圍的變形的非常線性的行爲（什麼材料也無妨？）：