前面博客中已经讲过MATLAB中常用的命令拟合 polyfit() , lsqcurvefit() ,nlinfit()  和 cftool等，这里简单介绍简单的回归分析的MATLAB和R语言实现。

例15.1

测得某种物质在不同温度下吸附另一种物质的重量如下表

  $\begin{array}{c|ccccccccc} x_i(^oC)&1.5& 1.8&2.4&3.0 &3.5 &3.9 &4.4 &4.8 &5.0 \\\hline y_i(mg)&4.8& 5.7&7.0&8.3&10.9&12.4&13.1&13.6&15.1 \end{array}$

由所给定的样本观测值，如果画出散点图，可以看出 $9$ 个点近乎在一条直线上。因此，我们可以假设（严格地讲，应该先检验这假设）吸附量 Y 与温度 $x$ 具有线性关系：$Y = a+b,x+varepsilon$

，求 $Y$ 对 $x$ 的线性回归。

按照书上公式，在MATLAB中输入：

x = [1.5 , 1.8,  2.4, 3.0, 3.5, 3.9, 4.4, 4.8,5.0];

y = [4.8, 5.7, 7.0, 8.3, 10.9, 12.4, 13.1, 13.6,15.1];

sx = sum(x); sy = sum(y);  sxy = sum(x.*y); sxx = sum(x.^2);

fprintf('sx = %6.2f,   sy = %6.2f, sxy = %6.2f,  sxx = %6.2fn', sx,sy,sxy,sxx)

mx =1/9* sum(x); my = 1/9*sum(y);

fprintf('mx = %6.4f,   my = %6.4fn', mx, my)

lxy = sxy - 9*mx*my; lxx = sxx - 9*mx^2;

bhat = lxy/lxx; ahat = my - bhat*mx;

fprintf('ahat = %6.4f,   bhat = %6.4fn', ahat, bhat)

即可得到

sx =  30.30,   sy =  90.90, sxy = 344.09,  sxx = 115.11

mx = 3.3667,   my = 10.1000

ahat = 0.3187,   bhat = 2.9053

R语言的程序代码为：

> x<-c(1.5,1.8,2.4,3.0,3.5,3.9,4.4,4.8,5.0)

> y<-c(4.8,5.7,7.0,8.3,10.9,12.4,13.1,13.6,15.1)

> slt<-lm(y~x)

> summary(slt)

运行后得到

Call:

lm(formula = y ~ x)

Residuals:

Min      1Q  Median      3Q     Max

-0.7347 -0.2915  0.1233  0.2546  0.7505

Coefficients:

Estimate Std. Error t value Pr(>|t|)    

(Intercept)  0.3187  0.5151  0.619 0.556    

x        2.9053   0.1440  20.170 1.84e-07 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5213 on 7 degrees of freedom

Multiple R-squared:  0.9831,    Adjusted R-squared:  0.9807

F-statistic: 406.8 on 1 and 7 DF,  p-value: 1.844e-07

显示参数的置信区间

>      confint(slt)

                2.5 %   97.5 %

(Intercept) -0.8994414 1.536795

x            2.5647360 3.245951

>      x0<-data.frame(x=32)

>      lm.pred<-predict(slt,x0,interval = "prediction",level = 0.95)

>      lm.pred

   fit      lwr      upr

1  93.28967 83.45075 103.1286

