曲线拟合的最小二乘法:最小二乘法 python

import numpy as np
def matDot(A,B):
    C=[]
    row=len(A)
    col=len(B[0])
    for i in range(row):
        oneRow=col*[0]
        for j in range(col):
            s=0.0
            for k in range(len(B)):
                s+=A[i][k]*B[k][j]
            oneRow[j]=s
        C.append(oneRow)
    return C                         #定义矩阵相乘函数



transpose=lambda X:list(map(list,zip(*X)))    #定义矩阵转置函数



def augmentMatrix(A, b):
    if(len(A) != len(b)):
        raise 'The number of rows is different'
    result = []
    for i in range(len(A)):
        row = []
        for j in range(len(A[i])):
            row.append(A[i][j])
        for j in range(len(b[i])):
            row.append(b[i][j])
        result.append(row)       
    return result                 #生成增广矩阵


def GJ(a):                      #定义高斯列主元素消去法
    row=a.shape[0]    
    for j in range(0,row):
        if j<row:
            b=FindLarge(a[j:,j])
        else:
            b=0
        b1=b+j                      # 主元和对角线所在行交换
        c= np.copy(a[b1,:])
        a[b1,:]=a[j,:]
        a[j,:]=c
        for i in range(0, row):
            if i==j:
                continue
            a[i,:]=a[i,:]-a[j,:]*a[i,j]/a[j,j]
    return a
def FindLarge(a0):          #寻找主元
    b0=np.argmax(a0)
    return b0  


x=[]
y=[]
C=[]
CT=[]
A=[]
b=[]
Y=[]
a=[]
n=int(input("请输入最高项的次数:"))
xxx=(input("请输入x的值,以空格分隔:"))
yyy=(input("请输入y的值,以空格分隔:"))

x=xxx.split(" ")
m=len(x)
y=yyy.split(" ")
m1=len(x)
x=[float(x) for x in x]
y=[float(x) for x in y]

for i in range(m):
    onerow=[]
    onerow1=[0]    
    for j in range(n+1):
        onerow.append(0)
    C.append(onerow)
    Y.append(onerow1)
for i in range(m):
    Y[i][0]=y[i]
    for j in range(n+1):
        C[i][j]=x[i]**(j)    
CT=transpose(C)
A=matDot(CT,C)
b=matDot(CT,Y)
AA1=augmentMatrix(A,b)
A1 = np.mat(AA1, dtype=float)
a1=GJ(A1)
for i in range(len(A)):
    xxx=a1[i,len(A)]/a1[i,i]
    yyy=float(xxx)            #结果保留四位小数
    a.append(yyy)
    print("a%d="%(i),yyy)
print('fai(x)=',a[0],end='')
for i in range(1,n+1):
    print('+',a[i],"x**%d"%(i),end='')

输入案例
请添加图片描述

  • 0
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值