本文實例為大家分享了Python曲線擬合的最小二乘法，供大家參考，具體內容如下

模塊導入

import numpy as npimport gaosi as gs

代碼

'''本函數通過創建增廣矩陣，并調用高斯列主元消去法模塊進行求解。'''import numpy as npimport gaosi as gsshape = 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 resultif __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)

以上就是本文的全部內容，希望對大家的學習有所幫助，也希望大家多多支持好吧啦網。