利用xtable包可以输出为 LaTeX 代码(见：使用 xtable 将 R 表格输出 LaTeX url{http://elegantlatex.org/2015/03/31/xtable-export-r-to-latex/})，

源代码为

% latex table generated in R 3.2.3 by xtable 1.8-2 package

% Tue Mar 08 14:47:45 2016

begin{table}[ht]

centering

caption{linear model regression}

begin{tabular}{rrrrr}

toprule

& Estimate & Std. Error & t value & Pr($>$$|$t$|$) \

midrule

(Intercept) & 0.3187 & 0.5151 & 0.62 & 0.5558 \

x & 2.9053 & 0.1440 & 20.17 & 0.0000 \

bottomrule

end{tabular}

end{table}

编译结果如下

注意：

lm( )为线性模型程序包，它可以完成参数估计、假设检验。本程序中将 lm(y~x) 结果存储于slt, 再用summary(), confint( ) 和 ~predict( ) 等函数显示或调用它。

例15.5

在某种产品表面进行腐蚀刻线试验，得到腐蚀深度~$Y$ 与腐蚀时间~$x$对应的一组数据：

$\begin{array}{c|ccccccccccc} x_i(s)& 5& 10&15&20 &30 &40 &50 &60 &70&90&120 \\\hline y_i(\mu g)&6&10&10&13&16&17&19&23&25&29&46 \end{array}$

（1）预测腐蚀时间为 $75s$ 时，腐蚀深度的范围（$1-alpha = 95%$ ）；

（2）若要求腐蚀深度在 $10sim20mu m$ 之间，问腐蚀时间应如何控制？

x = [5 10 15 20 30 40 50 60 70 90 120];

y = [6 10 10 13 16 17 19 23 25 29 46];

 n = length(x); lxx = sum(x.^2) - n*mean(x)^2; lxy = sum(x.*y)-n*mean(x)*mean(y);

lyy = sum(y.^2) - n*mean(y)^2; bhat = lxy/lxx; ahat = mean(y) - bhat*mean(x);

fprintf('ahat = %6.4f, bhat = %6.4fn', ahat, bhat)

disp(strcat('回归方程 y = ',num2str(ahat),' + ',' ',num2str(bhat),' x'))

%(1) x0 = 75

x0 = 75; y0hat = ahat + bhat*x0;sig2hat = 1/(n-2)*(lyy-lxy^2/lxx);sighat = sqrt(sig2hat)

deltax0 = sighat*tinv(0.975,n-2)*sqrt(1+1/n+(x0-mean(x))^2/lxx)

得到

ahat = 5.3444, bhat = 0.3043

回归方程 y =5.3444 +0.30434 x

sighat =

   2.2356

deltax0 =

   5.4315

x0 = 75; y0hat = ahat + bhat*x0;sig2hat = 1/(n-2)*(lyy-lxy^2/lxx);sighat = sqrt(sig2hat)

deltax0 = sighat*tinv(0.975,n-2)*sqrt(1+1/n+(x0-mean(x))^2/lxx)

fprintf('腐蚀时间为75s时，腐蚀深度Y0的预测区间为[%6.4f, %6.4f]n',y0hat-deltax0,y0hat+deltax0)

%(2) 当要求腐蚀深度在10～20 mu m 之间时，近似地有

x1 = 1/bhat*(10+sighat*norminv(0.975)-ahat); x2 = 1/bhat*(20-sighat*norminv(0.975)-ahat);

fprintf('即腐蚀时间应控制在%6.4f s到%6.4f s之间，就能得到10~20 之间的腐蚀深度。n',x1,x2)

得到

sighat =

   2.2356

deltax0 =

   5.4315

腐蚀时间为75s时，腐蚀深度Y0的预测区间为[22.7381, 33.6011]

即腐蚀时间应控制在29.6950 s到33.7584 s之间，就能得到10~20 之间的腐蚀深度。

例15.6

出钢时所用盛钢水的钢包，由于钢水对耐火材料的侵蚀，容积不断扩大，我们希望找出使用次数 $x$ 与增大容积 $Y$之间的关系。试验数据列于下表：

  $\begin{array}{c|cccccccc} \hline x_i&2& 3&4&5&6 &7 &8 &9\\ y_i&6.42& 8.20&9.58&9.50&9.70&10.00&9.93&9.99\\ \hline \hline x_i&10&11&12&13&14&15&16& \\ y_i&10.49&10.59&10.60&10.80&10.60&10.90&11.76&\\ \hline \end{array}$

x = 2:16; y = [6.42,8.20,9.58,9.50,9.70,10.00,9.93,...

   9.99,10.49,10.59,10.60,10.80,10.60,10.90,11.76];

z = 1./y'; t = 1./x'; n = length(x);

[abhat] = [ones(n,1) t]z;

fprintf('ahat = %6.4f, bhat = %6.4fn', abhat(1), abhat(2))

disp(strcat('回归方程 z = ',num2str(abhat(1)),' + ',' ',num2str(abhat(2)),' t'))

得到

ahat = 0.0812, bhat = 0.1349

回归方程 z =0.081193 +0.13491 t

R 语言程序如下

>  x<-2:16; y<-c(6.42,8.20,9.58,9.50,9.70,10.00,9.93,9.99,10.49,10.59,10.60,10.80,10.60,10.90,11.76)

>   x1<- 1/x

>   y1<-1/y

>   slt<-lm(y1~x1)

>   summary(slt)

得到

Call:

lm(formula = y1 ~ x1)

Residuals:

Min         1Q     Median         3Q        Max

-0.0105351 -0.0017476  0.0009717  0.0022768  0.0071176

Coefficients:

Estimate Std. Error t value Pr(>|t|)    

(Intercept) 0.081193 0.001915 42.40 2.52e-15 ***

x1          0.134905 0.009701 13.91 3.50e-09 ***

---

Signif. codes:0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.004408 on 13 degrees of freedom

Multiple R-squared:  0.937,     Adjusted R-squared:  0.9322

F-statistic: 193.4 on 1 and 13 DF,  p-value: 3.5e-09

>   confint(slt)

            2.5 %     97.5 %

(Intercept) 0.07705633 0.08532936

x1          0.11394772 0.15586327

非线性回归nlinfit()

  nlinfit( ) 函数可以做一般的非线性回归，其调用格式为

[beta,r,J]=nlinfit(x,y,’model’, beta0)

其中beta表示估计出的回归系数，r表示残差，J表示Jacobian矩阵，x, y表示输入数据，x、y分别为矩阵和n维列向量, 对一元非线性回归, x为n维列向量, model表示是事先定义的非线性函数, beta0表示回归系数的初值.

   相关命令还有产生交互界面的nlintool(x,y,’model’, beta0,alpha)、非线性模型置信区间预测的[Y,DELTA]=nlpredci(’model’, x,beta,r,J) ，表示nlinfit 或nlintool所得的回归函数在x处的预测值Y及预测值的显著性为1-alpha的置信区间Y±DELTA.

在 例15.6中，也可以直接用MATLAB解算过程如下，输入

>>x=2:16;

y=[6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76];

Fit1 = @(beta,x)(x./(beta(1)*x+beta(2))) ;% 非线性模型1

[beta,r ,J]=nlinfit(x,y,Fit1,[0.1,0.1])

得到

beta =

    0.084463      0.11524

r =

    -0.61815

    0.061728

   …………

    -0.14924

J =    -49.535      -24.768

     -66.231      -22.077

…………

     -119.01      -7.4382

即回归模型为 $y = dfrac{x}{0.084463x+0.11524}$，下面预测并作图，输入

>>  [yh,delta]=nlpredci(Fit1,x,beta,r ,J); plot(x,y,'k+',x,yh,'r')

得到

图15.13 散点图及非线性拟合曲线

若用例15.7模型，则继续输入

>>  Fit2=@(beta,x) (beta(1)*exp(beta(2)./x)); [beta,r ,J]=nlinfit(x,y,Fit2,[11,8])

得到

beta =       11.604      -1.0641

r =   -0.39589

    0.061455

…………

   -0.097033

J =

     0.58739       3.4079

     0.70138       2.7128

…………

     0.93566      0.67856

即拟合模型为 $y = 11.604 e^{-1.641/x}$.

>>   [yhh,delta]=nlpredci(Fit2,x,beta,r ,J); plot(x,y,'k+',x,yhh,'r')

可得到类似于图15.13 的散点图和新的非线性拟合曲线。