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