本段代码可实现OLS法的线性回归分析，并可对回归系数做出分析
1.代码
%%OLS法下的线性回归
function prodict = Linear_Regression(X,Y)
x = sym('x');
n = max(size(X));
%%定义画图窗格属性
h = figure;
set(h,'color','w');
%%回归相关值
XX_s_m = (X-Expection(X,1))*(X-Expection(X,1))';
XY_s_m = (X-Expection(X,1))*(Y-Expection(Y,1))';
YY_s_m = (Y-Expection(Y,1))*(Y-Expection(Y,1))';
%%回归系数
beta = XY_s_m/XX_s_m;
alpha = Expection(Y,1)-beta*Expection(X,1);
%%画出实际值与回归值的差
yyaxis left
for k = 1:n
delta_Y(k) = Y(k) - (alpha+beta*X(k));
end
width = (max(X)-min(X))/(2*n);
b = bar(X,delta_Y,width,'Facecolor',[0.9 0.9 0.9]);
set(b,'Edgecolor','none');
%标出差值
disp('是否标出差值？');
str = input('yes or no?','s');
Y_dis = (max(Y) - min(Y))/40;
if strcmp(str,'yes') == 1
if n<=10
for k = 1:length(delta_Y)
if delta_Y(k)>=0
text(X(k),delta_Y(k)+Y_dis,strcat(num2str(delta_Y(k),'%.3f')),...
'VerticalAlignment','middle','HorizontalAlignment','center');
else
text(X(k),delta_Y(k)-Y_dis,strcat(num2str(delta_Y(k),'%.3f')),...
'VerticalAlignment','middle','HorizontalAlignment','center');
end
end
end
end
xlabel('X');ylabel('差值');
hold on
%%画出数据点
yyaxis right
plot(X,Y,'g--*')
if n<=20
grid on;
end
hold on
ylabel('Y');
%%相关系数
R2 = (XY_s_m)^2/(XX_s_m*YY_s_m);
%%随机项方差无偏估计
sigma_2 = (Y*Y'-n*Expection(Y,1)^2-beta*(X*Y'-n*Expection(X,1)*Expection(Y,1)))/(n-2);
%%beta和alpha的方差
Var_beta = sigma_2/(X*X'-n*Expection(X,1)^2);
Var_alpha = sigma_2*(X*X')/(n*(X*X'-n*Expection(X,1)^2));
%%beta和alpha的协方差
cov_alpha_beta = -Expection(X,1)*sigma_2/(X*X'-n*Expection(X,1)^2);
%%画出回归直线
x_v = min(X)-Expection(X,1)/sqrt(n):(max(X)-min(X))/100:max(X)+Expection(X,1)/sqrt(n);
y_v = beta*x_v+alpha;
y_ave = ones(1,max(size(x_v)))*Expection(Y,1);
plot(x_v,y_v,'r',x_v,y_ave,'k','LineWidth',0.2,'LineStyle','-')
%%标出数据点
disp('是否标出数据点？');
str = input('yes or no?','s');
if strcmp(str,'yes') == 1
if n <=10
for k = 1:max(size(X))
text(X(k),Y(k),['(',num2str(X(k)),',',num2str(Y(k)),')'],...
'color','b','VerticalAlignment','middle','HorizontalAlignment','center');hold on;
end
end
end
disp('是否输出回归方程的置信区间？');
str = input('yes or no?','s');
if strcmp(str,'yes') == 1
str = input('单侧置信系数 或者 双侧置信系数？','s');
if strcmp(str(1),'单') == 1
ALPHA = input('输入单侧置信系数：');
a = 1-ALPHA*2;
elseif strcmp(str(1),'双') == 1
ALPHA = input('输入双侧置信系数：');
a = 1-ALPHA;
end
disp('t(n-2)为：');
t_value = nctinv(a,n-2,0)
addvalue = t_value*sqrt(sigma_2)*sqrt(1+1/n+(x-Expection(X,1))^2/(X*X'-n*Expection(X,1)^2));
disp('置信区间上限的方程为：');
vpa(alpha+beta*x+addvalue,4)
disp('置信区间下限的方程为：');
vpa(alpha+beta*x-addvalue,4)
y_ub = subs(alpha+beta*x+addvalue,x_v);
y_lb = subs(alpha+beta*x-addvalue,x_v);
plot(x_v,y_ub,'--',x_v,y_lb,'-.','color',[0.2 0.5 0.8]);
title('OLS法对数据的线性回归');
legend('差值','Y:实际值','y:拟合直线','平均值','置信上限','置信下限');
hold off
yyaxis left
text(min(X)-Expection(X,1)/sqrt(n),max(delta_Y)-Expection(delta_Y,1)/sqrt(n),['置信程度',num2str(100*a),'%'],'color',[0.3 0.3 0.3]);
hold off
else
legend('差值','Y:实际值','y:拟合直线','平均值');
title('OLS法对数据的线性回归');
hold off
end
disp('是否输出回归系数的置信区间？');
str = input('yes or no?','s');
if strcmp(str,'yes') == 1
str = input('单侧置信系数 或者 双侧置信系数？','s');
if strcmp(str(1),'单') == 1
ALPHA = input('输入单侧置信系数：');
a = 1-ALPHA*2;
elseif strcmp(str(1),'双') == 1
ALPHA = input('输入双侧置信系数：');
a = 1-ALPHA;
end
disp('t(n-2)为：');
t_value = nctinv(a,n-2,0)
disp('回归系数beta的置信区间为：');
[beta-t_value*sqrt(Var_beta),beta+t_value*sqrt(Var_beta)]
disp('回归系数alpha的置信区间为：');
[Expection(Y,1)-(beta+t_value*sqrt(Var_beta))*Expection(X,1),...
Expection(Y,1)-(beta-t_value*sqrt(Var_beta))*Expection(X,1)]
end
%%预测值输出
disp('是否需要求预测值？');
str = input('yes or no?','s');
if strcmp(str,'yes') == 1
X0 = input('输入预测向量：');
str = input('单侧置信系数 或者 双侧置信系数？','s');
if strcmp(str(1),'单') == 1
ALPHA = input('输入单侧置信系数：');
a = 1-ALPHA*2;
elseif strcmp(str(1),'双') == 1
ALPHA = input('输入双侧置信系数：');
a = 1-ALPHA;
end
disp('t(n-2)为：');
t_value = nctinv(a,n-2,0)
addvalue0 = t_value*sqrt(sigma_2).*sqrt(1+1/n+(X0-Expection(X,1)).^2/(X*X'-n*Expection(X,1)^2));
Y0 = alpha+beta*X0;
prodict{1,1} = {'预测点x'};
prodict{2,1} = {'预测值y'};
prodict{3,1} = {'置信区间'};
for k = 2:max(size(X0))+1
interval(k-1,:) = [Y0(k-1)-addvalue0(k-1),Y0(k-1)+addvalue0(k-1)];
prodict{1,k} = X0(k-1);
prodict{2,k} = Y0(k-1);
prodict{3,k} = interval(k-1,:);
Y0_lb(k-1) = Y0(k-1)-addvalue0(k-1);
Y0_ub(k-1) = Y0(k-1)+addvalue0(k-1);
end
h2 = figure(2);
set(h2,'color','w');
plot(X0,Y0,'r-*',X0,Y0_ub,'g-o',X0,Y0_lb,'b-o');hold on;
errorbar(X0,Y0,addvalue0,'color',[0.8 0.5 0.2]);
xlabel('X');ylabel('y');
grid on;
title('点预测');
legend('预测值','置信上限','置信下限','置信限');
Y0_dis = (max(Y0)-min(Y0))/40;
Y0_lb_dis = (max(Y0_lb)-min(Y0_lb))/40;
Y0_ub_dis = (max(Y0_ub)-min(Y0_ub))/40;
disp('是否显示数据？');
str = input('yes or no?','s');
if strcmp(str ,'yes') == 1
if max(size(X0)) <= 3
for k = 1:max(size(X0))
text(X0(k),Y0(k)+Y0_dis,['(',num2str(X0(k)),',',num2str(Y0(k)),')'],...
'color',[0.2 0.2 0.2]);hold on;
text(X0(k),Y0_lb(k)+Y0_lb_dis,['(',num2str(X0(k)),',',num2str(Y0_lb(k)),')'],...
'color',[0.2 0.2 0.2],'VerticalAlignment','middle','HorizontalAlignment','center');hold on;
text(X0(k),Y0_ub(k)+Y0_ub_dis,['(',num2str(X0(k)),',',num2str(Y0_ub(k)),')'],...
'color',[0.2 0.2 0.2],'VerticalAlignment','middle','HorizontalAlignment','center');hold on;
end
end
end
text(Expection(X0,1),(max(Y0)+max(Y0_ub))/2,['置信程度',num2str(100*a),'%'],'color',[0.3 0.3 0.3]);
else
prodict = [];
end
disp('回归方程为：');
vpa(alpha+beta*x,4)
disp('相关系数平方为：');
vpa(R2,6)
disp('随机项方差的无偏估计为：');
vpa(sigma_2,6)
disp('回归系数beta和alpha的方差分别为：');
vpa(Var_beta,6)
vpa(Var_alpha,6)
disp('回归系数beta和alpha的协方差为：');
vpa(cov_alpha_beta,6)
%%SST,SSE和SSR
disp('SST：');
SST = YY_s_m
Y_OLS = reshape(alpha+beta*X',1,n);
disp('SSR：');
SSR = (Y_OLS-Expection(Y,1))*(Y_OLS-Expection(Y,1))'
disp('SSE：');
SSE = (Y-Y_OLS)*(Y-Y_OLS)'
function En = Expection(M,n) %%n阶原点矩%%
En = sum(M.^n)/size(M,2);
end
end
2.下面以一个例子来说明程序运行结果
clear all clc X = [-2.269 -1.248 -0.269 1.248 2.006 4.598 5.264 7.002]; Y = [4582 4027 4958 5078 6072 4156 5948 5027]; S = Linear_Regression(X,Y)
是否标出差值？
yes or no?yes
是否标出数据点？
yes or no?no
是否输出回归方程的置信区间？
yes or no?yes
单侧置信系数 或者 双侧置信系数？双
输入双侧置信系数：0.05
t(n-2)为：
t_value =
1.9432
置信区间上限的方程为：
ans =
78.43*x + 1466.0*(0.013*(x - 2.041)^2 + 1.125)^(1/2) + 4821.0
置信区间下限的方程为：
ans =
78.43*x - 1466.0*(0.013*(x - 2.041)^2 + 1.125)^(1/2) + 4821.0
是否输出回归系数的置信区间？
yes or no?yes
单侧置信系数 或者 双侧置信系数？双
输入双侧置信系数：0.05
t(n-2)为：
t_value =
1.9432
回归系数beta的置信区间为：
ans =
-88.7367 245.5885
回归系数alpha的置信区间为：
ans =
1.0e+03 *
4.4796 5.1622
是否需要求预测值？
yes or no?yes
输入预测向量：[-2.156 -1.306 0.248 1.296]
单侧置信系数 或者 双侧置信系数？双
输入双侧置信系数：0.05
t(n-2)为：
t_value =
1.9432
是否显示数据？
yes or no?no
回归方程为：
ans =
78.43*x + 4821.0
相关系数平方为：
ans =
0.121668
随机项方差的无偏估计为：
ans =
569067.0
回归系数beta和alpha的方差分别为：
ans =
7400.35
ans =
101976.0
回归系数beta和alpha的协方差为：
ans =
-15107.8
SST：
SST =
3887366
SSR：
SSR =
4.7297e+05
SSE：
SSE =
3.4144e+06
OLS =
10×5 cell 数组
列 1 至 3
{'拟合斜率beta' } {'拟合截距alpha' } {'相关系数R^2' }
{[ 78.4259]} {[ 4.8209e+03]} {[ 0.1217]}
{'EX^2' } {'EY^2' } {'EXY' }
{[ 13.7799]} {[ 2.5296e+07]} {[ 1.0923e+04]}
{'min(X)' } {'min(Y)' } {'cov(X,Y)' }
{[ -2.2690]} {[ 4027]} {[ 753.8428]}
{'样本标准差SD(X)'} {'样本标准差SD(Y)'} {'总体标准差sigma(X)'}
{[ 3.3144]} {[ 745.2100]} {[ 9.6122]}
{'变异系数CV(Y)' } {'X极差' } {'Y极差' }
{[ 0.1496]} {[ 9.2710]} {[ 2045]}
列 4 至 5
{'X平均值' } {'Y平均值' }
{[ 2.0415]} {[ 4981]}
{'max(X)' } {'max(Y)' }
{[ 7.0020]} {[ 6072]}
{'X绝对值几何均值' } {'Y绝对值几何均值'}
{[ 2.0591]} {[ 4.9328e+03]}
{'总体标准差sigma(Y)'} {'变异系数CV(X)' }
{[ 4.8592e+05]} {[ 1.6235]}
{'n*E(Y-Y_ba)' } {'SSR' }
{[ -3.1238e+06]} {[ 1.2431e+12]}
prodict =
3×5 cell 数组
列 1 至 4
{1×1 cell} {[ -2.1560]} {[ -1.3060]} {[ 0.2480]}
{1×1 cell} {[4.6518e+03]} {[4.7185e+03]} {[4.8403e+03]}
{1×1 cell} {1×2 double } {1×2 double } {1×2 double }
列 5
{[ 1.2960]}
{[4.9225e+03]}
{1×2 double }
>>
posted on
浙公网安备 33010602011771号