本帖最后由 xiaowu55 于 2016-3-25 22:30 编辑
 
flag_go_on=1;
 
num_of_loop=0;
 
a = [];%存储剔除的变量下标
 
xx = 1:k;%存储最原始的变量下标,如果有5个变量x,则存储1,2,3,4,5
 
while flag_go_on
 
cij=inv(X'*X);
 
cii=diag(cij);
 
F_fenweidian_1=finv(1-F_alpha,1,n-k-1);
 
ci=sqrt(cii(2:end)*Se_square*F_fenweidian_1/(n-k-1));
 
format_str='%15.4f';
 
for ii=1:k-1
 
format_str=[format_str '%13.4f'];
 
end
 
fprintf(['\r第%d次检验：\rcii: ' format_str '%13.4f\r ci:              ' ...
 
format_str '\rβi：' format_str '%13.4f'],num_of_loop+1,cii,ci,beta_mao)
 
if ~all(abs(beta_mao(2:end))>ci')
 
flag_go_on=1;
 
beta_1tok=beta_mao;
 
beta_1tok(1)=[];
 
fi_xin=beta_1tok.^2./cii(1:end-1)';
 
min_fi=min(fi_xin);
 
beta_index=find(fi_xin==min_fi)+1;
 
%     这样就可以输出正确的下标了
 
index_xx = xx(beta_index-1);
 
fprintf('\r x%d 对y的线性影响最不显著(|β%d|=%0.4f)。 删除x%d， 进行第%d次计算：',...
 
index_xx,index_xx,abs(beta_all(beta_index)), index_xx,loop_num + 2);
 
xx(beta_index-1)=[];
 
%         fprintf('\rx%d对y的线性影响最不显著( |β%d|=%0.4f )。删除x%d，进行第%d次计算：',...
 
%         beta_index-1+num_of_loop,beta_index-1+num_of_loop,...
 
%         abs(beta_mao(beta_index)),beta_index-1+num_of_loop,...
 
%         num_of_loop+2)
 
else
 
fmt_str2='x%d';
 
index_of_xi=find(index_of_xi_array);
 
for i2=1:length(find(index_of_xi))-1
 
fmt_str2=[fmt_str2 '、x%d'];
 
end
 
fprintf(['\r\r经过检验，剩余所有变量:' ...
 
fmt_str2 '对y的线性影响均显著。检验结束。\r'],index_of_xi)
 
flag_go_on=0;
 
end
 
if flag_go_on
 
num_of_loop=num_of_loop+1;
 
k=k-1;
 
if ~k
 
fprintf('\r\r警告：通过一一对所有变量做显著性检验，已剔除所有变量！');
 
break;
 
end
 
beta_mao=beta_mao-beta_mao(beta_index)/cii(beta_index)*cij(beta_index,:);
 
beta_mao(beta_index)=[];
 
%    这样更改就没问题了
 
k1 = (beta_index-1)*10+1;
 
k2 = (beta_index)*10;
 
k = k1:k2;
 
format_str0(k)=[];
 
fprintf(format_str0,beta_mao);
 
%     fmt_str1='';
 
%     for i1=1:beta_index-2
 
%         fmt_str1=[fmt_str1 'β' num2str(i1) ' = %0.4f\r'];
 
%     end
 
%    for i1=beta_index:k+1
 
%         fmt_str1=[fmt_str1 'β' num2str(i1-1+num_of_loop) ' = %0.4f\r'];
 
%     end
 
%     fprintf(['\rβ0 = %0.4f\r' fmt_str1],beta_mao)
 
X(:,beta_index)=[];
 
index_of_xi_array(beta_index-1+num_of_loop-1)=0;
 
xi=X(:,2:end);
 
x_ba=mean(xi);
 
lxy=sum((xi-ones(n,1)*x_ba).*((Y-y_ba)*ones(1,k)));
 
Sr_square=sum(beta_mao(2:end).*lxy);
 
Se_square=St_square-Sr_square;
 
end
 
end
 
这样更改后显示显著性的系数就没有问题了。

来源于机器学习实战中p138，求解线性回归的回归系数w的最优解，涉及到矩阵求导等知识，推导过程中还对矩阵求导的分子、分母布局进行说明，部分参考链接如下： 
 1.https://en.wikipedia.org/wiki/Matrix_calculus#Denominator-layout_notation（wiki矩阵求导说明） 
 2.http://f.dataguru.cn/thread-598486-1-1.html（最小二乘估计、岭回归和lasso） 
 3.http://blog.csdn.net/shouhuxianjian/article/details/46669365（对wiki矩阵求导进行了翻译） 
 **
 
背景
 
** 
 问题的背景是从大量数据中求出回归方程，常用的方法就是找出误差最小的w（存放回归系数的向量）。一般采用的方法是下图即平方误差（这里X是n*p的矩阵，y是n*1的结果向量，w便是p*1的回归系数向量），求和公式是误差项的平方值，我们需要做的就是取求和公式最小值下的w
 
 
 
 
 
用矩阵表示还可以写作
 
 
 
 
 
下面也就是对上述公式进行求导，令其导数等于0，所得w即为最优回归参数
 
**
 
求导过程
 
** 
 敲公式太麻烦就放了张手写图
 
 
1-4步涉及到了矩阵求导知识和分子布局、分母布局两种不同的规定格式（具体见wiki）
 
1.先将分子展开成四项，然后就可以对四个部分分别进行求导（累加后求导等于求导后累加）
2.步骤1中的分子为标量（1*n向量与n*1向量），属于表格中第一种情形，即分子为标量分母为向量，因此求导结果为0
3.同样，步骤2，3中对应表格第3行，步骤4对应表格的第11行。但是涉及到不同的布局最后的表现形式也不同，如本次采用的是分母布局（也就是行向量和列向量的表示区别）
 
**
 
布局约定
 
**
 
这里引用第三个连接中对wiki中矩阵求导的翻译
 

一、多元线性回归
 
所谓的多元线性回归就是指在回归分析中，如果有两个或两个以上的自变量，就称为多元回归。
 
二、多元线性回归模型
 
1.建立模型
 以二元线性回归模型为例 ，二元线性回归模型如下：
 
 类似的使用最小二乘法进行参数估计 ：
 
 2.拟合优度指标
 标准误差：对y值与模型估计值之间的离差的一种度量。其计算公式为：
 
 3.置信范围
 置信区间的公式为：置信区间=
 其中， 是自由度为 的 统计量数值表中的数值, 是观察值的个数， 是包括因变量在内的变量的个数。
 
三、估值方法
 
1.普通最小二乘法
 普通最小二乘法(Ordinary Least Square, OLS)通过最小化误差的平方和寻找最佳函数。通过矩阵运算求解系数矩阵：
 
 2.广义最小二乘法
 广义最小二乘法(Generalized Least Square)是普通最小二乘法的拓展，它允许在误差项存在异方差或自相关，或二者皆有时获得有效的系数估计值。
 
 
四、推导过程
 

 
 
五、关于矩阵的计算的程序
 
我之前写了一个矩阵计算的轮子，详情请参考这篇文章：
 《python实现矩阵操作（自造轮子）》

https://editor.csdn.net/md/?articleId=105137945
 其实上一篇讨论的多元线性回归还不是很全面，而且训练和测试集都没有分，在这里继续讨论多元线性回归模型检验的相关问题。
 
只要有P值的出现，样本量不超过5000，比如线性回归和逻辑回归；搞清楚算法背后的逻辑才是比较重要的。
 
多元线性回归需要关注一些什么点？R2和模型稳定性，也就是那些β是不是稳定的，检验模型是不是稳定需要对模型进行诊断。
 
多元线性回归的输出变量是连续变量，输入变量中如果包含离散变量需要做哑变量或One-hot编码，连续变量可以直接用来建模。
 
 
多元线性回归假设解释
 
多元线性回归需要满足的假设其实是比较强的，但是在机器学习或者是数据挖掘领域，后3条针对误差项（其实就是残差）的假设基本上被忽略了。
 第1条： 看因变量y和自变量x之间的关系，可以通过绘制散点图，确定是线性、二次函数还是指数函数关系，根据这个来建立x和y之间的关系。后面的神经网络和SVM等模型就是为了方便找到x和y之间的关系。
 第2条： 解释变量和随机扰动项不存在线性关系。我们想象一下，如果他们之间存在线性关系的话，是不是会导致回归系数估计的不准确啊，举个例子解释变量y是收入，x是受教育程度，并假设回归方程是 y = 0.5x + e，设想如果扰动项里面包含父母收入，实际上父母收入会影响孩子的收入y，那么回归系数估计值0.5是不是偏高了啊。那怎么解决这个问题呢，那就多纳入一些变量来参与建模吧，这也是多元线性回归存在的必要性，同时这也引出了一个变量筛选的问题。
 第3条：解释变量之间不存在线性关系（或强相关）。在建模时，我们不但需要估计回归系数的均值，还需要估计回归系数的标准差：S（β）= S（e）/ |x|，那么如果解释变量x之间存在线性关系的话，分母趋向于0了，回归系数标准差趋于无穷大，所以多重共线性问题是需要去避免的。
 
多元线性回归诊断方法
 
如果扰动项是右偏，那么残差图肯定是异方差分布，取对数即可，所以下图中假设5和假设4可以说是一致的。
 
 
多元线性回归模型的诊断
 
（1）残差分析：实际上当残差不包含任何信息的时候是最好的，如果还包含一些信息，需要把这个信息提取出来。残差图的纵坐标是残差，横坐标可以是各个解释变量x，实际上在做单变量分析，解释变量x被解释变量y做相关性分析的时候就知道了；比如某个解释变量x和被解释变量y都是右偏，那么残差图肯定是异方差，同时取对数重新建模；如果某个解释变量x和被解释变量y存在抛物线关系，那么加入二次项重新进行建模；自相关一般在时间序列数据中比较常见。
 
 
 （2）强影响点分析：
 为什么要做强影响点分析？？？因为有了强影响点的存在之后，会把本来没有关系的数据带出关系来，而且这个关系特别不稳定。比如下图，本来数据点之间没有什么关系，但是因为强影响点的存在之后，给数据带出来了这么一个线性关系出来，但这个关系是非常不稳定的，不具有大众性。
 
 怎么解决强影响点分析问题？？？学生化残差（只做一次）。
 
|SR| = (残差 - 残差均值) / 残差标准差。
 |SR| > 2，剔除掉满足条件的记录（几百个样本）
 |SR| > 3，剔除掉满足条件的记录（几千个样本）
 
（3）共线性问题
 可以参考下面的链接：https://www.sohu.com/a/326904117_100103806
 共线性的判别指标：膨胀系数VIF、相关系数
 共线性的解决方法：根据VIF和相关系数手动剔除变量、逐步回归法、岭回归。