向中间靠拢——线性回归分析

---

联系数计加油站获取论文和代码等资源。
邮箱:MathComputerSharing@outlook.com
微信公众号:数计加油站
致谢:清风数学建模

---

何谓回归？

回归分析是一个很高大上的词。我们姑且可以这样理解它:数据是变化莫测的，一个因变量可能受多个因变量的影响，但最终这种影响能够回归于一种稳态。

举一个栗子。孩子的身高可能受父母身高的影响。假设有两个妈妈的身高相同，但一个爸爸很高，一个爸爸很矮，常常出现这样的情况:爸爸高的孩子比爸爸矮一点，爸爸矮的孩子比爸爸高一点，他们的身高都接近于一个平均身高。这就是回归现象。抽象而言，就是固定其他自变量，只改变一个自变量，因变量的值在一个平均值附近小幅度偏差波动。

回归分析，就是找到一种函数关系，来昭示回归现象，我们能够使用它来估计和预测特定自变量下因变量的平均值。

线性回归

最简单的函数关系当然是线性关系。在这种假定下，因变量可以分解为若干个自变量的线性组合。线性组合的系数称为回归系数。

$$y={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_n{x}_n+\mu$$

回归系数的含意是，固定其他自变量，该自变量改变一个单位，因变量的平均值被影响而改变的幅度。

应当注意，上式中的$\mu$即为导致因变量波动的未知因素。我们可以通过最小化$\mu$的2－范数的平方，来获得回归系数。这通常称为最小二乘法。

向中间靠拢

比较好的回归结果，应当使数据尽量接近回归方程，或者说，应当让数据向中间(回归方程)靠拢。这也是回归分析的本质。

一篇简短的数学建模论文

Matlab 代码

蒙特卡洛模拟

clc,clear

% 构造数据的函数
f=@(x1,x2) 2*x1+5*x2+0.5;
% 蒙特卡洛模拟次数
times = 10000;
n = 100;

% 记录相关系数和回归系数
r = ones(times,1);
k = ones(times,1);
for i=1:times
	% 随机数据x1
	x1 = -10+rand(n,1)*20;
	% 与x1相关的遗漏变量x2
	x2 = 0.3*x1+5*randn(size(x1,1),1);
	% 原扰动误差
	u0 = randn(n,1);
	y = f(x1,x2)+u0;
	% 遗漏变量x2进行回归
	theta = LinearRegression(x1,y);
	% 现扰动误差
	u1 = u0+5*x2;
	tmp = corrcoef(x1,u1);
	r(i) = tmp(1,2);
	k(i) = theta(2);
end

%% 绘图
scatter(r,k);
xlabel("r_{x_1,\mu}","Interpreter","tex");
ylabel("k",Interpreter="tex",Rotation=0);

自定义的线性回归函数

function [theta,values] = LinearRegression(X,y)
% 线性回归
    X = [ones(size(X,1),1),X];
    theta = (transpose(X)*X)\transpose(X)*y;
    values = X*theta;
end

线性回归简单示例

clc,clear

%% 读取数据
data = readmatrix("data\高数成绩.csv","NumHeaderLines",1);
X = data(:,1:3);
y = data(:,4);

%% 普通线性回归
[theta,values] = LinearRegression(X,y);
disp(theta);

%% 标准化线性回归
norm_X = normalize(X,1,"zscore");
norm_y = normalize(y,1,"zscore");
norm_theta = LinearRegression(norm_X,norm_y);
disp(norm_theta);

%% 比较真实值与回归值
plot(1:size(data,1),values);
hold on;
plot(1:size(data,1),y);
legend("回归值","真实值");

Python代码

初始化

# 导包
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler

plt.rcParams["font.sans-serif"] = ["SimHei"]  #设置字体
plt.rcParams["axes.unicode_minus"] = False  # 解决图像中的“-”负号的乱码问题

# 导入数据
data = pd.read_csv("../data/高数成绩.csv", encoding="gbk")
display(data)

# 数据初始化
D = data.to_numpy()
X = D[:,0:3]
y = D[:,3]
display(X)
display(y)

线性回归简单示例

# 普通线性回归
model = LinearRegression()
model.fit(X, y)
print(model.coef_)
print(model.intercept_)
print(model.score(X, y))

# 比较真实值与回归值
y_predict = model.predict(X)
plt.plot(list(range(1, 21)), y)
plt.plot(list(range(1, 21)), y_predict)
plt.xticks(list(range(1, 21)))
plt.legend(["真实值", "预测值"])

# 标准化线性回归
std = StandardScaler()
n_X = std.fit_transform(X)
n_y = std.fit_transform(y.reshape(-1, 1))
n_model = LinearRegression()
n_model.fit(n_X, n_y)
print(n_model.coef_)
print(n_model.intercept_)
print(n_model.score(n_X, n_y))

蒙特卡洛模拟探究内生性

# 构造数据的函数
f = lambda a, b: 2 * a + 5 * b + 0.5
# 蒙特卡洛模拟次数
times = 10000
n = 100

# 记录相关系数和回归系数
r = np.ones((times, 1))
k = np.ones((times, 1))
for i in range(times):
    # 随机数据x1
    x1 = -10 + np.random.rand(n, 1) * 20
    # x1相关的遗漏变量x2
    x2 = 0.3 * x1 + 5 * np.random.randn(x1.shape[0], 1)
    # 原扰动误差
    u0 = np.random.randn(n, 1)
    y = f(x1, x2) + u0
    # 遗漏变量x2进行回归
    m_model = LinearRegression()
    m_model.fit(x1, y)
    # 现扰动误差
    u1 = u0 + 5 * x2
    tmp = np.corrcoef(x1[:, 0], u1[:, 0])
    r[i] = tmp[0, 1]
    k[i] = m_model.coef_[0]

# 绘图
plt.scatter(r, k, 1)
plt.xlabel("$r_{x_1,\mu}$")
plt.ylabel("k")