一元二次曲线拟合Excel用

import numpy
import xlrd
import matplotlib.pyplot as plt

from scipy.optimize import leastsq

 

X = numpy.array(range(0, 19))#2001,2019 
# X = numpy.array(range(0, 30)) #x里面保存是0━29之间的所有整数

 

# 读取excel数据
def input_data(wbname):
    wb = xlrd.open_workbook(wbname)
    sh = wb.sheet_by_index(0)
                
    table = [] # 存放每年被引频次
    sumy = [] # 总被引频次
    sum5 = [] # 发表最初五年被引频次
    #sum6_9 = [] # 随后四年???
    for row in range(1, sh.nrows):#逐行读取
        Y = sh.row_values(row)[32: 32 + 19]
        # Y = sh.row_values(row)[21: 21 + 30] #y里面保存的是这30年每年文献的被引次数。
        table.append(Y)
        sumy.append(sh.row_values(row)[19]) # 合计引用次数
        sum5.append(sh.row_values(row)[32: 32 + 5]) # 发表最初五年被引频次
        #sum6_9.append(sh.row_values(row)[32 + 5: 32 + 5 + 4]) # 随后四年
    return table, sumy, sum5
    #, sum6_9

 

# 最小二乘法误差
def RSquare(X, Y, f):
    yba = 0.0
    for y in Y:
        yba += y
    yba /= len(Y)

    SStot = 0.0
    for y in Y:
        SStot += (y - yba) ** 2

    SSres = 0.0
    for x, y in zip(X, Y):
        fx = f(x)
        SSres += (y - fx) ** 2

    return 1 - SSres / SStot

 

if __name__ == '__main__':
    # table, sumlist = input_data("~/Desktop/毕业论文/引文/2001.xlsx")
    table, sumlist, sum5list = input_data("2001.xlsx")

    def residuals(p, x, y): # 定义计算误差函数
        polyfun = numpy.poly1d(p) # numpy.poly1d构造一个多项式
        return y - polyfun(x)
    
    def test_func(x, p1):
        a, b, c = p1
        return a*x**2+b*x+c
    
    #labels = []
    
    i = 2
    #j = 0
    for Y,sumy,sum5 in zip(table, sumlist,sum5list):
        #r = scipy.optimize.leastsq(residuals, [0, 0, 0], args=(X, numpy.array(Y))) # 最小二乘法
        r = leastsq(residuals, [0, 0, 0], args=(X, numpy.array(Y)))# 计算误差函数,拟合参数初始值,需要拟合的实验数据
        
        if r[1] not in [1, 2, 3, 4]:
            print("Error!")
            print(r)
            continue
        # print(r)
        p = r[0] # leastsq输出的第一个参数:拟合解组成的数组
        rsquare = RSquare(X, Y, numpy.poly1d(p))

        if rsquare >= 0.75 and sumy >= 20 and np.mean(sum5) <= 2: 
        # r方大于0.75,总被引频次大于20,发表最初五年年均被引频次不超过两次,随后四年至少为20次
            a, b, c = p

            t = -b / (2 * a)# 对称轴
            if 2.5 <= t <= 5.5: # 对称轴限定2.5 5.5
                if a > 0:# 开口向上
                    print("{}:   {}   {}".format(i, rsquare, t)) # 第几行,r方,对称轴
                    print("{}:   {}   ".format(i, p)) # 第几行,拟合的参数
                    #print(numpy.poly1d(p)) # numpy.poly1d(p)是一个函数
                                      
                    colorlist=plt.cm.cool(np.linspace(0,1,7))
                    markers=["D","x","h",".","^",">","v"]
                    #labels.append(r'$y =%ix**2 + %ix + %i$'% (p[0],p[1],p[2]))
                    
                    for n in range(0,4):# 几张图
                        plt.plot(X, numpy.array(Y),color=colorlist[2],marker=markers[0],linestyle="-") # 被引频次图,蓝色是真实的,橙色的是拟合出来的
                        plt.plot(X, test_func(X, p),color=colorlist[2],marker=markers[2],linestyle="--") # 一元二次函数
                        #plt.legend(labels)
                    plt.show()
                    #j+=1 # 改变形状颜色
                    
        i += 1