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)