Python实现曲线拟合的最小二乘法_Python

本文实例为大家分享了Python曲线拟合的最小二乘法，供大家参考，具体内容如下

模块导入

				?

									import numpy as np

									import gaosi as gs

代码

				?

									"""

									本函数通过创建增广矩阵，并调用高斯列主元消去法模块进行求解。

									"""

									import numpy as np

									import gaosi as gs

									shape = int(input('请输入拟合函数的次数:'))

									x = np.array([0.6,1.3,1.64,1.8,2.1,2.3,2.44])

									y = np.array([7.05,12.2,14.4,15.2,17.4,19.6,20.2])

									data = []

									for i in range(shape*2+1):

									 if i != 0:

									 data.append(np.sum(x**i))

									 else:

									 data.append(len(x))

									b = []

									for i in range(shape+1):

									 if i != 0:

									 b.append(np.sum(y*x**i))

									 else:

									 b.append(np.sum(y))

									b = np.array(b).reshape(shape+1,1)

									n = np.zeros([shape+1,shape+1])

									for i in range(shape+1):

									 for j in range(shape+1):

									 n[i][j] = data[i+j]

									result = gs.Handle(n,b)

									if not result:

									 print('增广矩阵求解失败！')

									 exit()

									fun='f(x) = '

									for i in range(len(result)):

									 if type(result[i]) == type(''):

									 print('存在自由变量！')

									 fun = fun + str(result[i])

									 elif i == 0:

									 fun = fun + '{:.3f}'.format(result[i])

									 else:

									 fun = fun + '+{0:.3f}*x^{1}'.format(result[i],i)

									print('求得{0}次拟合函数为：'.format(shape))

									print(fun)

高斯模块

				?

									# 导入 numpy 模块

									import numpy as np

									# 行交换

									def swap_row(matrix, i, j):

									 m, n = matrix.shape

									 if i >= m or j >= m:

									 print('错误! : 行交换超出范围 ...')

									 else:

									 matrix[i],matrix[j] = matrix[j].copy(),matrix[i].copy()

									 return matrix

									# 变成阶梯矩阵

									def matrix_change(matrix):

									 m, n = matrix.shape

									 main_factor = []

									 main_col = main_row = 0

									 while main_row < m and main_col < n:

									 # 选择进行下一次主元查找的列

									 main_row = len(main_factor)

									 # 寻找列中非零的元素

									 not_zeros = np.where(abs(matrix[main_row:,main_col]) > 0)[0]

									 # 如果该列向下全部数据为零，则直接跳过列

									 if len(not_zeros) == 0:

									 main_col += 1

									 continue

									 else:

									 # 将主元列号保存在列表中

									 main_factor.append(main_col)

									 # 将第一个非零行交换至最前

									 if not_zeros[0] != [0]:

									 matrix = swap_row(matrix,main_row,main_row+not_zeros[0])

									 # 将该列主元下方所有元素变为零

									 if main_row < m-1:

									 for k in range(main_row+1,m):

									 a = float(matrix[k, main_col] / matrix[main_row, main_col])

									 matrix[k] = matrix[k] - matrix[main_row] * matrix[k, main_col] / matrix[main_row, main_col]

									 main_col += 1

									 return matrix,main_factor

									# 回代求解

									def back_solve(matrix, main_factor):

									 # 判断是否有解

									 if len(main_factor) == 0:

									 print('主元错误,无主元！ ...')

									 return None

									 m, n = matrix.shape

									 if main_factor[-1] == n - 1:

									 print('无解！ ...')

									 return None

									 # 把所有的主元元素上方的元素变成0

									 for i in range(len(main_factor) - 1, -1, -1):

									 factor = matrix[i, main_factor[i]]

									 matrix[i] = matrix[i] / float(factor)

									 for j in range(i):

									 times = matrix[j, main_factor[i]]

									 matrix[j] = matrix[j] - float(times) * matrix[i]

									 # 先看看结果对不对

									 return matrix

									# 结果打印

									def print_result(matrix, main_factor):

									 if matrix is None:

									 print('阶梯矩阵为空！ ...')

									 return None

									 m, n = matrix.shape

									 result = [''] * (n - 1)

									 main_factor = list(main_factor)

									 for i in range(n - 1):

									 # 如果不是主元列，则为自由变量

									 if i not in main_factor:

									 result[i] = '(free var)'

									 # 否则是主元变量，从对应的行，将主元变量表示成非主元变量的线性组合

									 else:

									 # row_of_main表示该主元所在的行

									 row_of_main = main_factor.index(i)

									 result[i] = matrix[row_of_main, -1]

									 return result

									# 得到简化的阶梯矩阵和主元列

									def Handle(matrix_a, matrix_b):

									 # 拼接成增广矩阵

									 matrix_01 = np.hstack([matrix_a, matrix_b])

									 matrix_01, main_factor = matrix_change(matrix_01)

									 matrix_01 = back_solve(matrix_01, main_factor)

									 result = print_result(matrix_01, main_factor)

									 return result

									if __name__ == '__main__':

									 a = np.array([[2, 1, 1], [3, 1, 2], [1, 2, 2]],dtype=float)

									 b = np.array([[4],[6],[5]],dtype=float)

									 a = Handle(a, b)