Python实现曲线拟合的最小二乘法
作者:努力写代码的瀟 发布时间:2022-03-24 16:13:25
标签: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)
来源:https://blog.csdn.net/weixin_45779228/article/details/109318170
0
投稿
猜你喜欢
- 作者:Dmitry @ Usability Post 版权所有 Copyright.译者:明月星光 @ UCD翻译小组原文:ht
- [Hack] 意为”劈”、”砍”。 [Hacker] 意为”黑客”CSS Hack 是指针对不同的浏览器写不同的CSS code的过程,简单
- 最近我因为要安装 Firebug 1.4 导致我不得不安装了 Firefox 3.5 ,所以很不小心地接触到了Wordpress后台那漂亮的
- 当然可以,我们使用强大的fso对象来获取文件夹的大小请敲入如下代码即可:<%Set MyFileSize =&nb
- 在对二维数据进行 resize / mapping / 坐标转换等操作时,经常会将原本的整数坐标变换为小数坐标,对于非整数的坐标值一种直观有
- xhtml+css页面制作过程中问题的解决方案,说是解决方案应该有点过了,充其量只不过是给刚刚开始学标准页面制作的朋友们的一些小建议,如果讲
- 在PHP界谈模板引擎,必不可免的要拿Smarty开刀, 这个无比傻帽的却又带有一点点官方色彩的模板引擎, 如果没有我这样人富有正义感又富有创
- (注:在看到大家如此关注JS里头的这几个对象,我试着把原文再修改一下,力求能再详细的阐明个中意义 2007-05-21)在提到上述的概念之前
- 这两天终于忍不住的去实验了一下,为什么网页的字体有时会显示成超级无敌难看的宋体呢?其实宋体不难看,难看的只是把它放在Leopard下,没有点
- 有时会被问到“看看XXX网站如何?”之类的问题。谈到评估,通常都是指产品级的网站,如果模式很新,了解需要花一定时间。于是,很多人又问“那么你
- 见过很多网站,在设计的时候给了用户很大的自由度,我个人并不赞同这种做法。最简单的例子,圈网。我在研究圈网的时候注册完成后圈网给了我一个搜索框
- SELECT表名=case when a.colorder=1 then d.name else '' end,表说明=ca
- 是否看见大站的广告都是放在内容中间实现文字环绕的呢,一般普通小站广告只能放在内容开头或者结尾,也许大站的cms系统带这个功能吧,我们小站常用
- 富文本编辑器,Rich Text Editor, 简称 RTE, 它提供类似于 Microsoft Word 的编辑功能,容易被不会编写 H
- 最近在学习PHP,以下是看PHP100视频教程,做的学习笔记,在这里存放以便今后使用。apache--PHP--DB(mysql)一、apa
- 靓丽的网页是怎样生成的?也许您会脱口而出,当然是自己设计出来的。没错!不过这其中也有网页制作工具的一部分功劳,因为功能强大的网页制作工具可以
- 应用OpenCV和Python进行SIFT算法的实现如下图为进行测试的gakki101和gakki102,分别验证基于BFmatcher、F
- 一:操作redis1:redis拓展安装composer require predis/predis或者你也可以通过 PECL 安装&nbs
- 需要在两个文件中实现:首先,在talker.asp(在线名单)中做如下处理:<%p1=trim(application("v
- Python数据类型之间的转换函数描述int(x [,base])将x转换为一个整数long(x [,base] )将x转换为一个长整数fl