[转载]【MATLAB】MATLAB 线性拟合小结 ——&nb

原文地址：【MATLAB】MATLAB 线性拟合小结 —— REGRESS多元线性回归（用最小二乘估计法）作者：vera吧噗

 MATLAB 线性拟合小结 —— REGRESS多元线性回归（用最小二乘估计法）

 

http://blog.sina.com.cn/s/blog_491b86bf0100ny17.html

http://www.labfans.com/bbs/t6794/

http://wenku.baidu.com/view/0a0ea0de941ea76e59fa0418.html?re=view 

 

在回归分析中，如果有两个或两个以上的自变量，就称为多元回归。如果只有一个自变量即一元线性回归。一元线性回归方程：y=β₀+β₁x+ε    (ε是随机误差，满足E(ε)=0，var(ε)=σ²)其中选择β₀，β₁通常采用最小二乘法（用残差平方和来刻画所有观察值与回归直线的偏差程度）。

 

 基本调用方法：

>> [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个值分别为：R²（判定系数），F（总模型的F测验值），P（总模型F的概率值P(F>Fz)），MSq（离回归方差或误差方差的估计值）。

判定系数(the Coefficient of the Determination)R²：是判断回归模型拟合程度的一个指标，其取值范围为[0, 1]；判定系数越大说明回归模型的拟合程度越高，回归方程越显著。

F>F_(1-α)(k, n-k-1)时拒绝H₀，F越大，说明回归方程越显著。

与F对应的概率P<α时拒绝H₀，回归模型成立。

MSq：由于最小二乘法中不求误差方差σ²，其误差平方和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才可以。

 

 

举例：http://feng2002226.blog.163.com/blog/static/4703365020122711468498/

>> 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<α=0.05，回归模型成立。

 

利用下面命令画出原始数据散点与回归直线图，如图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