python 批量做线性拟合

1. 知识点

1.1 最小二乘法：scipy.optimize.leastsq

定义好误差函数 residual, 初始点 [1.0, 0.0]，然后给参数 np.array(y), np.array(x)，调用 leastsq：

                plsq = leastsq(residual, [1.0 , 0.0], args=(np.array(y),np.array(x)) )

返还值中，plsq[0]中包括拟合结果，即对应初值[1.0,0.0] 的最佳[k,b]值。

1.2 文件写入

    fout = open(fileresult, 'w')

上面这行以写入模式打开 fileresult，返还文件输出流对象fout

                print(plsq[0][0], "    ", plsq[0][1], file=fout)

上面这行中 file=fout 表示输出至 fout，而不是默认的sys.stdout。

1.3 文件读取

with open( filedata ) as file:# read data from txt file

上面这行表示以默认的读取模式，打开 filedata，并返回文件输入流对象 file。

            line = file.readline(); words = line.split();

上面这行表示从文件输入流对象 file 读入一行，然后将这一行内容打碎为（空格分开的）字符串列表。

2. 代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

def func(x, p): # 函数形式
    k,b = p
    return k*x + b

def residual(p,y,x): # 误差函数
    return y - func(x,p)

# This does not allow blank lines in the data file, unless in the end
# This assumes in each line, data is arranged as X1 Y1 X2 Y2 ..., separated by spaces
def linearfit( filedata, fileresult ):
    results = []
    #help(fout.write); exit(1)
    count = 0
    fout = open(fileresult, 'w')
    with open( filedata ) as file:# read data from txt file
        words = [111]
        while( len(words)>0 ):
            line = file.readline(); words = line.split(); count = count+1
            if( len(words)>0 ):
                x = []; y = []
                for i in range( (int)(len(words)/2) ):
                    x.append( float(words[2*i]) ); y.append( float(words[2*i+1]) )
                x = np.array(x); y = np.array(y)
                plsq = leastsq(residual, [1.0 , 0.0], args=(np.array(y),np.array(x)) )
                print(plsq[0][0], "    ", plsq[0][1], file=fout)
            if count == 200:
                plt.plot(x, y, 'o', x, func(x, plsq[0]), '--');
                plt.title(" k = %f, b = %f" % (plsq[0][0], plsq[0][1]))
                plt.show()
    fout.close()



linearfit( "CV-3.txt", "result-CV-3.txt")