本文实例为大家分享了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)
以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持。
广告合作:本站广告合作请联系QQ:858582 申请时备注:广告合作(否则不回)
免责声明:本站资源来自互联网收集,仅供用于学习和交流,请遵循相关法律法规,本站一切资源不代表本站立场,如有侵权、后门、不妥请联系本站删除!
免责声明:本站资源来自互联网收集,仅供用于学习和交流,请遵循相关法律法规,本站一切资源不代表本站立场,如有侵权、后门、不妥请联系本站删除!
暂无评论...
稳了!魔兽国服回归的3条重磅消息!官宣时间再确认!
昨天有一位朋友在大神群里分享,自己亚服账号被封号之后居然弹出了国服的封号信息对话框。
这里面让他访问的是一个国服的战网网址,com.cn和后面的zh都非常明白地表明这就是国服战网。
而他在复制这个网址并且进行登录之后,确实是网易的网址,也就是我们熟悉的停服之后国服发布的暴雪游戏产品运营到期开放退款的说明。这是一件比较奇怪的事情,因为以前都没有出现这样的情况,现在突然提示跳转到国服战网的网址,是不是说明了简体中文客户端已经开始进行更新了呢?
更新日志
2024年12月24日
2024年12月24日
- 小骆驼-《草原狼2(蓝光CD)》[原抓WAV+CUE]
- 群星《欢迎来到我身边 电影原声专辑》[320K/MP3][105.02MB]
- 群星《欢迎来到我身边 电影原声专辑》[FLAC/分轨][480.9MB]
- 雷婷《梦里蓝天HQⅡ》 2023头版限量编号低速原抓[WAV+CUE][463M]
- 群星《2024好听新歌42》AI调整音效【WAV分轨】
- 王思雨-《思念陪着鸿雁飞》WAV
- 王思雨《喜马拉雅HQ》头版限量编号[WAV+CUE]
- 李健《无时无刻》[WAV+CUE][590M]
- 陈奕迅《酝酿》[WAV分轨][502M]
- 卓依婷《化蝶》2CD[WAV+CUE][1.1G]
- 群星《吉他王(黑胶CD)》[WAV+CUE]
- 齐秦《穿乐(穿越)》[WAV+CUE]
- 发烧珍品《数位CD音响测试-动向效果(九)》【WAV+CUE】
- 邝美云《邝美云精装歌集》[DSF][1.6G]
- 吕方《爱一回伤一回》[WAV+CUE][454M]