广义线性回归拟合教程和源码

原地址：http://www.matlabsky.com/thread-9145-1-1.html

为了使得回归分析的结果更加直观，我调用regstats函数，编写了一个更为实用的函数：reglm，代码如下（代码中有调用方法和例子）。下面再给出一个例子： 【例】这里有这样一组数据，涉及三个变量：p，t 和z，要拟合出 z = f(p,t) 的关系式（非线性的）。   解：调用自编的reglm函数进行二次拟合，命令如下： >> z = [9.73875 20.75 36.59875 13.5725 29.6325 50.93875 18.97875 36.59875 80.13875 20.75125 38.22125 90.925 22.055 44.58 104.7725]; >> [p,t] = meshgrid([0.8 1 1.2],[60:60:300]); >> stats = reglm(z(:),[p(:), t(:)],'quadratic',{'p','t'}); >> [pnew,tnew] = meshgrid(linspace(0.8,1.2,20),linspace(60,300,20)); >> pp = pnew(:); >> tt = tnew(:); >> zhat = [ones(400,1) pp tt pp.*tt pp.^2 tt.^2]*stats.beta; >> mesh(pnew,tnew,reshape(zhat,[20,20])); >> hold on >> plot3(p(:),t(:),z(:),'k*') 拟合结果： ------------------------------------方差分析表------------------------------------ 方差来源    自由度            平方和             均方             F值          p值 回归       5.0000       11548.9147        2309.7829         93.4739      0.0000 残差       9.0000         222.3942          24.7105 总计      14.0000       11771.3089        均方根误差(Root MSE)         4.9710           判定系数(R-Square)    0.9811 因变量均值(Dependent Mean)        41.2168        调整的判定系数(Adj R-Sq)    0.9706 -----------------------------------参数估计-----------------------------------       变量               估计值            标准误             t值          p值      常数项            242.6188          69.0439           3.5140      0.0066          p           -513.7781         137.3777          -3.7399      0.0046          t             -0.3637           0.1212          -3.0002      0.0150        p*t              0.6022           0.0926           6.5010      0.0001        p*p            272.2625          68.0677           3.9999      0.0031        t*t             -0.0003           0.0002          -1.1946      0.2628

reglm函数的代码如下：

1. function stats = reglm(y,X,model,varnames)
   
2. % 多重线性回归分析或广义线性回归分析
   
3. %
   
4. %   reglm(y,X)，产生线性回归分析的方差分析表和参数估计结果，并以表格形式显示在屏幕上. 参
   
5. %   数X是自变量观测值矩阵，它是n行p列的矩阵. y是因变量观测值向量，它是n行1列的列向量.
   
6. %
   
7. %   stats = reglm(y,X)，还返回一个包括了回归分析的所有诊断统计量的结构体变量stats.
   
8. %
   
9. %   stats = reglm(y,X,model)，用可选的model参数来控制回归模型的类型. model是一个字符串，
   
10. %   其可用的字符串如下
    
11. %       'linear'          带有常数项的线性模型（默认情况）
    
12. %       'interaction'     带有常数项、线性项和交叉项的模型
    
13. %       'quadratic'       带有常数项、线性项、交叉项和平方项的模型
    
14. %       'purequadratic'   带有常数项、线性项和平方项的模型
    
15. %
    
16. %   stats = reglm(y,X,model,varnames)，用可选的varnames参数指定变量标签. varnames 
    
17. %   可以是字符矩阵或字符串元胞数组，它的每行的字符或每个元胞的字符串是一个变量的标签，它的行
    
18. %   数或元胞数应与X的列数相同. 默认情况下，用X1,X2,…作为变量标签.
    
19. %
    
20. %   例:
    
21. %   x = [215 250 180 250 180 215 180 215 250 215 215
    
22. %       136.5 136.5 136.5 138.5 139.5 138.5 140.5 140.5 140.5 138.5 138.5]';
    
23. %   y = [6.2 7.5 4.8 5.1 4.6 4.6 2.8 3.1 4.3 4.9 4.1]';
    
24. %   reglm(y,x,'quadratic')
    
25. %
    
26. %   ----------------------------------方差分析表----------------------------------
    
27. %   方差来源    自由度          平方和           均方             F值          p值
    
28. %   回归       5.0000        15.0277         3.0055          7.6122      0.0219
    
29. %   残差       5.0000         1.9742         0.3948
    
30. %   总计      10.0000        17.0018
    
31. %
    
32. %         均方根误差(Root MSE)     0.6284             判定系数(R-Square)    0.8839
    
