【Python】一元线性回归的分析

在科学实验中，人们研究各物理量之间的关系，通常是通过实验测得一批实验数据，然后求得自变量 x 与因变量y之间的近似函数关系 y=f(x) , 这就是按给定的 n 个点的数据拟合的问题。

线性回归（最小二乘方法）

1、今有实验数据：(共有N= 6组数据 )
| 序号 | 0 |1 | 2 |3 |4 |5|
|Xi |1 |2 |3 4 5 6
Yi 2.5 3.51 4.45 5.52 6.47 7.51

序号	0	1	2	3	4	5
Xi	1	2	3	4	5	6
Yi	2.5	3.51	4.45	5.52	6.47	7.51

假设 y = ax + b, 求 a, b=?

统计学上，为了评估实验数据x与y的线性相关性，定义相关系数：

当 r=0时，x与y完全线性无关，当 r=1时，x与y完全线性相关。

Python实现

import math
#线性拟合
#入口:
#    x = 实验数据的x向量
#    y = 实验数据的y向量
#返回:
#   a,b,r        a=斜率 b=截距 r=相关系数
#模型 y=ax + b
def linefit( x,y):N = len(x)sx,sy,sxx,syy,sxy=0,0,0,0,0for i in range(0,N):sx  += x[i]sy  += y[i]sxx += x[i]*x[i]syy += y[i]*y[i]sxy += x[i]*y[i]a = (sy*sx/N -sxy)/( sx*sx/N -sxx)b = (sy - a*sx)/Nr = abs(sy*sx/N-sxy)/math.sqrt( (sxx-sx*sx/N)*(syy-sy*sy/N))return a,b,r
X=[ 1    ,2    ,3  ,4  ,5  ,6]
Y=[ 2.5    ,3.51   ,4.45   ,5.52   ,6.47   ,7.51]
a,b,r=linefit(X,Y)
print("实验数据:")
print("X=",X)
print("Y=",Y)
print("拟合结果: y= %10.5f x + %10.5f , r=%10.5f" % (a,b,r) )

2、实验测得某段铜导线在不同温度下的电阻数据如下：
温度 T(℃) = 19.1 ,25.0 ,30.1 ,36.0,40.0 ,45.1 ,50.0
电阻 R(Ω) = 76.30,77.80 ,79.25 ,80.80 ,82.35 ,83.90 ,85.10
试按 R = a * T + b 拟合出该关系中的a,b。

#温度 T(℃) = 19.1 ,25.0 ,30.1 ,36.0,40.0 ,45.1 ,50.0
#电阻 R(Ω) = 76.30,77.80 ,79.25 ,80.80 ,82.35 ,83.90 ,85.10
import math
#线性拟合
#入口:
#    x = 实验数据的x向量
#    y = 实验数据的y向量
#返回:
#   a,b,r
#   a=斜率 b=截距 r=相关系数
#模型 y=ax + b
def linefit( x,y):N = len(x)sx,sy,sxx,syy,sxy=0,0,0,0,0for i in range(0,N):sx  += x[i]sy  += y[i]sxx += x[i]*x[i]syy += y[i]*y[i]sxy += x[i]*y[i]a = (sy*sx/N -sxy)/( sx*sx/N -sxx)b = (sy - a*sx)/Nr = abs(sy*sx/N-sxy)/math.sqrt( (sxx-sx*sx/N)*(syy-sy*sy/N))return a,b,r
T=[19.1 ,25.0 ,30.1 ,36.0,40.0 ,45.1 ,50.0 ]
R=[76.30,77.80 ,79.25 ,80.80 ,82.35 ,83.90 ,85.10]
a,b,r=linefit(T,R)
print("实验数据:")
print("T=",T)
print("R=",R)
print("拟合结果: R= %10.5f T = %10.5f , r=%10.5f" % (a,b,r) )

曲线拟合

（以前的数据）

1).可以转化为线性回归的曲线拟合方法
二十世纪六十年代，世界人口情况如下：
年份： 1960,1961,1962,1963,1964,1965,1966,1967,1968
人口(亿): 29.72,30.61,31.51,32.13,32.34,32.85,33.56,34.20,34.83
马尔萨斯人口模型:

（上式中，S代表人口, t代表年份。）
试根据上面的数据，求出 α、β的值，并推测2020年的世界人口的数量。
两边同时取自然对数，变换为线性模型:

令 a=β ， b=lnα ， y=lnS , 则转化为 y= a * t + b的线性形式。

#二十世纪六十年代，世界人口情况如下：
#年份：    1960,1961,1962,1963,1964,1965,1966,1967,1968
#人口(亿):   29.72,30.61,31.51,32.13,32.34,32.85,33.56,34.20,34.83
#马尔萨斯人口模型: S = αexp(βt)import math
#线性拟合
#入口:
#    x = 实验数据的x向量
#    y = 实验数据的y向量
#返回:
#   a,b,r
#   a=斜率 b=截距 r=相关系数
#模型 y=ax + b
def linefit( x,y):N = len(x)sx,sy,sxx,syy,sxy=0,0,0,0,0for i in range(0,N):sx  += x[i]sy  += y[i]sxx += x[i]*x[i]syy += y[i]*y[i]sxy += x[i]*y[i]a = (sy*sx/N -sxy)/( sx*sx/N -sxx)b = (sy - a*sx)/Nr = abs(sy*sx/N-sxy)/math.sqrt( (sxx-sx*sx/N)*(syy-sy*sy/N))return a,b,r
T = [1960,1961,1962,1963,1964,1965,1966,1967,1968]
S = [29.72,30.61,31.51,32.13,32.34,32.85,33.56,34.20,34.83]
Y = [ math.log(S[i]) for i in range(0,len(S))]
a,b,r=linefit(T,Y)    #注：此处的 a 是模型中的β, b是 模型中的lnα
print("实验数据:")
print("T=",T)
print("S=",S)
print("Y=",Y)
print("拟合结果: lnS= %10.5f T + %10.5f , r=%10.5f" % (a,b,r) )
y2020 = a*2020+b
s2020 = math.exp(y2020)
print("推测2020年世界人口= %10.2f 亿" % s2020 )

本题 n=9 , n-2=7, 查相关系数检验表，rr=0.66 , r = 0.9912>0.660,即

显著相关，说明马氏模型

显著成立。
说明：马氏公式仅适用于较短的时间段（如1 00年内），如用此式计算3000年的人口，将是7509206502.52亿，显然这是不可能的。如果使用2961~2968的世界人口，得出的公式，推测3000年的人口，应该比较准确。