MATLAB
線性拟合小結 —— REGRESS多元線性回歸(用最小二乘估計法)
http://wenku.baidu.com/view/0a0ea0de941ea76e59fa0418.html?re=view
在回歸分析中，如果有兩個或兩個以上的自變量，就稱為多元回歸。如果隻有一個自變量即一進制線性回歸。一進制線性回歸方程：y=β0+β1x+ε (ε是随機誤差，滿足E(ε)=0，var(ε)=σ2)其中選擇β0，β1通常采用最小二乘法(用殘差平方和來刻畫所有觀察值與回歸直線的偏差程度)。
基本調用方法：
>> [b, bint, r, rint, s] =
regress(y, xdata); % 調用regress函數作一進制線性回歸
>> yhat = xdata*b; %
計算y的估計值
B = REGRESS(Y,X)
傳回值為線性模型Y = X*B的回歸系數向量
X ：n-by-p 矩陣，行對應于觀測值，列對應于預測變量
Y ：n-by-1 向量，觀測值的響應(即因變量，譯者注)
[B,BINT] =
REGRESS(Y,X)
BINT：B的95%的置信區間矩陣。Bint
置信區間不大，說明有效性較好；若含零點，說明結果無效。
[B,BINT,R] =
REGRESS(Y,X)
R：殘差向量(因變量的真實值減去估計值)
[B,BINT,R,RINT] =
REGRESS(Y,X)
RINT：傳回殘差的95%置信區間，它是一個2×n的矩陣，第1列為置信下限，第2列為置信上限。該矩陣可以用來診斷異常(即發現奇異觀測值，譯者注)。
如果第i組觀測的殘差的置信區間RINT(i，：)所定區間沒有包含0，則第i個殘差在預設的5%的顯著性水準比我們所預期的要大，這可說明第i個觀測值是個奇異點(即說明該點可能是錯誤而無意義的，如記錄錯誤等，譯者注)
[B,BINT,R,RINT,STATS] =
REGRESS(Y,X)
STATS：向量，STATS中的4個值分别為：R2(判定系數)，F(總模型的F測驗值)，P(總模型F的機率值P(F>Fz))，MSq(離回歸方差或誤差方差的估計值)。
判定系數(the Coefficient of the
Determination)R2：是判斷回歸模型拟合程度的一個名額，其取值範圍為[0,
1]；判定系數越大說明回歸模型的拟合程度越高，回歸方程越顯著。
F>F(1-α)(k,
n-k-1)時拒絕H0，F越大，說明回歸方程越顯著。
與F對應的機率P
MSq：由于最小二乘法中不求誤差方差σ2，其誤差平方和Msq定義為SSR/自由度，其中SSR為Regression Sum of Squares。
[B,BINT,R,RINT,STATS]=regress(Y,X,
ALPHA)
用100*(1-ALPHA)%的置信水準來計算BINT，用(100*ALPHA)%的顯著性水準來計算RINT。
PS. X應該包含一個全“1”的列，這樣則該模型包含常數項。利用它實作X=[ones(size(x,1))
x]
Fz(總模型的F測驗值)和p(總模型F的機率值P(F>Fz))是在模型有常數項的假設下計算的，如果模型沒有常數項，則計算得的F統計量和p值是不正确的。
若模型沒有常數項，則這個值可以為負值，這也表明這個模型對資料是不合适的。(即資料不适合用多元線性模型，譯者注)
如果X的列是線性相關的，則REGRESS将使B的元素中“0”的數量盡量多，以此獲得一個“基本解”，并且使B中元素“0”所對應的BINT元素為“0”。
REGRESS 将X或者Y中的NaNs當作缺失值處理，并且移除它們。
PPS.
rcoplot(r,rint)
這是個畫殘差的函數，紅色的表示超出期望值的資料。圓圈代表殘差的值，豎線代表置信區間的範圍。
PPPS.polyfit函數中第一個多項式系數p([p,S] =
polyfit(x,y,n))是按高次到低次排列的行向量；而regress中的b(b =
regress(y,X))是從低到高，而且是個列向量。若想得到與polyfit相同的多項式系數向量，需要将b順時針旋轉90度，即p=rot90(b,-1)。注意，僅僅将b進行轉置還不夠，還需要fliplr才可以。
>> ClimateData =
xlsread('examp01_01.xls'); % 從Excel檔案讀取資料
>> x =
ClimateData(:, 1); % 提取ClimateData的第1列，即年平均氣溫資料
>> y =
ClimateData(:, 5); % 提取ClimateData的第5列，即全年日照時數資料
>> xdata =
[ones(size(x, 1), 1), x]; % 在原始資料x的左邊加一列1，即模型包含常數項
>> [b, bint, r,
rint, s] = regress(y, xdata); % 調用regress函數作一進制線性回歸
>> yhat = xdata*b;
% 計算y的估計值
% 定義元胞數組，以元胞數組形式顯示系數的估計值和估計值的95%置信區間
>> head1 =
{'系數的估計值','估計值的95%置信下限','估計值的95%置信上限'};
>> [head1; num2cell([b, bint])]
ans =
'系數的估計值'
'估計值的95%置信下限'
'估計值的95%置信上限'
[ 3.1154e+003] [ 2.6592e+003] [ 3.5716e+003]
[ -76.9616] [ -105.9970] [ -47.9261]
% 定義元胞數組，以元胞數組形式顯示y的真實值、y的估計值、殘差和殘差的95%置信區間
>> head2 =
{'y的真實值','y的估計值','殘差','殘差的95%置信下限','殘差的95%置信上限'};
% 同時顯示y的真實值、y的估計值、殘差和殘差的95%置信區間
>> [head2; num2cell([y, yhat, r, rint])]
ans =
'y的真實值'
'y的估計值'
'殘差'
'殘差的95%置信下限'
'殘差的95%置信上限'
[2.3511e+003] [2.0379e+003] [ 313.1847] [ -461.6917]
[1.0881e+003]
[2.1654e+003] [2.0687e+003] [ 96.7001] [ -686.1726] [
879.5727]
[2.1677e+003] [1.9686e+003] [ 199.0501] [ -581.9400] [
980.0403]
[2.1746e+003] [2.2380e+003] [ -63.4154] [ -840.7773] [
713.9466]
[2.6478e+003] [2.4227e+003] [ 225.0768] [ -534.8155] [
984.9692]
[2.3609e+003] [2.4227e+003] [ -61.8232] [ -826.3056] [
702.6593]
[2.5336e+003] [2.5228e+003] [ 10.8268] [ -744.1662] [
765.8198]
[2.3592e+003] [2.6074e+003] [-248.2309] [ -987.0492] [
490.5873]
[1.5222e+003] [1.6916e+003] [-169.3882] [ -944.3372] [
605.5608]
[1.6803e+003] [1.7762e+003] [ -95.9459] [ -876.4774] [
684.5855]
[1.4729e+003] [1.6993e+003] [-226.3844] [ -999.5520] [
546.7832]
[1.8146e+003] [1.7762e+003] [ 38.3541] [ -742.9172] [
819.6253]
[1.5438e+003] [1.4992e+003] [ 44.6157] [ -719.2939] [
808.5253]
[ 2102] [1.6377e+003] [ 464.2849] [ -289.2774]
[1.2178e+003]
[1.8198e+003] [1.9610e+003] [-141.1537] [ -924.0249] [
641.7175]
[1.7472e+003] [1.8840e+003] [-136.7921] [ -919.1600] [
645.5757]
[1.9342e+003] [1.6839e+003] [ 250.3079] [ -520.9526]
[1.0216e+003]
[1.7422e+003] [1.6685e+003] [ 73.7003] [ -702.2379] [
849.6384]
[ 1616] [1.3299e+003] [ 286.1312] [ -451.5246]
[1.0238e+003]
[ 1614] [1.4453e+003] [ 168.6888] [ -587.4691] [
924.8467]
[1.6691e+003] [1.2606e+003] [ 408.4966] [ -311.0565]
[1.1280e+003]
[ 856.2000] [1.6531e+003] [-796.9074] [-1.5087e+003] [
-85.1229]
[ 935.6000] [1.8224e+003] [-886.8229] [-1.5907e+003] [
-182.9957]
[1.0148e+003] [1.9686e+003] [-953.8499] [-1.6466e+003] [
-261.0701]
[2.0386e+003] [1.9148e+003] [ 123.8232] [ -659.2486] [
906.8951]
[ 3181] [2.3612e+003] [ 819.8461] [ 118.1779]
[1.5215e+003]
[1.8936e+003] [1.9148e+003] [ -21.1768] [ -805.6671] [
763.3135]
[2.2141e+003] [2.2611e+003] [ -47.0039] [ -823.2940] [
729.2863]
[2.3647e+003] [2.6459e+003] [-281.2117] [-1.0132e+003] [
450.7301]
[2.5298e+003] [2.3150e+003] [ 214.8230] [ -553.8992] [
983.5453]
[2.8534e+003] [2.4612e+003] [ 392.1961] [ -353.8710]
[1.1383e+003]
% 定義元胞數組，以元胞數組形式顯示判定系數、F統計量的觀測值、檢驗的p值和誤差方差的估計值
>> head3 =
{'判定系數','F統計量的觀測值','檢驗的p值','誤差方差的估計值'};
>> [head3; num2cell(s)]
'判定系數'
'F統計量的觀測值'
'檢驗的p值'
'誤差方差的估計值'
[ 0.5033] [ 29.3884] [7.8739e-006] [ 1.4689e+005]
從輸出的結果看，常數項和回歸系數的估計值分别為3.1154e+003和-76.9616，進而可以寫出線性回歸方程為
y = 3.1154e+003 *
x -76.9616
p=7.8739e-006，即p
利用下面指令畫出原始資料散點與回歸直線圖，如圖1-2所示。
>> plot(x, y, 'k.',
'Markersize', 15) % 畫原始資料散點
>> hold on
>> plot(x, yhat,
'linewidth', 3) % 畫回歸直線
>>
xlabel('年平均氣溫(x)')
% 給X軸加标簽
>>
ylabel('全年日照時數(y)')
% 給Y軸加标簽
>>
legend('原始散點','回歸直線');
% 加标注框
一個有用的函數
function stats =
reglm(y,X,model,varnames)
%
多重線性回歸分析或廣義線性回歸分析
%
%
reglm(y,X)，産生線性回歸分析的方差分析表和參數估計結果，并以表格形式顯示在螢幕上.
參
%
數X是自變量觀測值矩陣，它是n行p列的矩陣.
y是因變量觀測值向量，它是n行1列的列向量.
%
%
stats = reglm(y,X)，還傳回一個包括了回歸分析的所有診斷統計量的結構體變量stats.
%
%
stats = reglm(y,X,model)，用可選的model參數來控制回歸模型的類型.
model是一個字元串，
%
其可用的字元串如下
%
'linear' 帶有常數項的線性模型(預設情況)
%
'interaction' 帶有常數項、線性項和交叉項的模型
%
'quadratic' 帶有常數項、線性項、交叉項和平方項的模型
%
'purequadratic' 帶有常數項、線性項和平方項的模型
%
%
stats = reglm(y,X,model,varnames)，用可選的varnames參數指定變量标簽.
varnames
%
可以是字元矩陣或字元串元胞數組，它的每行的字元或每個元胞的字元串是一個變量的标簽，它的行
%
數或元胞數應與X的列數相同.
預設情況下，用X1,X2,…作為變量标簽.
%
%
例:
% x =
[215 250 180 250 180 215 180 215 250 215 215
%
136.5 136.5 136.5 138.5 139.5 138.5 140.5 140.5 140.5 138.5
138.5]';
% y =
[6.2 7.5 4.8 5.1 4.6 4.6 2.8 3.1 4.3 4.9 4.1]';
%
reglm(y,x,'quadratic')
%
%
----------------------------------方差分析表----------------------------------%
方差來源自由度 平方和 均方 F值
p值
%
回歸 5.0000 15.0277 3.0055 7.6122 0.0219
%
殘差 5.0000 1.9742 0.3948
%
總計 10.0000 17.0018
%
%
均方根誤差(Root MSE) 0.6284 判定系數(R-Square) 0.8839
%
因變量均值(Dependent Mean) 4.7273 調整的判定系數(Adj R-Sq) 0.7678
%
%
-----------------------------------參數估計-----------------------------------
%
變量估計值 标準誤 t值
p值
%
常數項 30.9428 2011.1117 0.0154 0.9883
% X1
0.7040 0.6405 1.0992 0.3218
% X2
-0.8487 29.1537 -0.0291 0.9779
%
X1*X2 -0.0058 0.0044 -1.3132 0.2461
%
X1*X1 0.0003 0.0003 0.8384 0.4400
%
X2*X2 0.0052 0.1055 0.0492 0.9626
%
%
Copyright 2009 - 2010 xiezhh.
%
$Revision: 1.0.0.0 $ $Date: 2009/12/22 21:41:00 $
if
nargin < 2
error('至少需要兩個輸入參數');
end
p =
size(X,2); % X的列數，即變量個數
if
nargin < 3 || isempty(model)
model
= 'linear'; % model參數的預設值
end
%
生成變量标簽varnames
if
nargin < 4 || isempty(varnames)
varname1 =
strcat({'X'},num2str([1:p]'));
varnames = makevarnames(varname1,model); %
預設的變量标簽
else
if
ischar(varnames)
varname1 = cellstr(varnames);
elseif
iscell(varnames)
varname1 = varnames(:);
else
error('varnames 必須是字元矩陣或字元串元胞數組');
end
if
size(varname1,1) ~= p
error('變量标簽數與X的列數不一緻');
else
varnames = makevarnames(varname1,model); %
指定的變量标簽
end
end
ST =
regstats(y,X,model); % 調用regstats函數進行線性回歸分析，傳回結構體變量ST
f =
ST.fstat; % F檢驗相關結果
t =
ST.tstat; % t檢驗相關結果
%
顯示方差分析表
fprintf('n');
fprintf('------------------------------方差分析表------------------------------');
fprintf('n');
fprintf('%s%7sssss','方差來源','自由度','平方和','均方','F值','p值');
fprintf('n');
fmt =
'%s.4f.4f.4f.4f.4f';
fprintf(fmt,'回歸',f.dfr,f.ssr,f.ssr/f.dfr,f.f,f.pval);
fprintf('n');
fmt =
'%s.4f.4f.4f';
fprintf(fmt,'殘差',f.dfe,f.sse,f.sse/f.dfe);
fprintf('n');
fmt =
'%s.4f.4f';
fprintf(fmt,'總計',f.dfe+f.dfr,f.sse+f.ssr);
fprintf('n');
fprintf('n');
%
顯示判定系數等統計量
fmt =
'"s.4f%s.4f';
fprintf(fmt,'均方根誤差(Root
MSE)',sqrt(ST.mse),'判定系數(R-Square)',ST.rsquare);
fprintf('n');
fprintf(fmt,'因變量均值(Dependent
Mean)',mean(y),'調整的判定系數(Adj
R-Sq)',...
ST.adjrsquare);
fprintf('n');
fprintf('n');
%
顯示參數估計及t檢驗相關結果
fprintf('-------------------------------參數估計-------------------------------');
fprintf('n');
fprintf('%8sssss','變量','估計值','标準誤','t值','p值');
fprintf('n');
for i
= 1:size(t.beta,1)
if i
== 1
fmt =
'%8s .4f.4f.4f.4fn';
fprintf(fmt,'常數項',t.beta(i),t.se(i),t.t(i),t.pval(i));
else
fmt =
's .4f.4f.4f.4fn';
fprintf(fmt,varnames{i-1},t.beta(i),t.se(i),t.t(i),t.pval(i));
end
end
if
nargout == 1
stats
= ST; % 傳回一個包括了回歸分析的所有診斷統計量的結構體變量
end
%
-----------------子函數-----------------------
function varnames =
makevarnames(varname1,model)
%
生成指定模型的變量标簽
p =
size(varname1,1);
varname2 = [];
for i
= 1:p-1
varname2 =
[varname2;strcat(varname1(i),'*',varname1(i+1:end))];
end
varname3 =
strcat(varname1,'*',varname1);
switch
model
case
'linear'
varnames = varname1;
case
'interaction'
varnames = [varname1;varname2];
case
'quadratic'
varnames =
[varname1;varname2;varname3];
case
'purequadratic'
varnames = [varname1;varname3];
end