33. %   因变量均值(Dependent Mean)     4.7273        调整的判定系数(Adj R-Sq)    0.7678
    
34. %
    
35. %   -----------------------------------参数估计-----------------------------------
    
36. %      变量               估计值            标准误             t值          p值
    
37. %     常数项             30.9428        2011.1117           0.0154      0.9883
    
38. %        X1              0.7040           0.6405           1.0992      0.3218
    
39. %        X2             -0.8487          29.1537          -0.0291      0.9779
    
40. %     X1*X2             -0.0058           0.0044          -1.3132      0.2461
    
41. %     X1*X1              0.0003           0.0003           0.8384      0.4400
    
42. %     X2*X2              0.0052           0.1055           0.0492      0.9626
    
43. %
    
44. %   Copyright 2009 - 2010 xiezhh. 
    
45. %   $Revision: 1.0.0.0 $  $Date: 2009/12/22 21:41:00 $
    
46. 
    
47. if nargin < 2
    
48.    error('至少需要两个输入参数');
    
49. end
    
50. p = size(X,2);    % X的列数，即变量个数
    
51. if nargin < 3 || isempty(model)
    
52.    model = 'linear';    % model参数的默认值
    
53. end
    
54. 
    
55. % 生成变量标签varnames
    
56. if nargin < 4 || isempty(varnames)
    
