多元线性回归代码

多元线性回归代码

%% 代码说明
% 该程序为多元线性回归的实现，仅供学习参考，切勿抄袭
% 输入：
%       X：为p个变量的已知值，若样本数为n，则X是一个n行p列的矩阵
%       Y：为1个变量的已知值，是一个列向量
%       alpha：显著水平，不输入则默认为0.05
% 输出：
%       coef：一元线性方程的系数
%       R_square：判定系数，越接近1，拟合越好；越接近0，拟合越差
%       F：F检验值
%       F_result：检验结果，若拒绝H0，表明存在显著线性关系

%% 代码实现
function [coef, R_square, F, F_result] = MultipleLinearRegression(X, Y, alpha)
% 判断是否输入显著水平
if nargin == 2
    alpha = 0.05;
end

% 一些中间值
[n, p] = size(X); % 样本数
y_bar = sum(Y) / n; % y的平均值
X = [ones(n, 1), X];

% 求解回归系数coef
coef = inv(X' * X) * (X' * Y);

% 求解判定系数
SST = sum((y_bar * ones(n, 1) - Y) .* (y_bar * ones(n, 1) - Y));
Y_hat = X * coef;
SSR = sum((Y_hat-y_bar) .* (Y_hat-y_bar));
R_square = SSR / SST; % 判定系数

% 线性关系检验
SSE = sum((Y - Y_hat) .* (Y - Y_hat));
F = (SSR / p) / ((SSE) / (n - p - 1));
v = finv(1 - alpha, p, n - p - 1);
if F > v
    F_result = "拒绝H0";
else
    F_result = "接受H0";
end

% 残差检验
E = Y - Y_hat; % 残差
figure
scatter(1:n, E');
hold on
plot(1:n, zeros(1, n), '--')
title("残差图")

end