import numpy as np
from sympy import symbols, diff, solve
from scipy import stats
# 这里是在计算损失函数
def function_2(a, b, x): # 定义一个函数,拿来储存a+bx,这里a和b是参数要最后求取的内容
return a+b*x
a, b, x, y = symbols('a b x y') # 利用symbols函数来表示变量a,b,x
y1 = np.array([0.5153, 0.5638, 0.6697, 0.7716, 0.8428, 0.8980, 0.9338, 0.9979, 1.0738, 1.2110, 1.3146, 1.6068, 1.9252,
2.4108, 3.1976, 3.7771, 4.6893, 5.6561, 6.0460, 6.9405, 8.0277, 8.8620, 9.5578]) # 因变量数据
x1 = np.array([1.8668, 2.1781, 2.6923, 3.5334, 4.8198, 6.0794, 7.1177, 7.8973, 8.4402, 8.9677, 9.9215, 10.9655, 12.0333,
13.5823, 15.9878, 18.4937, 21.6314, 26.5810, 31.4045, 34.0903, 40.1513, 47.3104, 51.8942]) #自变量数据
e1 = function_2(a, b, x1) # 调用上文函数function_2用于储存a+b*x
average_y = sum(y1)/len(y1)
average_X = sum(x1)/len(x1)
loss1 = (y1 - e1)**2 # 计算损失函数
for term in loss1:
print(term)
loss_fun = loss1.sum() # 将前文的每一项的损失函数全部求和加起来
# 这里是在求取偏导数并求得关于a和b的解
y11 = diff(loss_fun, a) # 损失函数求导得到关于a的方程
print(y11)
y22 = diff(loss_fun, b) # 损失函数求导得到关于b的方程
print(y22)
solutions = solve([y11, y22], (a, b)) # 调用solve函数,求解得到关于a和b的值
print("a =", solutions[a])
print("b =", solutions[b])
print('线性回归方程为')
print(f"y = {solutions[b]}*x+{solutions[a]}") # 使用了Python中的f-string(格式化字符串)来输出线性回归方程。在这个字符串中,我们使用花括号 {} 来表示占位符,其中的表达式将会被替换为实际的值。
# 这里进行显著性检验
head_y = solutions[b]*x1+solutions[a] # 将x1带入回归方程,计算预测值
print(head_y)
SSR = sum((head_y-average_y)**2) # 表示回归平方和,用预测值减去平均值
print('计算SSR')
print(SSR)
SSE = sum((y1-head_y)**2) # 表示残差平方和,用真实值减去预测值
print('计算SSE')
print(SSE)
SST = SSE+SSR
R = SSR/SST # 回归方程拟合度
print('计算回归方程的拟合度')
print(R)
if solutions[b]>0:
if R > 0.8:
print("正相关,拟合程度高")
elif 0.8 > R > 0.5:
print('正相关,拟合程度一般')
else:
print('正相关,拟合程度较低')
elif solutions[b] < 0:
if R > 0.8:
print("负相关,拟合程度高")
elif 0.8 > R > 0.5:
print('负相关,拟合程度一般')
else:
print('负相关,拟合程度较低')
df = len(x1)-2 # 计算自由度
t1 = sum((x1-average_X)**2)
print('计算t1')
print(t1)
Lxx = pow(t1, 0.5)
print('计算Lxx')
print(Lxx)
t_b = Lxx*solutions[b] / pow(SSE / df,0.5)
print('计算回归系数b的t统计量')
print(t_b)
alpha = 0.05
t_value = stats.t.isf(alpha/2, df) # 进行t检验,寻找对应的t值
print('计算t检验上alpha/2分位点所在的值为')
print(t_value)
if abs(t_b) > t_value:
print('拒绝原假设,即样本回归系数显著,表明存在线性关系')
else:
print('接受原假设,即样本回归系数不显著,不存在线性关系')