57.     varname1 = strcat({'X'},num2str([1:p]'));
    
58.     varnames = makevarnames(varname1,model);    % 默认的变量标签
    
59. else
    
60.     if ischar(varnames)
    
61.         varname1 = cellstr(varnames);
    
62.     elseif iscell(varnames)
    
63.         varname1 = varnames(:);
    
64.     else
    
65.         error('varnames 必须是字符矩阵或字符串元胞数组');
    
66.     end
    
67.     if size(varname1,1) ~= p
    
68.         error('变量标签数与X的列数不一致');
    
69.     else
    
70.         varnames = makevarnames(varname1,model);    % 指定的变量标签
    
71.     end
    
72. end
    
73. 
    
74. ST = regstats(y,X,model);    % 调用regstats函数进行线性回归分析，返回结构体变量ST
    
75. f = ST.fstat;    % F检验相关结果
    
76. t = ST.tstat;    % t检验相关结果
    
77. 
    
78. % 显示方差分析表
    
79. fprintf('\n');
    
80. fprintf('------------------------------方差分析表------------------------------');
    
81. fprintf('\n');
    
82. fprintf('%s%7s%15s%15s%15s%12s','方差来源','自由度','平方和','均方','F值','p值');
    
83. fprintf('\n');
    
84. fmt = '%s%13.4f%17.4f%17.4f%16.4f%12.4f';
    
85. fprintf(fmt,'回归',f.dfr,f.ssr,f.ssr/f.dfr,f.f,f.pval);
    
86. fprintf('\n');
    
87. fmt = '%s%13.4f%17.4f%17.4f';
    
88. fprintf(fmt,'残差',f.dfe,f.sse,f.sse/f.dfe);
    
89. fprintf('\n');
    
90. fmt = '%s%13.4f%17.4f';
    
91. fprintf(fmt,'总计',f.dfe+f.dfr,f.sse+f.ssr);
    
92. fprintf('\n');
    
93. fprintf('\n');
    
94. 
    
95. % 显示判定系数等统计量
    
96. fmt = '%22s%15.4f%25s%10.4f';
    
97. fprintf(fmt,'均方根误差(Root MSE)',sqrt(ST.mse),'判定系数(R-Square)',ST.rsquare);
    
98. fprintf('\n');
    
99. fprintf(fmt,'因变量均值(Dependent Mean)',mean(y),'调整的判定系数(Adj R-Sq)',...
    
100.         ST.adjrsquare); 
     
101. fprintf('\n');
     
102. fprintf('\n');
     
103. 
     
104. % 显示参数估计及t检验相关结果
     
105. fprintf('-------------------------------参数估计-------------------------------');
     
106. fprintf('\n');
     
107. fprintf('%8s%18s%15s%15s%12s','变量','估计值','标准误','t值','p值');
     
108. fprintf('\n');
     
109. for i = 1:size(t.beta,1)
     
110.     if i == 1
     
111.         fmt = '%8s%20.4f%17.4f%17.4f%12.4f\n';
     
112.         fprintf(fmt,'常数项',t.beta(i),t.se(i),t.t(i),t.pval(i));
     
113.     else
     
114.         fmt = '%10s%20.4f%17.4f%17.4f%12.4f\n';
     
115.         fprintf(fmt,varnames{i-1},t.beta(i),t.se(i),t.t(i),t.pval(i));
     
116.     end
     
117. end
     
118. 
     
119. if nargout == 1
     
120.     stats = ST;    % 返回一个包括了回归分析的所有诊断统计量的结构体变量
     
121. end
     
122. 
     
123. % -----------------子函数-----------------------
     
124. function varnames = makevarnames(varname1,model)
     
125. % 生成指定模型的变量标签
     
126. p = size(varname1,1);
     
127. varname2 = [];
     
128. for i = 1:p-1
     
129.     varname2 = [varname2;strcat(varname1(i),'*',varname1(i+1:end))];
     
130. end
     
131. varname3 = strcat(varname1,'*',varname1);
     
132. switch model
     
133.     case 'linear'
     
134.         varnames = varname1;
     
135.     case 'interaction'
     
136.         varnames = [varname1;varname2];
     
137.     case 'quadratic'
     
138.         varnames = [varname1;varname2;varname3];
     
139.     case 'purequadratic'
     
140.         varnames = [varname1;varname3];
     
141. end

1. function stats = reglm(y,X,model,varnames)
   
2. % 多重线性回归分析或广义线性回归分析
   
3. %
   
4. %   reglm(y,X)，产生线性回归分析的方差分析表和参数估计结果，并以表格形式显示在屏幕上. 参
   
5. %   数X是自变量观测值矩阵，它是n行p列的矩阵. y是因变量观测值向量，它是n行1列的列向量.
   
6. %
   
7. %   stats = reglm(y,X)，还返回一个包括了回归分析的所有诊断统计量的结构体变量stats.
   
8. %
   
9. %   stats = reglm(y,X,model)，用可选的model参数来控制回归模型的类型. model是一个字符串，
   
10. %   其可用的字符串如下
    
11. %       'linear'          带有常数项的线性模型（默认情况）
    
12. %       'interaction'     带有常数项、线性项和交叉项的模型
    
13. %       'quadratic'       带有常数项、线性项、交叉项和平方项的模型
    
14. %       'purequadratic'   带有常数项、线性项和平方项的模型
    
15. %
    
16. %   stats = reglm(y,X,model,varnames)，用可选的varnames参数指定变量标签. varnames 
    
17. %   可以是字符矩阵或字符串元胞数组，它的每行的字符或每个元胞的字符串是一个变量的标签，它的行
    
18. %   数或元胞数应与X的列数相同. 默认情况下，用X1,X2,…作为变量标签.
    
19. %
    
20. %   例:
    
21. %   x = [215 250 180 250 180 215 180 215 250 215 215
    
22. %       136.5 136.5 136.5 138.5 139.5 138.5 140.5 140.5 140.5 138.5 138.5]';
    
23. %   y = [6.2 7.5 4.8 5.1 4.6 4.6 2.8 3.1 4.3 4.9 4.1]';
    
24. %   reglm(y,x,'quadratic')
    
25. %
    
26. %   ----------------------------------方差分析表----------------------------------
    
27. %   方差来源    自由度          平方和           均方             F值          p值
    
28. %   回归       5.0000        15.0277         3.0055          7.6122      0.0219
    
29. %   残差       5.0000         1.9742         0.3948
    
30. %   总计      10.0000        17.0018
    
31. %
    
32. %         均方根误差(Root MSE)     0.6284             判定系数(R-Square)    0.8839
    
33. %   因变量均值(Dependent Mean)     4.7273        调整的判定系数(Adj R-Sq)    0.7678
    
34. %
    
35. %   -----------------------------------参数估计-----------------------------------
    
36. %      变量               估计值            标准误             t值          p值
    
37. %     常数项             30.9428        2011.1117           0.0154      0.9883
    
38. %        X1              0.7040           0.6405           1.0992      0.3218
    
39. %        X2             -0.8487          29.1537          -0.0291      0.9779
    
40. %     X1*X2             -0.0058           0.0044          -1.3132      0.2461
    
41. %     X1*X1              0.0003           0.0003           0.8384      0.4400
    
42. %     X2*X2              0.0052           0.1055           0.0492      0.9626
    
43. %
    
44. %   Copyright 2009 - 2010 xiezhh. 
    
45. %   $Revision: 1.0.0.0 $  $Date: 2009/12/22 21:41:00 $
    
46. 
    
47. if nargin < 2
    
48.    error('至少需要两个输入参数');
    
49. end
    
50. p = size(X,2);    % X的列数，即变量个数
    
51. if nargin < 3 || isempty(model)
    
52.    model = 'linear';    % model参数的默认值
    
53. end
    
54. 
    
55. % 生成变量标签varnames
    
56. if nargin < 4 || isempty(varnames)
    
57.     varname1 = strcat({'X'},num2str([1:p]'));
    
58.     varnames = makevarnames(varname1,model);    % 默认的变量标签
    
59. else
    
60.     if ischar(varnames)
    
61.         varname1 = cellstr(varnames);
    
62.     elseif iscell(varnames)
    
63.         varname1 = varnames(:);
    
64.     else
    
65.         error('varnames 必须是字符矩阵或字符串元胞数组');
    
66.     end
    
67.     if size(varname1,1) ~= p
    
68.         error('变量标签数与X的列数不一致');
    
69.     else
    
70.         varnames = makevarnames(varname1,model);    % 指定的变量标签
    
71.     end
    
72. end
    
73. 
    
74. ST = regstats(y,X,model);    % 调用regstats函数进行线性回归分析，返回结构体变量ST
    
75. f = ST.fstat;    % F检验相关结果
    
76. t = ST.tstat;    % t检验相关结果
    
77. 
    
78. % 显示方差分析表
    
79. fprintf('\n');
    
80. fprintf('------------------------------方差分析表------------------------------');
    
81. fprintf('\n');
    
82. fprintf('%s%7s%15s%15s%15s%12s','方差来源','自由度','平方和','均方','F值','p值');
    
83. fprintf('\n');
    
84. fmt = '%s%13.4f%17.4f%17.4f%16.4f%12.4f';
    
85. fprintf(fmt,'回归',f.dfr,f.ssr,f.ssr/f.dfr,f.f,f.pval);
    
86. fprintf('\n');
    
87. fmt = '%s%13.4f%17.4f%17.4f';
    
88. fprintf(fmt,'残差',f.dfe,f.sse,f.sse/f.dfe);
    
89. fprintf('\n');
    
90. fmt = '%s%13.4f%17.4f';
    
91. fprintf(fmt,'总计',f.dfe+f.dfr,f.sse+f.ssr);
    
92. fprintf('\n');
    
93. fprintf('\n');
    
94. 
    
95. % 显示判定系数等统计量
    
96. fmt = '%22s%15.4f%25s%10.4f';
    
97. fprintf(fmt,'均方根误差(Root MSE)',sqrt(ST.mse),'判定系数(R-Square)',ST.rsquare);
    
98. fprintf('\n');
    
99. fprintf(fmt,'因变量均值(Dependent Mean)',mean(y),'调整的判定系数(Adj R-Sq)',...
    
100.         ST.adjrsquare); 
     
101. fprintf('\n');
     
102. fprintf('\n');
     
103. 
     
104. % 显示参数估计及t检验相关结果
     
105. fprintf('-------------------------------参数估计-------------------------------');
     
106. fprintf('\n');
     
107. fprintf('%8s%18s%15s%15s%12s','变量','估计值','标准误','t值','p值');
     
108. fprintf('\n');
     
109. for i = 1:size(t.beta,1)
     
110.     if i == 1
     
111.         fmt = '%8s%20.4f%17.4f%17.4f%12.4f\n';
     
112.         fprintf(fmt,'常数项',t.beta(i),t.se(i),t.t(i),t.pval(i));
     
113.     else
     
114.         fmt = '%10s%20.4f%17.4f%17.4f%12.4f\n';
     
115.         fprintf(fmt,varnames{i-1},t.beta(i),t.se(i),t.t(i),t.pval(i));
     
116.     end
     
117. end
     
118. 
     
119. if nargout == 1
     
120.     stats = ST;    % 返回一个包括了回归分析的所有诊断统计量的结构体变量
     
121. end
     
122. 
     
123. % -----------------子函数-----------------------
     
124. function varnames = makevarnames(varname1,model)
     
125. % 生成指定模型的变量标签
     
126. p = size(varname1,1);
     
127. varname2 = [];
     
128. for i = 1:p-1
     
129.     varname2 = [varname2;strcat(varname1(i),'*',varname1(i+1:end))];
     
130. end
     
131. varname3 = strcat(varname1,'*',varname1);
     
132. switch model
     
133.     case 'linear'
     
134.         varnames = varname1;
     
135.     case 'interaction'
     
136.         varnames = [varname1;varname2];
     
137.     case 'quadratic'
     
138.         varnames = [varname1;varname2;varname3];
     
139.     case 'purequadratic'
     
140.         varnames = [varname1;varname3];
     
141. end