最小二乘法

最新推荐文章于 2024-10-06 12:05:39 发布

San Chang

最新推荐文章于 2024-10-06 12:05:39 发布

阅读量64

点赞数

文章标签：线性代数矩阵

本文链接：https://blog.csdn.net/qq_26657001/article/details/87989113

版权

import numpy as np
from scipy.optimize import leastsq
#from scipy.integrate import quad,dblquad,nquad
import math
import xlrd
import matplotlib as mat
mat.use("TkAgg")
import matplotlib.pyplot as plt


#t=[]
t=np.linspace(0,4,50)
x=[]
y=[]

data = xlrd.open_workbook(r'D:/2.xlsx')
table = data.sheet_by_name(u'Sheet1')
num=37#需要变化
for i in range(num):
    #t.append(table.cell_value(i,0))
    x.append(table.cell_value(i+1,9))
    y.append(table.cell_value(i+1,11))
#y=y/np.max(y)
t=np.asarray(t)
x=np.asarray(x)
y=np.asarray(y)
print(t,x,y,len(t))
'''    
rows = table.row_values(0)
cols = table.col_values(0)
cell_A1 =
cell_C4 = table.cell_value(3,2)



G_t=un_G*(math.e**(-un_beta*t))#关于t的一个序列
juanji=np.convolve(G_t,x)#一个关于t的序列长度为t的二倍-1
sigmod_e=2*un_mill*((1+x)**(un_alpha-1)-(1+x)**((-un_alpha/2)-1))/un_alpha#一个序列
fun_sum=sigmod_e+juanji#sigmod的一个序列
'''
def func(p,t,x):
    un_alpha, un_mill, un_G, un_beta =p #[0, 0, 0, 0]
    #print(2*un_mill*((1+x)**(un_alpha-1)-(1+x)**((-un_alpha/2)-1))/un_alpha + np.convolve(un_G*(math.e**(-un_beta*t)),x)[0:num])
    return 2*un_mill*((1+x)**(un_alpha-1)-(1+x)**((-un_alpha/2)-1))/un_alpha + np.convolve(un_G*(math.e**(-un_beta*t)),x)[0:num]

def error(p,t,x,y,s):
    print (s)
    return func(p,t,x)-y

p0=[12.32,0.01,0.5,50]

s="Test the number of iteration" #试验最小二乘法函数leastsq得调用几次error函数才能找到使得均方误差之和最小的a~c
Para=leastsq(error,p0,args=(t,x,y,s)) #把error函数中除了p以外的参数打包到args中
un_alpha, un_mill, un_G, un_beta=Para[0]
print("un_alpha:", un_alpha ,"un_mill:",un_mill, "un_G:",un_G, "un_beta:",un_beta)
plt.figure()
plt.plot(x,y,"red")
plt.plot(x,func(Para[0],t,x),"green")
plt.show()