关于matlab最小二乘法的误差估计的方法

recbio

Matlab让线性最小二乘拟合问题简化为只用一个backslash operator。

比如  Y = f(x) = A1 + A2 * x + A3 * X ^2 + A4 * X ^3 方程，假设

x = rand( 100 , 1);
y =  0.7 + 33 * x + 7 * x.^2 +  2.7 * x.^3 + randn (100,1) * 0.05;

到这一步，我们假设了 A1 ~ A4 分别是 0.7, 33, 7, 2.7；x, y 都是 100个数的一维矩阵；并且在y上加了一定的误差。

如何通过拟合得到A呢？
只要假设一个矩阵 X 包含所有的拟合参数对自变量的一阶导数，如果是线性问题，参数的一阶导数就是参数对应的自变量这一项。
具体的：

Y 对 A1 的一阶导数是 1，Y 对 A2 的一阶导数是 x， Y 对 A3 的一阶导数是 x ^2，Y 对 A4 的一阶导数是  x ^3，于是我们构建如下的X：

X(:,1)=ones(100 , 1); X(:,2)=x; X(:,3)=x.^2; X(:,4)=x.^3;
并求得参数

a = X \ y                     。。。。。。。。。。（这一行就是最关键的“最小二乘拟合”）

执行结果如下（在Matlab中，如果命令不用;结尾，比如这里的a = X \ y，可以直接看到命令执行结果）
a =
    0.7190
   32.8645
    7.1925
    2.6268
这时候，我们自然要问，如何得到这些参数的误差呢？

要回答这个问题，我们要构建一个误差矩阵，但是在构建误差矩阵之前，先要计算我们拟合的模型的误差：

y_predict = X*a;         .................... （通过这一步，我们得到拟合的y值）

然后根据定义：
（模型的误差）^2 = sum（（实际y值 - 拟合的y值）^2）/ （自变量总数 - 参数总数） 求得：

s = sqrt( sum((y - y_predict).^2)  / (100 -4) );

于是，a 的误差可以这样表示：

XTX_inverse = (X'*X)^-1;

da = s * sqrt( diag( XTX_inverse , 0) )

结果如下：
da =
    0.0188
    0.1561
    0.3645
    0.2395

意思是：
A1 =  0.7190   +/- 0.0188
A2 =  32.8645 +/- 0.1561
A3 =  7.1925   +/- 0.3645
A3 =  2.6268   +/- 0.2395

这里的误差是标准误差，涵盖了大约68%的置信度，假如需要更大的置信区间，可以查表得到t值，
比如 N= 100 - 4 的 95% 的置信区间 t=1.98，然后用
误差 = t * da得到。如：

da_95= 1.98 * da

da_95 =
    0.0372
    0.3090
    0.7216
    0.4742

最后，补充几点：
1） 这里的 100 - 4 表示自由度，我们的自变量有100个，参数有4个，
自由度 = 自变量总数 - 参数总数

2）关于 XTX_inverse = (X'*X)^-1;
注意，这里不是 “.^” ，Matrix ^-1  这个运算在Matlab矩阵里面也是一个求矩阵逆的运算。

还有什么问题，欢迎回贴提问。

谢谢！