前面博客中已经讲过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
> y
> slt
> 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.0