python基于最小二乘法计算拟合

Posted 2021-06-29 _雪辉_

导入csv数据，计算拟合曲线，拟合度

#!/usr/bin/python3
# -*- coding: utf-8 -*-

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

#计算拟合度函数
def __sst(y_no_fitting):
    """
    计算SST(total sum of squares) 总平方和
    """
    y_mean = sum(y_no_fitting) / len(y_no_fitting)
    s_list =[(y - y_mean)**2 for y in y_no_fitting]
    sst = sum(s_list)
    return sst


def __ssr(y_fitting, y_no_fitting):
    """
    计算SSR(regression sum of squares) 回归平方和
    """
    y_mean = sum(y_no_fitting) / len(y_no_fitting)
    s_list =[(y - y_mean)**2 for y in y_fitting]
    ssr = sum(s_list)
    return ssr


def __sse(y_fitting, y_no_fitting):
    """
    计算SSE(error sum of squares) 残差平方和
    """
    s_list = [(y_fitting[i] - y_no_fitting[i])**2 for i in range(len(y_fitting))]
    sse = sum(s_list)
    return sse


def goodness_of_fit(y_fitting, y_no_fitting):
    """
    计算拟合优度R^2
    """
    SSR = __ssr(y_fitting, y_no_fitting)
    SST = __sst(y_no_fitting)

    rr = SSR /SST
    return rr

plt.style.use('ggplot')
#原数据
df = pd.read_csv('./test.csv',encoding='gbk')
#转置原数据
ori = df.T
#真实时间字段
xcordori = ori.loc['x']
#真实输出电压
ycordori = ori.loc['y']

#转换真实数据为时间输出电压字典
oritimelist = (list(xcordori))
oridatalist = (list(ycordori))
oridict = {}
for i in range(len(oritimelist)):
    oridict[oritimelist[i]] = oridatalist[i]

#剔除重复数据（生成宏观函数）
dfd = df.drop_duplicates(subset='y', keep='last')
#转置宏观数据
dft = dfd.T
#宏观时间字段
xcord = dft.loc['x']
#宏观输出电压
ycord = dft.loc['y']

#真实图像
plt.scatter(xcordori,ycordori,s=6,c='blue',marker='s')
# 宏观数据
#plt.scatter(xcord,ycord,s=6,c='green',marker='s')

##最小二乘法计算拟合
## xy 的均值
(xcord*ycord).mean()
## x 的均值乘以 y 的均值
xcord.mean()* ycord.mean()
## x 的平方均值
pow(xcord,2).mean()
## x 的均值的平方
pow(xcord.mean(),2)
# m 分子是 xy 的均值减去 x 的均值乘以 y 的均值；
# m 分母是 x 平方的均值 减去 x 的均值的平方
m = ((xcord*ycord).mean() - xcord.mean()* ycord.mean())/(pow(xcord,2).mean()-pow(xcord.mean(),2))
# c 等于 y 的均值 - m 乘以 x 的均值
c = ycord.mean() - m*xcord.mean()
x=np.arange(min(xcord),max(xcord),0.025)
y=m*x+c
plt.plot(x,y)
plt.show()
print("拟合函数：y = %sx + %s" %(m,c))

##求相交点
y2 = []
for i in range(len(y)):
    y2.append(float('%.2f' % y[i]))

ycord2 = []
for j in range(len(list(ycord))):
    ycord2.append(float('%.2f' % list(ycord)[j]))
virdata = set(y2)
dirdata = set(ycord2)
result = virdata & dirdata
minresult = min(result)
maxresult = max(result)
print("最小相交点：",minresult)
print("最大相交点：",maxresult)
ymin = minresult
ymax = maxresult
xmax = (ymin-c)/m
xmin = (ymax-c)/m 

#拟合数据
res = []
for i in x:
    if i >= xmin and i <= xmax:
        i = (float('%.3f' % i))
        yres=m*i+c
        res.append(float('%.2f' % yres))

#真实数据
tot = []
for i in oritimelist:
    if i >= xmin and i <= xmax:       
        tot.append(float('%.2f' % oridict[i]))

#计算拟合
rr = goodness_of_fit(res,tot)
print("拟合度："+str(rr))