爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 7370|回复: 16

[程序设计] 非线性拟合matlab

[复制链接]

新浪微博达人勋

发表于 2017-4-15 15:42:22 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
提供给需要的人。
%非线性回归
%一个自变量的形式
%指数模型(I)线性回归
x=[0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5];
y=[26300 25100 19900 15500 11500 9800 5200 4600        3200 2300 1700 1200        900        700        600        500];
plot(x,y,'r.');                             %绘制散点图
xlabel('Distance');                         %横轴标签(到城市中心的距离)
ylabel('Average density');                  %纵轴标签(人口平均密度)
hold on                                     %保持图形
X=[ones(length(y),1),x'];                   %自变量矩阵
Y=log(y');                                  %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=exp(B(1));                                %模型常数还原
b=-B(2);                                    %回归系数
f=a*exp(-b*x);                              %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,R2,s                                    %输出主要结果

%指数模型(I)的非线性拟合
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       yhat = b1*exp(b2*x);
%以上模型通过编辑窗口保存在Matlab的work文件夹中

%指数模型(I)的非线性拟合
x=[0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5];
y=[26300 25100 19900 15500 11500 9800 5200 4600        3200 2300 1700 1200        900        700        600        500];
plot(x,y,'r.');                             %绘制散点图
xlabel('Distance');                         %横轴标签(到城市中心的距离)
ylabel('Average density');                  %纵轴标签(人口平均密度)
hold on                                     %保持图形
beta0=[0 0];                                %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
a=B(1);                                     %模型常数
b=-B(2);                                    %回归系数
f=a*exp(-b*x);                              %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,s                                       %输出主要结果

%对数模型的线性回归
x=[100.6 103.5 231.3 120.4 230.4 234.3 162.7 236.3 158.7 145.2 207 203.3 433.5 372.9 525.3 629.2 963.4 608.8 876.7 832 703.1 872.6 2196.2 2422.4 2230.5        1117.1 2558.6 1190 1750.2 3710 6050        5760 4460 6618.2 6272.8        3840 6926.3        3580        5817.2        6610];
y=[2.6 4 6.7 8.9 10.2 12.6 14.6        18 21 22 26.1 28 31        32 34 36.4 38.6        40.8 43        44.5 46.4 48 50        52.7 54.3 59 60.1 62 64.4 67 69        70 72 74 76        78 80.7        82 86.4        88];
plot(x,y,'r.');                             %绘制散点图
xlabel('Per capita income');                %横轴标签(人均收入)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
X=[ones(length(y),1),log(x')];              %自变量矩阵
Y=y';                                       %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数
f=a+b*log(x);                               %模型表达
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
f=a+b*log(sort(x));                         %模型表达
plot(sort(x),f,'b-');                       %添加趋势线
hold off                                    %绘图结束
a,b,R2,s                                    %输出主要结果
x=sort(x);                                  %自变量排序

%对数模型的非线性拟合
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       yhat = b1 + b2*log(x);
%以上模型通过编辑窗口保存在Matlab的work文件夹中
x=[100.6 103.5 231.3 120.4 230.4 234.3 162.7 236.3 158.7 145.2 207 203.3 433.5 372.9 525.3 629.2 963.4 608.8 876.7 832 703.1 872.6 2196.2 2422.4 2230.5        1117.1 2558.6 1190 1750.2 3710 6050        5760 4460 6618.2 6272.8        3840 6926.3        3580        5817.2        6610];
y=[2.6 4 6.7 8.9 10.2 12.6 14.6        18 21 22 26.1 28 31        32 34 36.4 38.6        40.8 43        44.5 46.4 48 50        52.7 54.3 59 60.1 62 64.4 67 69        70 72 74 76        78 80.7        82 86.4        88];
plot(x,y,'r.');                             %绘制散点图
xlabel('Per capita income');                %横轴标签(人均收入)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
beta0=[0 0];                                %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数
f=a+b*log(x);                               %模型表达
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
f=a+b*log(sort(x));                         %模型表达
plot(sort(x),f,'b-');                       %添加趋势线
hold off                                    %绘图结束
a,b,s                                       %输出主要结果

%幂指数模型的线性回归
x=[3.48        3.69 1.43 2.05 3.05        4.19 2.43 2.4 2.72 2.99        4.78 1.33 1.67 3.14        2.04 1.77 0.59 0.69        0.5        0.69 0.63 0.58 0.86        0.41 1.23];
y=[38.83 43.92 9.14        16.66 36.16        38.66 17.74        19.46 23 29.75 51.19 6.6 9.04 34.27        17.61 13.37        2.04 2.22 1.46 1.92        1.86 1.69 3.31 1.13        6.74];
plot(x,y,'r.');                             %绘制散点图
xlabel('Perimeter');                        %横轴标签(周长)
ylabel('Area');                             %纵轴标签(面积)
hold on                                     %保持图形
X=[ones(length(y),1),log(x')];              %自变量矩阵
Y=log(y');                                  %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=exp(B(1));                                %模型常数还原
b=B(2);                                     %回归系数
f=a*x.^b;                                   %模型表达
s=sqrt(sumsqr(y-f)/(length(y)-2));          %计算标准误差
f=a*(sort(x)).^b;                           %模型表达
plot(sort(x),f,'b-');                       %添加趋势线
hold off                                    %绘图结束
a,b,R2,s                                    %输出主要结果

%幂指数模型的非线性拟合
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       yhat = b1*x.^b2;
%以上模型通过编辑窗口保存在Matlab的work文件夹中
x=[3.48        3.69 1.43 2.05 3.05        4.19 2.43 2.4 2.72 2.99        4.78 1.33 1.67 3.14        2.04 1.77 0.59 0.69        0.5        0.69 0.63 0.58 0.86        0.41 1.23];
y=[38.83 43.92 9.14        16.66 36.16        38.66 17.74        19.46 23 29.75 51.19 6.6 9.04 34.27        17.61 13.37        2.04 2.22 1.46 1.92        1.86 1.69 3.31 1.13        6.74];
plot(x,y,'r.');                             %绘制散点图
xlabel('Perimeter');                        %横轴标签(周长)
ylabel('Area');                             %纵轴标签(面积)
hold on                                     %保持图形
beta0=[0 0];                                %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数
f=a*x.^b;                                   %模型表达
s=sqrt(sumsqr(y-f)/(length(y)-2));          %计算标准误差
f=a*(sort(x)).^b;                           %模型表达
plot(sort(x),f,'b-');                       %添加趋势线
hold off                                    %绘图结束
a,b,s                                       %输出主要结果

%双曲线模型的线性回归
t=[1650        1700 1750 1800 1850        1900 1950 1960 1965        1970 1975 1980 1985        1990 1995 2000];
y=[0.510 0.625 0.710 0.910 1.130 1.600 2.525 3.307 3.354 3.696 4.066 4.432 4.822 5.318 5.660 6.060];
x=t-t(1);                                   %将年份转换为时序
figure(1)                                   %创建第一个图形窗口
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Population');                       %纵轴标签(人口)
hold on                                     %保持图形
X=[ones(length(y),1),x'];                   %自变量矩阵
Y=(y.^(-1))';                               %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数还原
f=(a+b*x).^(-1);                            %模型表达
plot(x,f,'b-');                             %添加趋势线
hold on                                     %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,R2,s                                    %输出主要结果

%绘制散点与直线的拟合图
figure(2)                                   %创建第二个图形窗口
plot(x,y.^(-1),'r.');                       %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Population Reciprocal');            %纵轴标签(人口倒数)
hold on                                     %保持图形
f=a+b*x;                                    %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束

%双曲线模型的非线性拟合
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       yhat = (b1+b2*x).^(-1);
%以上模型通过编辑窗口保存在Matlab的work文件夹中
t=[1650        1700 1750 1800 1850        1900 1950 1960 1965        1970 1975 1980 1985        1990 1995 2000];
y=[0.510 0.625 0.710 0.910 1.130 1.600 2.525 3.307 3.354 3.696 4.066 4.432 4.822 5.318 5.660 6.060];
x=t-t(1);                                   %将年份转换为时序
figure(1)                                   %创建第一个图形窗口
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Population');                       %纵轴标签(人口)
hold on                                     %保持图形
beta0=[1 1];                                %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数
x=sort(x);                                  %自变量排序
f=(a+b*x).^(-1);                            %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,s                                       %输出主要结果

%绘制散点与直线的拟合图
figure(2)                                   %创建第二个图形窗口
plot(x,y.^(-1),'r.');                       %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Population Reciprocal');            %纵轴标签(人口倒数)
hold on                                     %保持图形
f=a+b*x;                                    %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束

%Logistic模型的线性回归
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
x=t-t(1);                                   %将年份转换为时序
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
X=[ones(length(y),1),x'];                   %自变量矩阵
Y=(log(y.^(-1)-1))';                        %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=exp(B(1));                                %模型常数还原
b=-B(2);                                    %回归系数
f=(1+a*exp(-b*x)).^(-1);                    %模型表达
plot(x,f,'b-');                             %添加趋势线
hold on                                     %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,R2,s                                    %输出主要结果

%Logistic模型的非线性拟合(二参数)
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       yhat = (1+b1*exp(b2*x)).^(-1);
%以上模型通过编辑窗口保存在Matlab的work文件夹中
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
x=t-t(1);                                   %将年份转换为时序
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
beta0=[1 0];                                %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
a=B(1);                                     %模型常数
b=-B(2);                                    %回归系数
f=(1+a*exp(-b*x)).^(-1);                    %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,s                                       %输出主要结果

%指数模型(II)线性回归
x=[1.380 1.444 1.401 1.526 1.674 1.446 1.474 1.366 1.544 1.742 1.307 1.565 1.432 1.463 1.472 1.471 1.502 1.494 1.340 1.305 1.422 1.278 1.466 1.538 1.456 1.356 1.441 1.494 1.366 1.426 1.493];
y=[0.350 0.127 0.276 0.204 0.080 0.198 0.158 0.346 0.118 0.111 0.450 0.172 0.248 0.239         0.209 0.195 0.241 0.152 0.424 0.318 0.160 0.421 0.247 0.131 0.234 0.260 0.168 0.135 0.288 0.247 0.249];
plot(x,y,'r.');                             %绘制散点图
xlabel('Dimension');                        %横轴标签(边界维数)
ylabel('Compactness');                      %纵轴标签(紧凑度)
hold on                                     %保持图形
X=[ones(length(y),1),(x.^(-1))'];           %自变量矩阵
Y=log(y');                                  %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=exp(B(1));                                %模型常数还原
b=B(2);                                     %回归系数
f=a*exp(b*x.^(-1));                         %模型表达
s=sqrt(sumsqr(y-f)/(length(y)-2));          %计算标准误差
f=a*exp(b*(sort(x)).^(-1));                 %模型表达
plot(sort(x),f,'b-');                       %添加趋势线
hold off                                    %绘图结束
a,b,R2,s                                    %输出主要结果

%指数模型(II)的非线性拟合
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       yhat = b1*exp(b2*x.^(-1));
%以上模型通过编辑窗口保存在Matlab的work文件夹中
x=[1.380 1.444 1.401 1.526 1.674 1.446 1.474 1.366 1.544 1.742 1.307 1.565 1.432 1.463 1.472 1.471 1.502 1.494 1.340 1.305 1.422 1.278 1.466 1.538 1.456 1.356 1.441 1.494 1.366 1.426 1.493];
y=[0.350 0.127 0.276 0.204 0.080 0.198 0.158 0.346 0.118 0.111 0.450 0.172 0.248 0.239         0.209 0.195 0.241 0.152 0.424 0.318 0.160 0.421 0.247 0.131 0.234 0.260 0.168 0.135 0.288 0.247 0.249];
plot(x,y,'r.');                             %绘制散点图
xlabel('Dimension');                        %横轴标签(边界维数)
ylabel('Compactness');                      %纵轴标签(紧凑度)
hold on                                     %保持图形
beta0=[0 0];                                %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数
f=a*exp(b*x.^(-1));                         %模型表达
s=sqrt(sumsqr(y-f)/(length(y)-2));          %计算标准误差
f=a*exp(b*(sort(x)).^(-1));                 %模型表达
plot(sort(x),f,'b-');                       %添加趋势线
hold off                                    %绘图结束
a,b,s                                       %输出主要结果

%Logistic模型化为指数模型的线性回归
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
x=t-t(1);                                   %将年份转换为时序
y=y./(1-y);                                 %城市化水平化为城乡人口比
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Urban-Rural Ratio');                %纵轴标签(城乡人口比)
hold on                                     %保持图形
X=[ones(length(y),1),x'];                   %自变量矩阵
Y=(log(y))';                                %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=exp(B(1));                                %模型常数还原
b=-B(2);                                    %回归系数
f=a*exp(-b*x);                              %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,R2,s                                    %输出主要结果

%Logistic模型化为指数模型的非线性拟合
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       yhat = b1*exp(b2*x);
%以上模型通过编辑窗口保存在Matlab的work文件夹中
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
x=t-t(1);                                   %将年份转换为时序
y=y./(1-y);                                 %城市化水平化为城乡人口比
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Urban-Rural Ratio');                %纵轴标签(城乡人口比)
hold on                                     %保持图形
beta0=[0 0];                                %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数
f=a*exp(b*x);                               %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,s                                       %输出主要结果

%一个自变量化为多个自变量的形式
%多项式线性回归
x=[3.4 1.8 4.6 2.3 3.1 5.5 0.7 3.0 2.6 4.3 2.1 1.1 6.1 4.8 3.8];
y=[26.2 17.8 31.3 23.1 27.5 36.0 14.1 22.3 19.6 31.3 24.0 17.3 43.2 36.4 26.1];
plot(x,y,'r.');                             %绘制散点图
xlabel('Distance');                         %横轴标签(到消防队的距离)
ylabel('Loss');                             %纵轴标签(火灾引起的损失)
hold on                                     %保持图形
X=[ones(length(y),1),x',x'.^2];             %自变量矩阵
Y=y';                                       %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
f=B(1)+B(2)*x+B(3)*x.^2;                    %模型表达
s=sqrt(sumsqr(y-f)/(length(f)-3));          %计算标准误差
x=sort(x);                                  %自变量排序
f=B(1)+B(2)*x+B(3)*x.^2;                    %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
B,R2,s                                      %输出主要结果

%多项式拟合
x=[3.4 1.8 4.6 2.3 3.1 5.5 0.7 3.0 2.6 4.3 2.1 1.1 6.1 4.8 3.8];
y=[26.2 17.8 31.3 23.1 27.5 36.0 14.1 22.3 19.6 31.3 24.0 17.3 43.2 36.4 26.1];
plot(x,y,'r.');                             %绘制散点图
xlabel('Distance');                         %横轴标签(到消防队的距离)
ylabel('Loss');                             %纵轴标签(火灾引起的损失)
hold on                                     %保持图形
N=2;                                        %确定多项式阶次
[P,S]=polyfit(x,y,N);                       %多项式拟合
lx=linspace(min(x),max(x));                 %限定趋势线的长度范围
z=polyval(P,lx);                            %计算分割点上多项式的函数值
plot(lx,z,'b-');                            %将趋势线添加到散点图中
hold off                                    %绘图结束
e=y-polyval(P,x);                           %计算残差
s=sqrt(sumsqr(e)/(length(y)-N-1));          %计算标准误差
P,s                                         %输出主要结果

%二次指数模型
x=[0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5];
y=[26300 25100 19900 15500 11500 9800 5200 4600        3200 2300 1700 1200        900        700        600        500];
plot(x,y,'r.');                             %绘制散点图
xlabel('Distance');                         %横轴标签(到城市中心的距离)
ylabel('Average density');                  %纵轴标签(人口平均密度)
hold on                                     %保持图形
P=polyfit(x,log(y),2);                      %多项式拟合
lx=linspace(min(x),max(x));                 %限定趋势线的长度范围
z=exp(polyval(P,lx));                       %计算分割点上多项式的函数值
plot(lx,z,'b-');                            %将趋势线添加到散点图中
hold off                                    %绘图结束
e=y-polyval(P,x);                           %计算残差
s=sqrt(sumsqr(e)/(length(y)-N-1));          %计算标准误差
P,s                                         %输出主要结果

%Logistic模型的非线性自回归(三参数)1——利用自回归计算
%绘制抛物线图
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
n=length(y);                                %计算样本路径长度
x=t-t(1);                                   %将年份转换为时序
figure(1)                                   %创造第一个图形窗口
plot(y(1:n-1),y(2:n),'r.');                 %抛物线散点图
hold on                                     %保持图形
P=polyfit(y(1:n-1),y(2:n),2);               %多项式拟合
lx=linspace(min(y(1:n-1)),max(y(1:n-1)));   %限定趋势线的长度范围
z=polyval(P,lx);                            %计算分割点上多项式的函数值
plot(lx,z,'b-');                            %将趋势线添加到散点图中
hold off                                    %绘图结束
Pv=polyval(P,y(1:n-1));                     %计算抛物线的预测值
G=sumsqr(Pv)/sumsqr(y(2:n));                %计算抛物线拟合优度(常数不为0)

%参数的初步估计(承接前面的程序)
X=[y(1:n-1);y(1:n-1).^2]';                  %自变量矩阵
Y=y(2:n)';                                  %因变量向量
B=inv(X'*X)*X'*Y;                           %最小二乘计算
G=sumsqr(X*B)/sumsqr(Y);                    %计算抛物线拟合优度(常数为0)
c=(1-B(1))/B(2);                            %估计承载量参数
figure(2)                                   %创造第二个图形窗口
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
a=c/y(1)-1;                                 %估计常数值
b=(B(1)-1)/10;                              %计算内生增长率
f=c*(1+a*exp(-b*x)).^(-1);                  %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,c,s                                     %输出主要结果

%第一种解决办法(承接初步估计程序)
figure(3)                                   %创造第三个图形窗口
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
at=exp(log(c*y.^(-1)-1)+b*x);               %自变量矩阵
a=mean(at);                                 %计算系数值
f=c*(1+a*exp(-b*x)).^(-1);                  %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(n-2));                  %计算标准误差
a,b,c,s                                     %输出主要结果

%第二种解决办法(承接初步估计程序)
X=[ones(length(y),1),x'];                   %自变量矩阵
Y=(log(c*y.^(-1)'-1));                      %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=exp(B(1));                                %模型常数还原
b=-B(2);                                    %回归系数
f=c*(1+a*exp(-b*x)).^(-1);                  %模型表达
figure(3)                                   %创造第三个图形窗口
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(n-2));                  %计算标准误差
a,b,c,R2,s                                  %输出主要结果

%Logistic模型的非线性自回归(三参数)2——利用差分自回归计算
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
n=length(y);                                %计算样本路径长度
x=t-t(1);                                   %将年份转换为时序
dt=mean(t(2:n)-t(1:n-1));                   %计算时间间隔
X=[y(1:n-1);y(1:n-1).^2]';                  %自变量矩阵
Y=[(y(2:n)-y(1:n-1))/dt]';                  %因变量向量
A=inv(X'*X)*X'*Y;                           %最小二乘计算
G=sumsqr(X*A)/sumsqr(Y);                    %计算抛物线拟合优度(常数为0)
c=-A(1)/A(2);                               %估计承载量参数
figure(1)                                   %创建第一个图形窗口
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
X=[ones(length(y),1),x'];                   %自变量矩阵
Y=(log(c*y.^(-1)'-1));                      %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=exp(B(1));                                %模型常数还原
b=-B(2);                                    %回归系数
f=c*(1+a*exp(-b*x)).^(-1);                  %模型表达
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
figure(2)                                   %创建第二个图形窗口
x=y(1:n-1);                                 %城市化水平序列去尾
y=(y(2:n)-y(1:n-1))/dt;                     %城市化水平增长率
plot(x,y,'r.');                             %绘制散点图
xlabel('Percent Urban');                    %横轴标签(城市化水平)
ylabel('Growth Rate');                      %纵轴标签(城市化水平增长率)
hold on                                     %保持图形
f=A(1)*x+A(2)*x.^2;                         %抛物线模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
a,b,c,R2,s                                  %输出主要结果

%Logistic模型的非线性拟合(三参数)
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       b3 = beta(3);
       yhat = b1*(1+b2*exp(b3*x)).^(-1);
%以上模型通过编辑窗口保存在Matlab的work文件夹中
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
x=t-t(1);                                   %将年份转换为时序
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
beta0=[0 1 0];                              %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
a=B(2);                                     %模型常数
b=-B(3);                                    %回归系数
c=B(1);                                     %承载量参数
f=c*(1+a*exp(-b*x)).^(-1);                  %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,c,s                                     %输出主要结果

%对于本例,初始值可以设为[0 1 0]或者[1 0 0]或者[1 1 0],但不宜设为[0 0 0]或者[0 0 1]或者[0 1 1]或者[1 0 1]或者[1 1 1]。否则计算过程不能收敛到正确的位置。

%Gamma模型
%Gamma模型的线性回归
x=[0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5];
y=[26300 25100 19900 15500 11500 9800 5200 4600        3200 2300 1700 1200        900        700        600        500];
plot(x,y,'r.');                             %绘制散点图
xlabel('Distance');                         %横轴标签(到城市中心的距离)
ylabel('Average density');                  %纵轴标签(人口平均密度)
hold on                                     %保持图形
X=[ones(length(y),1),log(x'),x'];           %自变量矩阵
Y=log(y');                                  %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=B(1);                                     %模型常数
b=-B(2);                                    %回归系数1
c=-B(3);                                    %回归系数2
k=exp(a);                                   %模型常数还原
f=k*(x.^(-b)).*exp(-c*x);                   %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
D=2-b;                                      %估计城市形态的分维
x0=1/c;                                     %估计城市特征半径
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
k,D,x0,R2,s                                 %输出主要结果

%Gamma模型的非线性拟合
%myfun.m
function  yhat = myfun(beta, x)
       b1 = beta(1);
       b2 = beta(2);
       b3 = beta(3);
       yhat = b1*(x.^b2).*exp(b3*x);
%以上模型通过编辑窗口保存在Matlab的work文件夹中
x=[0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5];
y=[26300 25100 19900 15500 11500 9800 5200 4600        3200 2300 1700 1200        900        700        600        500];
plot(x,y,'r.');                             %绘制散点图
xlabel('Distance');                         %横轴标签(到城市中心的距离)
ylabel('Average density');                  %纵轴标签(人口平均密度)
hold on                                     %保持图形
beta0=[0 0 1];                              %设定迭代初始值
O=statset('MaxIter',200);                   %设定最大迭代次数
[B,E,J]=nlinfit(x,y,'myfun',beta0,O);       %非线性拟合
k=B(1);                                     %模型常数
D=2+B(2);                                   %估计分维值
x0=-1/B(3);                                 %估计分维值
f=B(1)*(x.^B(2)).*exp(B(3)*x);              %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
k,D,x0,s                                    %输出主要结果

%对于本例,初始值可以设为[0 0 0]、[0 0 1]、[0 1 0]、[0 1 1]、[1 1 0],但不宜设为[1 0 0]、[1 0 1]、[1 1 1]。否则计算过程不能收敛到正确的位置。

%淄博旅游业Cobb-Douglas生产函数的线性回归
x=[350.97 387.43 337.00 340.00 420.00 444.85 508.46 583.61 595.84 771.54
0.52 0.79 0.70 0.98 1.22 1.45 2.10 2.68 2.30 3.00];
y=[9.19 12.26 10.38 11.15 18.31 22.00 28.20 34.06 34.05 48.08];
[m,n]=size(x);                              %计算自变量数组的行列数
X=[ones(length(y),1),log(x')];              %自变量矩阵
Y=log(y');                                  %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %提前拟合优度
DW=sumsqr(E(2:n)-E(1:n-1))/sumsqr(E);       %计算Durbin-Watson统计量
rcoplot(E,Eint);                            %绘制残差序列的置信区间图
a=exp(B(1));                                %模型常数还原
b1=B(2);                                    %回归系数1
b2=B(3);                                    %回归系数2
k=exp(a);                                   %模型常数还原
f=a*(x(1,:).^b1).*x(2,:).^b2;               %模型表达
s=sqrt(sumsqr(y-f)/(length(f)-m-1));        %计算非线性模型的标准误差
a,b1,b2,R2,DW,s                             %输出主要结果

%淄博旅游业问题的线性回归
x=[350.97 387.43 337.00 340.00 420.00 444.85 508.46 583.61 595.84 771.54
0.52 0.79 0.70 0.98 1.22 1.45 2.10 2.68 2.30 3.00];
y=[9.19 12.26 10.38 11.15 18.31 22.00 28.20 34.06 34.05 48.08];
[m,n]=size(x);                              %计算自变量数组的行列数
X=[ones(length(y),1),x'];                   %自变量矩阵
Y=y';                                       %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %提前拟合优度
DW=sumsqr(E(2:n)-E(1:n-1))/sumsqr(E);       %计算Durbin-Watson统计量
a=B(1);                                     %模型常数
b1=B(2);                                    %回归系数1
b2=B(3);                                    %回归系数2
k=exp(a);                                   %模型常数还原
s=Stats(4)^0.5;                             %标准误差
rcoplot(E,Eint);                            %绘制残差序列的置信区间图
a,b1,b2,R2,DW,s                             %输出主要结果

%广义模型拟合
%最大积雪深度与灌溉面积之间的关系
x=[15.2        10.4 21.2 18.6 26.4        23.4 13.5 16.7 24 19.1];
y=[28.6        19.3 40.5 35.6 48.9        45 29.2        34.1 46.7 37.4];
scatter(x,y,'r.');                          %绘制散点图
xlabel('最大积雪深度x');                     %添加横轴标签
ylabel('灌溉面积y');                         %添加纵轴标签
hold on                                     %保持图形
[B, Dev, Stats]=glmfit(x ,y);               %广义线性拟合
x=sort(x);                                  %自变量排序
[yhat,dl,du]=glmval(B,x,'identity',Stats);  %计算预测值
plot(x',yhat,'b-');                         %添加趋势线
hold on                                     %保持图形
plot(x',yhat+du,'g--');                     %添加预测误差上限
hold on                                     %保持图形
plot(x',yhat-dl,'g--');                     %添加预测误差下限
hold off                                    %绘图结束
B,Dev,Stats                                 %输出主要结果

%工农业发展对交通运输业的影响
x1=[57.82 58.05 59.15 63.83 65.36 67.26 66.92 67.79 75.65 80.57 79.02 80.52 86.88 95.48 109.71 126.5 138.89 160.56];
x2=[27.05 28.89 33.02 35.23        24.94 32.95 30.35 38.70 47.99 54.18 58.73 59.85 64.57 70.97 81.54 94.01 103.23 119.33];
x3=[14.54 16.83 12.26 12.87        11.65 12.87 10.80 10.93 14.71 17.56 20.32 18.67 25.34 25.06 29.69 43.86 48.90 60.98];
y=[3.09 3.40 3.88 3.90 3.22 3.76 3.59 4.03 4.34 4.65 4.78 5.04 5.59 6.01 7.03 10.03 10.83 12.9];
x=[x1;x2;x3]';                              %自变量向量合并为矩阵
y=y';                                       %因变量向量转置
[B, Dev, Stats]=glmfit(x ,y);               %广义线性拟合
[yhat,dl,du]=glmval(B,x,'identity',Stats);  %计算预测值
B,Dev,Stats                                 %输出主要结果
%常数为0的广义线性拟合采用如下语句代替上面的相应语句
[B, Dev, Stats]=glmfit(x ,y, [], 'identity',[],[],[],'off');  

%指数模型(I)广义线性模型拟合(因变量取对数)
x=[0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5];
y=[26300 25100 19900 15500 11500 9800 5200 4600        3200 2300 1700 1200        900        700        600        500];
scatter(x,y,'r.');                          %绘制散点图
xlabel('Distance');                         %横轴标签(到城市中心的距离)
ylabel('Average density');                  %纵轴标签(人口平均密度)
hold on                                     %保持图形
[B,Dev,Stats]=glmfit(x,log(y));             %广义线性模型拟合
a=exp(B(1));                                %模型常数还原
b=-B(2);                                    %回归系数
f=a*exp(-b*x);                              %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,s                                       %输出主要结果

%指数模型(I)广义线性模型拟合(调用连接'log')
x=[0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5];
y=[26300 25100 19900 15500 11500 9800 5200 4600        3200 2300 1700 1200        900        700        600        500];
scatter(x,y,'r.');                          %绘制散点图
xlabel('Distance');                         %横轴标签(到城市中心的距离)
ylabel('Average density');                  %纵轴标签(人口平均密度)
hold on                                     %保持图形
[B,Dev,Stats]=glmfit(x,y,[],'log');         %广义线性模型拟合
[yhat,dl,du]=glmval(B,x,'log',Stats);       %计算预测值
a=exp(B(1));                                %模型常数还原
b=-B(2);                                    %回归系数
plot(x,yhat,'b-');                          %添加趋势线
hold on                                     %保持图形
plot(x',yhat+du,'g--');                     %添加预测误差上限
hold on                                     %保持图形
plot(x',yhat-dl,'g--');                     %添加预测误差下限
hold off                                    %绘图结束
s=sqrt(sumsqr(y'-yhat)/(length(y)-2));      %计算标准误差
a,b,s                                       %输出主要结果

%对数模型的广义线性模型拟合
x=[100.6 103.5 231.3 120.4 230.4 234.3 162.7 236.3 158.7 145.2 207 203.3 433.5 372.9 525.3 629.2 963.4 608.8 876.7 832 703.1 872.6 2196.2 2422.4 2230.5        1117.1 2558.6 1190 1750.2 3710 6050        5760 4460 6618.2 6272.8        3840 6926.3        3580        5817.2        6610];
y=[2.6 4 6.7 8.9 10.2 12.6 14.6        18 21 22 26.1 28 31        32 34 36.4 38.6        40.8 43        44.5 46.4 48 50        52.7 54.3 59 60.1 62 64.4 67 69        70 72 74 76        78 80.7        82 86.4        88];
plot(x,y,'r.');                             %绘制散点图
xlabel('Per capita income');                %横轴标签(人均收入)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
B=glmfit(log(x),y);                         %广义线性模型拟合
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数
f=a+b*log(x);                               %模型表达
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
x=sort(x);                                  %自变量排序
f=a+b*log(x);                               %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
a,b,s                                       %输出主要结果

%幂指数模型的广义线性拟合(自变量和因变量取对数)
x=[3.48        3.69 1.43 2.05 3.05        4.19 2.43 2.4 2.72 2.99        4.78 1.33 1.67 3.14        2.04 1.77 0.59 0.69        0.5        0.69 0.63 0.58 0.86        0.41 1.23];
y=[38.83 43.92 9.14        16.66 36.16        38.66 17.74        19.46 23 29.75 51.19 6.6 9.04 34.27        17.61 13.37        2.04 2.22 1.46 1.92        1.86 1.69 3.31 1.13        6.74];
plot(x,y,'r.');                             %绘制散点图
xlabel('Perimeter');                        %横轴标签(周长)
ylabel('Area');                             %纵轴标签(面积)
hold on                                     %保持图形
B=glmfit(log(x),log(y));                    %广义线性模型拟合
a=exp(B(1));                                %模型常数还原
b=B(2);                                     %回归系数
f=a*x.^b;                                   %模型表达
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
x=sort(x);                                  %自变量排序
f=a*x.^b;                                   %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
a,b,s                                       %输出主要结果

%双曲线模型的广义线性模型拟合(因变量取倒数)
clear
t=[1650        1700 1750 1800 1850        1900 1950 1960 1965        1970 1975 1980 1985        1990 1995 2000];
y=[0.510 0.625 0.710 0.910 1.130 1.600 2.525 3.307 3.354 3.696 4.066 4.432 4.822 5.318 5.660 6.060];
x=t-t(1);                                   %将年份转换为时序
figure(1)                                   %创建第一个图形窗口
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Population');                       %纵轴标签(人口)
hold on                                     %保持图形
[B,Dev,Stats]=glmfit(x,y.^(-1));            %广义线性模型拟合
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数还原
f=(a+b*x).^(-1);                            %模型表达
plot(x,f,'b-');                             %添加趋势线
hold on                                     %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,s                                       %输出主要结果

%绘制散点与直线的拟合图
figure(2)                                   %创建第二个图形窗口
plot(x,y.^(-1),'r.');                       %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Population Reciprocal');            %纵轴标签(人口倒数)
hold on                                     %保持图形
f=a+b*x;                                    %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束

%双曲线模型的广义线性模型拟合(调用连接'reciprocal')
clear
t=[1650        1700 1750 1800 1850        1900 1950 1960 1965        1970 1975 1980 1985        1990 1995 2000];
y=[0.510 0.625 0.710 0.910 1.130 1.600 2.525 3.307 3.354 3.696 4.066 4.432 4.822 5.318 5.660 6.060];
x=t-t(1);                                   %将年份转换为时序
figure(1)                                   %创建第一个图形窗口
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Population');                       %纵轴标签(人口)
hold on                                     %保持图形
[B,Dev,Stats]=glmfit(x,y,[],'reciprocal');  %广义线性模型拟合
yhat=glmval(B,x,'reciprocal');              %计算预测值
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数还原
plot(x',yhat,'b-');                         %添加趋势线
hold on                                     %绘图结束
s=sqrt(sumsqr(y'-yhat)/(length(y)-2));      %计算标准误差
a,b,s                                       %输出主要结果

%绘制散点与直线的拟合图
figure(2)                                   %创建第二个图形窗口
plot(x,y.^(-1),'r.');                       %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Population Reciprocal');            %纵轴标签(人口倒数)
hold on                                     %保持图形
plot(x',yhat.^(-1),'b-');                   %添加趋势线
hold off                                    %绘图结束

%Logistic模型的广义线性模型拟合(线性化处理)
clear
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
x=t-t(1);                                   %将年份转换为时序
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
[B,Dev,Stats]=glmfit(x,log(y.^(-1)-1));     %广义线性模型拟合
a=exp(B(1));                                %模型常数
b=-B(2);                                    %回归系数
f=(1+a*exp(-b*x)).^(-1);                    %模型表达
plot(x,f,'b-');                             %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,s                                       %输出主要结果

%Logistic模型的广义线性模型拟合(调用连接'logit')
clear
t=[1790        1800 1810 1820 1830        1840 1850 1860 1870        1880 1890 1900 1910        1920 1930 1940 1950        1960];
y=[0.0513 0.0607 0.0726        0.0719 0.0877 0.1081 0.1541        0.1977 0.2568 0.2815 0.351 0.3965 0.4561 0.5117        0.5614 0.5652 0.5956 0.6305];
x=t-t(1);                                   %将年份转换为时序
plot(x,y,'r.');                             %绘制散点图
xlabel('Time');                             %横轴标签(时序)
ylabel('Percent Urban');                    %纵轴标签(城市化水平)
hold on                                     %保持图形
[B,Dev,Stats]=glmfit(x,y,[],'logit');       %广义线性模型拟合
yhat=glmval(B,x,'logit');                   %计算预测值
a=1/exp(B(1));                              %模型常数
b=B(2);                                     %回归系数
f=(1+a*exp(-b*x)).^(-1);                    %模型表达
plot(x',yhat,'b-');                         %添加趋势线
hold off                                    %绘图结束
s=sqrt(sumsqr(y'-yhat)/(length(y)-2));      %计算标准误差
a,b,s                                       %输出主要结果

%其他
%双曲线模型的线性回归(标准双曲线函数)
t=[1650        1700 1750 1800 1850        1900 1950 1960 1965        1970 1975 1980 1985        1990 1995 2000];
y=[0.510 0.625 0.710 0.910 1.130 1.600 2.525 3.307 3.354 3.696 4.066 4.432 4.822 5.318 5.660 6.060];
x=t-t(1);                                   %将年份转换为时序
figure(1)                                   %创建第一个图形窗口
plot(t,y,'r.');                             %绘制散点图
xlabel('Year');                             %横轴标签(年份)
ylabel('Population');                       %纵轴标签(人口)
hold on                                     %保持图形
X=[ones(length(y),1),t.^(-1)'];             %自变量矩阵
Y=(y.^(-1))';                               %因变量向量
[B,Bint,E,Eint,Stats]=regress(Y,X);         %回归分析
R2=Stats(1);                                %拟合优度
a=B(1);                                     %模型常数
b=B(2);                                     %回归系数还原
f=(a+b*t.^(-1)).^(-1);                      %模型表达
plot(t,f,'b-');                             %添加趋势线
hold on                                     %绘图结束
s=sqrt(sumsqr(y-f)/(length(f)-2));          %计算标准误差
a,b,R2,s                                    %输出主要结果
%绘制散点与直线的拟合图
figure(2)                                   %创建第二个图形窗口
plot(t.^(-1),y.^(-1),'r.');                 %绘制散点图
xlabel('Year Reciprocal');                  %横轴标签(年份倒数)
ylabel('Population Reciprocal');            %纵轴标签(人口倒数)
hold on                                     %保持图形
f=a+b*t.^(-1);                              %模型表达
plot(t.^(-1),f,'b-');                       %添加趋势线
hold off                                    %绘图结束

%Logistic回归
x=[18        21        23        23        28        31        36        42        46        48        55        56        58        18        20        25        27        28        30        32        33        33        38        41        45        48        52        56
850        1200        850        950        1200        850        1500        1000        950        1200        1800        2100        1800        850        1000        1200        1300        1500        950        1000        1800        1000        1200        1500        1800        1000        1500        1800
1        1        1        1        1        1        1        1        1        1        1        1        1        0        0        0        0        0        0        0        0        0        0        0        0        0        0        0]';
y=[0        0        1        1        1        0        1        1        1        0        1        1        1        0        0        0        0        0        1        0        0        0        0        0        1        0        1        1]';
[B, Dev, Stats]=glmfit(x ,y,'probit', 'logit')  %常数为0的广义线性拟合

%最大积雪深度与灌溉面积之间的关系的回归统计
x=[15.2        10.4 21.2 18.6 26.4        23.4 13.5 16.7 24 19.1];
y=[28.6        19.3 40.5 35.6 48.9        45 29.2        34.1 46.7 37.4];
scatter(x,y,'r.')
[B,Dev,Stats]=glmfit(x,y,[],'reciprocal');  %广义线性模型拟合
yhat=glmval(B,x,'reciprocal');              %计算预测值
r=y'-yhat;                                  %计算残差
regstats(y,x,'linear');                     %回归统计
scatter(yhat,r);                            %散点图
xlabel('Fitted Values');                    %横轴标签
ylabel('Residuals');                        %纵轴标签

评分

参与人数 1金钱 +2 贡献 +1 收起 理由
wxzy1314 + 2 + 1 很给力!

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-4-15 15:42:45 | 显示全部楼层
自己顶一下
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-4-15 15:43:02 | 显示全部楼层
111111111111111111
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-4-15 15:55:50 | 显示全部楼层
{:eb511:}{:eb511:}
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2017-12-18 12:21:40 | 显示全部楼层
赞,太感谢楼主了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-12-18 13:45:29 | 显示全部楼层
大赞楼主,非常实用
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-12-18 14:55:36 | 显示全部楼层
赞一个,很全面
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-12-18 15:46:49 | 显示全部楼层
学习了,感谢楼主。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-1-27 16:21:10 | 显示全部楼层
感谢楼主分享
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-3-29 11:04:16 | 显示全部楼层
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表