本段代碼可實現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 }
>>
