layout: post
title: matlab稳健回归函数文档
categories:
- 科
tags:
- matlab
math: true
regress函数和regstats函数利用普通最小二乘法估计模型中的参数, 参数的估计值受异常值的影响比较大. robustfit函数采用加权最小二乘法估计模型中的参数, 受异常值的影响就比较小. robustfit函数可用来作稳健的多重线性或广义线性回归分析, 下面介绍robustfit函数的用法.
调用方法
<div class="highlight"><pre style="line-height:125%"><span></span>b = robustfit(X,y)
b = robustfit(X,y,wfun,tune)
b = robustfit(X,y,wfun,tune,const)
[b,stats] = robustfit(<span style="color: #008800; font-style: italic">...)</span>
</pre></div>使用说明
1. b = robustfit(X, y)
通过执行稳健回归来分析多元线性回归模型 $y= X \b$, 并返回系数向量 $\b$ 的估计值 $b$.
输入参数 $X$ 为 $n \times p$ 的自变量矩阵(或称预测变量矩阵, 设计矩阵), 对应 $p$ 个预测因子对 $n$ 个观测值中每个的贡献. $y$ 是 $n \times 1$ 观测值向量(或称响应向量), 输出的 $b$ 为 $(p + 1) \times 1$ 向量.
缺省情况下, 算法使用基于bisquare加权函数的迭代重加权最小二乘法.
注意 与regress函数不同的是, 默认情况下, robustfit函数会自动在 $X$ 第1列元素的左边加入一列1, 而不需要用户自己添加. 此列向量对应于模型中的常量项. 不要直接为 $X$ 添加一个全1的列向量, 你可以通过更改变量const的值来改变robustfit的默认行为.
robustfit会把 $X$ 或 $y$ 中的不确定数据NaN作为缺失数据, 并将其移除.
2. b = robustfit(X, y, wfun, tune)
指定加权方法wfun和调节常数tune. 在计算权重之前tune会被划分到残差向量. 如果指定了wfun, 那么tune必不可少.
加权方法wfun为字符串, 可以取下表中的任何一个.
<table id='tab-0'><caption>robustfit支持的加权方法</caption>
<tr>
<th rowspan="1" colspan="1" style="text-align:center;">加权方法</th>
<th rowspan="1" colspan="1" style="text-align:center;">权重函数</th>
<th rowspan="1" colspan="1" style="text-align:center;">默认调节常数</th>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'andrews'</td>
<td rowspan="1" colspan="1" style="text-align:center;">$w = \sin(r)/r \;\;\text{if}\; \abs r <\p$</td>
<td rowspan="1" colspan="1" style="text-align:center;">1.339</td>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'bisquare'(默认)</td>
<td rowspan="1" colspan="1" style="text-align:center;">$w = (1-r^2)^2 \;\;\text{if} \abs r<1$</td>
<td rowspan="1" colspan="1" style="text-align:center;">4.685</td>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'cauchy'</td>
<td rowspan="1" colspan="1" style="text-align:center;">$w = 1/(1 + r^2)$</td>
<td rowspan="1" colspan="1" style="text-align:center;">2.385</td>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'fair'</td>
<td rowspan="1" colspan="1" style="text-align:center;">$w = 1/(1 + \abs r)$</td>
<td rowspan="1" colspan="1" style="text-align:center;">1.400</td>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'huber'</td>
<td rowspan="1" colspan="1" style="text-align:center;">$w = 1/\max(1, \abs r)$</td>
<td rowspan="1" colspan="1" style="text-align:center;">1.345</td>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'logistic'</td>
<td rowspan="1" colspan="1" style="text-align:center;">$w = \tanh(r)/r$</td>
<td rowspan="1" colspan="1" style="text-align:center;">1.205</td>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'ols'</td>
<td rowspan="1" colspan="1" style="text-align:center;">普通最小二乘法(无加权)</td>
<td rowspan="1" colspan="1" style="text-align:center;">无</td>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'talwar'</td>
<td rowspan="1" colspan="1" style="text-align:center;">$w = 1 \;\; \text{if}\; \abs r<1$</td>
<td rowspan="1" colspan="1" style="text-align:center;">2.795</td>
</tr>
<tr>
<td rowspan="1" colspan="1" style="text-align:center;">'welsch'</td>
<td rowspan="1" colspan="1" style="text-align:center;">$w = \exp(-r^2)$</td>
<td rowspan="1" colspan="1" style="text-align:center;">2.985</td>
</tr>
</table>若调用时没有指定调节常数tune, 则使用表中的默认值. 默认调节常数给出的系数估计值约为普通最小二乘估计的95%, 前提是响应服从正态分布且无异常值. 减小调节常数会增加分配给大残差的downweight; 增加调节常数会减少分配给大残差的downweight.
权重函数中的 $r$ 值通过下式计算
$$r = {\text{resid} \over \text{tune}s\sqrt{1-h}}$$
其中 $\text{resid}$ 是上一次迭代中残差的向量, $h$ 是由最小二乘拟合得到的中心化杠杆值向量, $s$ 是误差项标准偏差的估计值, 计算公式为
$s = \text{MAD}/0.6745$
其中 $\text{MAD}$ 为残差绝对值的中位数. 常数0.6745保证了在正态分布下估计是无偏的. 如果 $X$ 中有 $p$ 列, 则在计算 $\text{MAD}$ 时会将残差绝对值向量的前 $p$ 个最小值舍去.
用户可以自己定义权重函数. 该函数的输入必须是残差向量, 并输出权重向量. 在这种情况下, 调用robustfit函数时把自定义权重函数的句柄(形如@myfun)作为wfun参数传递给robustfit函数, 此时必须指定tune参数
3. b = robustfit(X, y, wfun, tune, const)
用参数const来控制模型中是否包含常数项. 若const取值为'on'(默认值), 则模型中包含常数项, 此时会自动在 $X$ 第1列的左边加入一列1, $b$ 变为 $(p + 1) \times 1$ 向量. 若const取值为'off', 模型中不包含常数项, 此时不改变X的值. 则 $b$ 为 $p \times 1$ 向量.
4. [b, stats] = robustfit(...)
返回一个结构体变量stats, 其字段包含了用于模型诊断的统计量. stats的字段为:
ols_s: 普通最小二乘法得出的σ的估计(RMSE)robust_s: σ的稳健估计mad_s: 用残差绝对值的中位数计算σ的估计; 用于在迭代拟合过程中缩放残差s: σ的最终估计, 是ols_s和robust_s加权平均值与robust_s值中的较大者resid: 残差rstud: 残差的t检验值(详细信息参见regress)se: 系数估计的标准误差covb: 系数向量协方差矩阵的估计coeffcorr: 系数向量相关系数矩阵的估计t: b与se的比值p t检验的p值w: 稳健拟合的权重向量R: 矩阵X的QR分解中的R因子dfe: 误差的自由度h: 最小二乘拟合的中心化杠杆值向量
robustfit函数使用 inv(X'*X)*stats.s^2来计算系数估计值的方差协方差矩阵. 标准误差和相关性由此得出.
参考译文