最小二乘法是对参数估计的一种算法，其优点是不需要系统的先验知识，所以是一种系统辨识很重要的算法。下面对其进行简单的介绍：
 
一般最小二乘法
这里假设状态x为n*1维的列向量，而测量值为状态x的线性组合，假设进行了L次观测，则局部观测模型有：
 
 
利用所有的观测信息合成观测模型则有：
 
 
其中有：
 
 
根据上述模型，给出最小二乘的极小化性能指标：
 
 
对最小化性能指标J求偏导，可得：
 
 
则有一般最小二乘估计，这样估计的误差均方差阵为：
 
 
加权最小二乘法
在没有测量值误差的统计信息时，采用一般最小二乘法，但是一般最小二乘法没有考虑两种测量仪器测量精度不同的特点，如果在已知测量值误差统计信息时，我们就可以采用加权最小二乘法，进一步提高状态估计精度。
 
加权最小二乘法与一般最小二乘法不同的地方是最小化性能指标发生改变，加入了权重矩阵，这里我们选择让测量精度高的数据权重值大一些，而测量精度低的数据权重值小些，构成如下最小化性能指标：
 
 
其中权重矩阵为：
 
 
对最小化指标求偏导可得：
 
 
同理求误差方差阵如下：
 
 
可以证明估计后状态误差小于每个传感器的误差。

解决异方差问题--加权最小二乘法

参考文章：
 
（1）解决异方差问题--加权最小二乘法
 
（2）https://www.cnblogs.com/jin-liang/p/9537539.html
 

备忘一下。

 [回归分析][9]--加权最小二乘法
 

 
 
 

   这一节会讲一下关于加权最小二乘法，这种方法是用来处理 ”异方差“ 的。如下图
 

   
 

 
 
 

 关于处理的方法：
 

 
 

 
 
 

 
 
 

 
 
 

 例子：我们来看一个直观的例子
 

 
 
 

 1.构造一组数据
 

 
 
x = Table[i, {i, 1, 1000}];
y = Table[2*i + RandomInteger[PoissonDistribution[Ceiling[i/100.0]]], {i, 1,1000}];
 
 
 

 2.对数据拟合
 
 

 
data = Transpose[{x, y}];
Clear[x];
lm = LinearModelFit[data, x, x]
 

 
 
 

 3.画出残差
 

 
cancha = lm["FitResiduals"];
p1 = ListPlot[cancha]
 

 
 
 
 

 
 
 

 4.按类别，看一下彩色的图
 

 
temp = {};
For[i = 1, i <= 10, i++,
  temp =  AppendTo[temp, cancha[[100*i - 99 ;; 100*i]]]
  ];
dataList = 
  Transpose[{Range @@ #[[1]], #[[2]]}] & /@ 
   Transpose@{# + {1, 0} & /@ 
      Partition[Accumulate[Length /@ Prepend[temp, {}]], 2, 1], temp};
p4 = ListPlot[dataList]
 

 
 

 可以看到残差随着x的增大而增大。
 
 
 
 
 

 5.求出加权的系数
 

 
 
 

 计算总的残差
 
 

 
canchaTotal = Total[Abs[cancha]]/Length[cancha]
 

 
 
 

 计算每一组的残差
 

 
canchazu = Table[0, {i, 1, 10}];
For[i = 1, i <= 10, i++,
 canchazu[[i]] =  Total[Abs[cancha[[100*i - 99 ;; 100*i]]]]/99
 ]
 

 
 
 

 计算系数
 

 
xishu = canchazu/canchaTotal
 
 
 

 计算权重
 

 
weight = Table[xishu, {i, 1, 100}];
weight = Flatten@Transpose[weight];
 

 
 6.对加权后的重新拟合
 

 
lmn = LinearModelFit[data/weight, x, x]
 

 查看其残差
 

 
temp = {};
For[i = 1, i <= 10, i++,
 temp =  AppendTo[temp, canchan[[100*i - 99 ;; 100*i]]]
 ]
dataList = 
  Transpose[{Range @@ #[[1]], #[[2]]}] & /@ 
   Transpose@{# + {1, 0} & /@ 
      Partition[Accumulate[Length /@ Prepend[temp, {}]], 2, 1], temp};
p3 = ListPlot[dataList, ImageSize -> Medium]
 

 
 
 可以看到残差比之前要好一点
 

 
 
 

 可以放在一起比较一下
 

 
Row[{p4, p3}]
 

 
 
 有较为明显的区别。
 

 
 
 

 以上，所有
 

 2016/11/20
 
 
 
 
 

 
 
 

 
 
 

 
 
 

 
 

 
 
 
  
 
   
 一.前言
 
   
 在往前可以看得见的历史里，我们漫长的一生中不知道做了多少个回归，然而并不是每一个回归都尽如人意，究其原因就很多了，可能是回归的方程选择的不好，也可能是参数估计的方法不合适。回归的本质是在探寻因变量Y和自变量X之间的影响关系(PS：国外的论文里常常叫因变量为响应变量，自变量为解释变量)，如何来描述这种相关关系呢？我们可以假设Y的值由两部分组成，一部分是X能决定的，记为
     
      f(X)
     ,另一部分由其它众多未加考虑的因素（如随机因素）所组成，记为随机误差
     
      e
     ,并且我们有理由认为
     
      E(e)=0
     .于是我们得到
 
   
 
    
     
      Y=f(X)+e
     
   
 特别地，当
   
   
    
     f(X)
    是线性函数时，我们便得到了众多回归组成的王国里最平民也最重要的回归模型——线性回归。回归参数的估计方法最基本的就是最小二乘回归。尽管长江后浪推前浪，我们有了更多回归参数的估计方法（详见优矿量化实验室），但依然不影响我们对它的喜爱，为了介绍后续的回归方法，我依然要先把普通最小二乘法的回归模型摆出来致敬一下:
    
   
 
    
     
      Y=Xβ+e
     
   
 这里的Y就是$n
   1
     
      的变量观测向量，X为
     np
   
   
    
     的已知设计矩阵，
    \beta
   
   
    
     为
    p*1
   
   
    
     未知参数向量，
    e$为随机误差向量。（这里
   需要注意的是由于有常数项，所以自变量个数其实是p-1个） 
   
 
   
 
 
   
 
    
     
      (y1 y2 ⋮ yn)=(1x11⋯x1,p−1 1x21⋯x2,p−1 ⋮⋮⋮⋮ 1xn1⋯xn,p−1)(β0 β1 ⋮ βp−1)+(e1 e2 ⋮ en)
     
   
 
   
 
   
 Gauss-Markov假设可以简写为
     
      E(e)=0,Cov(e)=σ2In，Cov(e,x)=0
     
 
   
 普通最小二乘法就是使得残差平方和
 
   
 
    
     
      RSS(β)=||Y−Xβ||2=(Y−Xβ)T(Y−Xβ)
     
   
 最小,通过对矩阵的求导，我们就得到了
   
   
    
     β
    的估计
   
   
    
     (XTX)−1XTY
    ，在假设下，我们可以证明该估计是
   
   
    
     β
    所有线性无偏估计中方差最小的。 
   
 
   
 与之相关的各种分布和检验就不再赘述，社区里也有很多帖子可以学习。既然模型有假设，那么局限性就出来了，我们常常会发现残差项并不满足假设，尤其当变量是时间序列时，非平稳性和自相关性常常会造成异方差的问题，那怎么处理异方差的问题呢？方法也挺多的，本帖主要讨论一种加权最小二乘估计的方法和机器学习里局部加权线性回归的方法。
 
   
 (注：本帖的线性回归是普适的多元线性回归，所以均用矩阵来表示，更方便)
 
  
 
 
 
 
 
 
 
  
 
   
 二.加权最小二乘法
 
   
 加权最小二乘法其实是广义最小二乘法的一种特殊情形，而普通最小二乘法也是一种特殊的加权最小二乘法。为了保证知识的完整性，不妨把广义最小二乘介绍一下：
 
   
 1.广义最小二乘法
 
   
 刚才说到我们的残差项项不满足Gauss-Markov假设，那么我们就把假设放宽一些：考虑以下模型：
 
   
 
    
     
      Y=Xβ+e,E(e)=0,Var(e)=σ2Σ,Cov(e,x)=0
     
   
 这里的
   
   
    
     Σ
    是我们已知的一个
   
   
    
     n∗n
    正定对称矩阵，其中
   
   
    
     σ2
    不一定是已知的。也就是说不要求误差项互不相关了。 
   
 
   
 这里我们的广义最小二乘法就是使得广义残差平方和
 
   
 
    
     
      RSS(β)=(Y−Xβ)TΣ−1(Y−Xβ)
     
   
 最小，最后
   
   
    
     β
    的估计为
   
   
    
     (XTΣ−1X)−1XTΣ−1Y
     
   
 
   
 实际上，
     
      Σ
     也常常是未知的，但当我们知道
     
      Σ
     的某种形式时，我们可以去估计它。举个例子，如果出于某种原因，我们样本的数据来源的地方不一样，怎样把他们整合在一起呢，我们可以假设那些来源相同的数据样本的残差项方差是一样的，如：
 
   
 
    
     
      Cov(e)=(σ21In1⋯00 0σ22In2⋯0 ⋮⋮⋮⋮ 00⋯σ2kInk)
     
   
 然后我们再通过迭代的方法（第一步就是普通的最小二乘回归）去估计
   
   
    
     σ2i
    ，直到相邻两次迭代求得的
   
   
    
     β
    的估计差不多为止。 
   
 
   
 从上面可以看出，不论你对自变量和因变量作了何等变换，最终都可以用最小二乘的模型去解决回归的问题。但普通最小二乘估计里的那些美好的性质，分布和检验等在广义最小二乘里还存在么？还一致么？其实，我们可以通过一些变换，把广义最小二乘的模型转换为满足普通最小二乘法假设的模型，关于这一点，我们以下面的加权最小二乘法举例说明。
 
  
 
 
 
 
 
 
 
  
 
   
 2.加权最小二乘法
 
   
 加权最小二乘法，就是对上述的
     
      Σ
     取一种特殊的矩阵--对角阵，而这个对角阵的对角元都是常数，也就是权重的倒数，如下：
 
   
 
    
     
      Cov(e)=σ2Σ=σ2(1/ω1⋯00 01/ω2⋯0 ⋮⋮⋮⋮ 00⋯1/ωn)
     
   
 
   
   
    
     ωi
    表示的就是第i个样本在回归里的权重，从上式可以看出来，具有较大权的样本具有较小的方差，它在回归问题里显得更加重要。不妨用W来表示权重矩阵，那么
   
   
    
     Σ−1=W
    ,这时我们用广义最小二乘的方法来求系数的估计，即最小化广义残差平方和
    
   
 
    
     
      RSS(β)=(Y−Xβ)TW(Y−Xβ)
     
   
 
   
   
    
     β
    最后的估计结果为：
   
   
    
     β^=(XTWX)−1XTWY
     
   
 
   
 我们并不满足于求出广义的最小二乘估计，我们还要研究它的很多性质，下面上面提过的广义最小二乘转换的方法把这个模型转化为满足普通最小二乘回归假设的模型：
 
   
 首先我们找一下
     
      Σ−1
     的平方根C，在这个问题里，很容易得到C：
 
   
 
    
     
      C=(ω−−√1⋯00 0ω−−√2⋯0 ⋮⋮⋮⋮ 00⋯ω−−√n)
     
   
 对回归模型每一项乘以C，得到：
    
   
 
    
     
      CY=CXβ+Ce
     
   
 
   
 
   
 这时候
     
      Ce
     的协方差阵为：
 
   
 
    
     
      Var(Ce)=σ2CC−1(C−1)TCT=σ2In
     
   
 这就是满足Gauss-Markov假设的普通线性回归模型了 
   
 
   
 不妨重新对变量命名：
 
   
 
    
     
      Z=Mβ+d
     
   
 其中，
    
   
 
    
     
      Z=CY=(ω−−√1y1 ω−−√2y2 ⋮ ω−−√nyn)
     
   
 
    
   
 
    
     
      M=CX=(ω−−√1ω−−√1x11⋯ω−−√1x1,p−1 ω−−√2ω−−√2x21⋯ω−−√2x2,p−1 ⋮⋮⋮⋮ ω−−√nω−−√nxn1⋯ω−−√nxn,p−1)
     
   
 
    
   
 
    
     
      d=Ce
     
   
 感兴趣的读者可以验证一下，这个新模型用普通最小二乘所估计出来的
   
   
    
     β
    和原模型是一样的，而且线性无偏方差最小的性质和分布，检验等都可以用起来了，如R2及显著性检验等来看拟合的好坏。 
   
 
  
 
 
 
 
 
 
 
  
 
   
 三.局部加权线性回归
 
   
 局部加权线性回归是机器学习里的一种经典的方法，弥补了普通线性回归模型欠拟合或者过拟合的问题。机器学习里分为无监督学习和有监督学习，线性回归里是属于有监督的学习。普通的线性回归属于参数学习算法(parametric learning algorithm)；而局部加权线性回归属于非参数学习算法(non-parametric learning algorithm)。所谓参数学习算法它有固定的明确的参数，参数 一旦确定，就不会改变了，我们不需要在保留训练集中的训练样本。而非参数学习算法，每进行一次预测，就需要重新学习一组 ， 是变化的，所以需要一直保留训练样本。也就是说，当训练集的容量较大时，非参数学习算法需要占用更多的存储空间。
 
   
 上面的话好像不够直观，下面我们来直观地看一下局部加权线性回归到底是怎么样的。
 
   
 局部加权线性回归就是在给待测点附近的每个点赋予一定的权重，也就是加一个核函数矩阵W，最后需要最小化的目标函数如下：
 
   
 
    
     
      ∑iwi(yi−y^i)2
     
   
 这是什么，，，这不就刚刚的加权最小二乘法么？再来看看系数的估计项：
    
   
 
    
     
      β^=(XTWX)−1XTWY
     
   
 
    
   
 
    
     
      W=(ω1⋯00 0ω2⋯0 ⋮⋮⋮⋮ 00⋯ωn)
     
   
 好吧，这就是加权最小二乘法，当然，如果这就是局部加权线性回归的全部，也不用发明它了，我们接着看。 
   
 
   
 这里出现了一个问题，权重是如何确定的，也就是这个所谓的核函数矩阵的形式，百度了一下，通常使用的都是高斯核,形式如下：
 
   
 
    
     
      ω(i,i)=exp(−(xi−x)2/2k2)
     
   
 
   
   
    
     xi
    代表的是第i个样本点，
   
   
    
     x
    是我们预测点，对于金融数据来说，完全可以用时间t来衡量，也就是说时间越近的样本数据越重要，这样的预测对想研究的对象而言更准确。 
   
 
   
 介绍到这里似乎没有体现出机器学习的意思，仔细观察就会发现Ｋ是一个很重要的东西，不信，我们举个例子来看看：
 
   
 例子：
 
  
 
 
 
 
 
 
 
  
 
   
 
    
 
     
 
      
     
     

      

       

        

         

          

           

           

            

           
          
          

          

          

           

          
          

           

            

             

              1
             
            
            
import numpy as np
           
           

            

             

              2
             
            
            
import pandas as pd
           
           

            

             

              3
             
            
            
import scipy.stats as stats
           
           

            

             

              4
             
            
            
from math import *
           
          
         
        
       
      
      

      

       

       

      
     
    
   
  
 
 查看全部


 

  

   
为了简单起见，我们就用一元的线性回归来看看，方便数据可视化。这里我们取一些非线性拟合的数据，因为局部加权线性回归的优势就在于处理非线性关系的异方差问题。
  
 


 

  

   

    

     

      
     
 
     
 
      
 
       
 
        
 
         
 
          
 
           
 
           
 
            
 
           
 
          
 
          
 
          
 
          
 
          
 
           
 
           
 
          
 
          
 
           
 
            
 
             

               1
             
 
            
 
            
#x取1~100,y取某一个非线性函数
 
           
 
           
 
            
 
             

               2
             
 
            
 
            
x = np.arange(1,101)
 
           
 
           
 
            
 
             

               3
             
 
            
 
            
x=np.array([float(i) for i in x])
 
           
 
           
 
            
 
             

               4
             
 
            
 
            
y =x+[10*sin(0.3*i)for i in x]+stats.norm.rvs(size=100, loc=0, scale=1.5)
 
           
 
          
 
         
 
        
 
       
 
      
 
      
 
      
 
      
 
       
 
       
 
       
 
       
 
      
 
     
 
    
 
   
 
  
 
 
 
 查看全部
 
 
 
 
  
 
   
 
    
 
     
 
      
     
     

      

       

        

         

          

           

           

            

           
          
          

          

          

           

          
          

           

            

             

              1
             
            
            
#先加载画图需要用的包
           
           

            

             

              2
             
            
            
import matplotlib as mpl
           
           

            

             

              3
             
            
            
import matplotlib.pyplot as plt
           
           

            

             

              4
             
            
            
mpl.style.use('ggplot')
           
           

            

             

              5
             
            
            
import seaborn as sns
           
          
         
        
       
      
      

      

       

       

      
     
    
   
  
 
 查看全部


 

  

   
看一下我们需要拟合的样本数据：
  
 


 

  

   

    

     

      
     
 
     
 
      
 
       
 
        
 
         
 
          
 
           
 
           
 
            
 
           
 
          
 
          
 
          
 
          
 
          
 
           
 
           
 
          
 
          
 
           
 
            
 
             

               1
             
 
            
 
            
plt.figure(figsize=(12,6))
 
           
 
           
 
            
 
             

               2
             
 
            
 
            
plt.scatter(x,y)
 
           
 
          
 
         
 
        
 
       
 
      
 
      
 
      
 
      
 
       
 
       
 
       
 
       
 
      
 
     
 
    
 
   
 
  
 
 
 
 查看全部 
 
 
  
 
   

    <matplotlib.collections.PathCollection at 0x8bc4d10>
   
 
   

    
   
 
  
 
 
 
 
 
 
 
  
 
   
 很明显，这是一个非线性关系的样本数据，我们先用普通最小二乘回归来处理这个问题：
 
  
 
 
 
 
 
 
 
  
 
   
 
    
 
     
 
      
     
     

      

       

        

         

          

           

           

            

           
          
          

          

          

           

          
          

           

            

             

              1
             
            
            
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
           
          
         
        
       
      
      

      

       

       

      
     
    
   
  
 
 查看全部


 

  

   

    

     

      
     
 
     
 
      
 
       
 
        
 
         
 
          
 
           
 
           
 
            
 
           
 
          
 
          
 
          
 
          
 
          
 
           
 
           
 
          
 
          
 
           
 
            
 
             

               1
             
 
            
 
            
plt.figure(figsize=(12,6))
 
           
 
           
 
            
 
             

               2
             
 
            
 
            
yHatLinear=intercept+slope*x
 
           
 
           
 
            
 
             

               3
             
 
            
 
            
plt.plot(x,yHatLinear,'r')
 
           
 
           
 
            
 
             

               4
             
 
            
 
            
plt.scatter(x,y)
 
           
 
           
 
            
 
             

               5
             
 
            
 
            
print 'y='+str(intercept)+'+'+str(slope)+'x'
 
           
 
          
 
         
 
        
 
       
 
      
 
      
 
      
 
      
 
       
 
       
 
       
 
       
 
      
 
     
 
    
 
   
 
  
 
 
 
 查看全部 
 
 
  
 
   

    y=1.34568655273+0.979532165941x
   
 
   
 
   

    
   
 
  
 
 
 
 
 
 
 
  
 
   
 可以看到，要用直线来拟合非线性关系略有牵强，这个例子还算举的不错，金融数据里很多时间序列的关系都是非线性的，回归的结果往往不好．
 
   
 下面我们用刚才介绍的局部加权线性回归来拟合一下这个模型，简单回顾一下过程：
 
   
 1.用高斯核函数计算出第ｉ个样本处，其它所有样本点的权重Ｗ
 
   
 2.用权重ｗ对第ｉ个样本作加权线性回归，得到回归方程，即拟合的直线方程
 
   
 3.用刚才得到的经验回归直线计算出
     
      xi
     处的估计值
     
      y^i
     
 
   
 4.重复一至三步，得到每个样本点的估计值
 
   
 这里作加权线性回归时，我使用的是把加权最小二乘转换为普通最小二乘的方法，也就是本帖第二部分内容，网上的做法大多是直接用上面的公式算出
     
      β
     的估计值。
 
  
 
 
 
 
 
 
 
  
 
   
 
    
 
     
 
      
     
     

      

       

        

         

          

           

           

            

           
          
          

          

          

           

          
          

           

            

             

              1
             
            
            
#局部加权线性回归里的权重的平方根
           
           

            

             

              2
             
            
            
def get_sqrtW(x0,k):    #x0处其它样本点的权重
           
           

            

             

              3
             
            
            
    w=np.zeros(len(x))
           
           

            

             

              4
             
            
            
    for i in range(len(x)):
           
           

            

             

              5
             
            
            
        w[i]=exp(-(x[i]-x0)**2/(2*k*k))
           
           

            

             

              6
             
            
            
    w=np.array([sqrt(i) for i in w])
           
           

            

             

              7
             
            
            
    return w
           
          
         
        
       
      
      

      

       

       

      
     
    
   
  
 
 查看全部


 

  

   

    

     

      
     
 
     
 
      
 
       
 
        
 
         
 
          
 
           
 
           
 
            
 
           
 
          
 
          
 
          
 
          
 
          
 
           
 
           
 
          
 
          
 
           
 
            
 
             

               1
             
 
            
 
            
import statsmodels.api as sm
 
           
 
          
 
         
 
        
 
       
 
      
 
      
 
      
 
      
 
       
 
       
 
       
 
       
 
      
 
     
 
    
 
   
 
  
 
 
 
 查看全部
 
 
 
 
  
 
   
 刚才说到，ｋ是一个很关键的参数，我们从高斯函数的形式可以看出，ｋ取非常大的时候，每个样本点的权重都趋近于１，我们可以先取ｋ很大，检验一下是否正确
 
  
 
 
 
 
 
 
 
  
 
   
 
    
 
     
 
      
     
     

      

       

        

         

          

           

           

            

           
          
          

          

          

           

          
          

           

            

             

              1
             
            
            
#把整个过程定义成一个函数，方便改参数进行研究
           
           

            

             

              2
             
            
            
def get_yHat2(k):
           
           

            

             

              3
             
            
            
    yHat2=np.zeros(len(x))
           
           

            

             

              4
             
            
            
    for i in range(len(x)):
           
           

            

             

              5
             
            
            
        #把加权最小二乘转化为普通最小二乘，详情请看第二部分
           
           

            

             

              6
             
            
            
        w=get_sqrtW(x[i],k)
           
           

            

             

              7
             
            
            
        x2=w*x
           
           

            

             

              8
             
            
            
        x2=x2[x2>0]      #去掉样本权重为0的样本
           
           

            

             

              9
             
            
            
        y2=w*y
           
           

            

             

              10
             
            
            
        y2=y2[y2>0]
           
           

            

             

              11
             
            
            
        X=np.zeros((1,len(x2)))
           
           

            

             

              12
             
            
            
        X[0]=x2
           
           

            

             

              13
             
            
            
        X=X.T
           
           

            

             

              14
             
            
            
        X = sm.add_constant(X,has_constant='skip')
           
           

            

             

              15
             
            
            
        X[:,0]=w[w>0]
           
           

            

             

              16
             
            
            
        Y=y2
           
           

            

             

              17
             
            
            
        model = sm.OLS(Y, X)
           
           

            

             

              18
             
            
            
        results = model.fit()
           
           

            

             

              19
             
            
            
        a = results.params[0]
           
           

            

             

              20
             
            
            
        b = results.params[1]
           
           

            

             

              21
             
            
            
        yHat2[i]=a+b*x[i]      #得到xi处的估计值
           
           

            

             

              22
             
            
            
    return yHat2
           
          
         
        
       
      
      

      

       

       

      
     
    
   
  
 
 查看全部


 

  

   

    

     

      
     
 
     
 
      
 
       
 
        
 
         
 
          
 
           
 
           
 
            
 
           
 
          
 
          
 
          
 
          
 
          
 
           
 
           
 
          
 
          
 
           
 
            
 
             

               1
             
 
            
 
            
yHat2=get_yHat2(100000)  #ｋ取100000
 
           
 
           
 
            
 
             

               2
             
 
            
 
            
plt.figure(figsize=(12,6))
 
           
 
           
 
            
 
             

               3
             
 
            
 
            
plt.plot(x,yHat2,'r')
 
           
 
           
 
            
 
             

               4
             
 
            
 
            
plt.scatter(x,y)
 
           
 
          
 
         
 
        
 
       
 
      
 
      
 
      
 
      
 
       
 
       
 
       
 
       
 
      
 
     
 
    
 
   
 
  
 
 
 
 查看全部 
 
 
  
 
   

    <matplotlib.collections.PathCollection at 0xa96edd0>
   
 
   

    
   
 
  
 
 
 
 
 
 
 
  
 
   
 
    
 
     
 
      
     
     

      

       

        

         

          

           

           

            

           
          
          

          

          

           

          
          

           

            

             

              1
             
            
            
data=pd.DataFrame()
           
           

            

             

              2
             
            
            
data['y']=y
           
           

            

             

              3
             
            
            
data['yHatLinear']=yHatLinear
           
           

            

             

              4
             
            
            
data['yHat2']=yHat2
           
           

            

             

              5
             
            
            
data.head()
           
          
         
        
       
      
      

      

       

       

      
     
    
   
  
 
 查看全部
 

  

   

    

     
 
y
yHatLinear
yHat2
0
3.144993
2.325219
2.325219
1
7.656654
3.304751
3.304751
2
10.793202
4.284283
4.284283
3
13.490299
5.263815
5.263815
4
15.914143
6.243347
6.243347
    
   
  
 


 

  

   
可以看到用普通最小二乘估计出来的值和我们用局部加权估计出来的值非常的一致，说明我们的逻辑是对的。下面调整一下ｋ的值，来看看各种ｋ值下的拟合状况：
  
 


 

  

   

    

     

      
     
 
     
 
      
 
       
 
        
 
         
 
          
 
           
 
           
 
            
 
           
 
          
 
          
 
          
 
          
 
          
 
           
 
           
 
          
 
          
 
           
 
            
 
             

               1
             
 
            
 
            
yHat2=get_yHat2(10)  
 
           
 
           
 
            
 
             

               2
             
 
            
 
            
plt.figure(figsize=(12,6))
 
           
 
           
 
            
 
             

               3
             
 
            
 
            
plt.plot(x,yHat2,'r')
 
           
 
           
 
            
 
             

               4
             
 
            
 
            
plt.scatter(x,y)
 
           
 
           
 
            
 
             

               5
             
 
            
 
            
plt.title('k=10')
 
           
 
           
 
            
 
             

               6
             
 
            
 
            
​
 
           
 
           
 
            
 
             

               7
             
 
            
 
            
yHat2=get_yHat2(1)  
 
           
 
           
 
            
 
             

               8
             
 
            
 
            
plt.figure(figsize=(12,6))
 
           
 
           
 
            
 
             

               9
             
 
            
 
            
plt.plot(x,yHat2,'r')
 
           
 
           
 
            
 
             

               10
             
 
            
 
            
plt.scatter(x,y)
 
           
 
           
 
            
 
             

               11
             
 
            
 
            
plt.title('k=1')
 
           
 
           
 
            
 
             

               12
             
 
            
 
            
​
 
           
 
           
 
            
 
             

               13
             
 
            
 
            
yHat2=get_yHat2(0.1)  
 
           
 
           
 
            
 
             

               14
             
 
            
 
            
plt.figure(figsize=(12,6))
 
           
 
           
 
            
 
             

               15
             
 
            
 
            
plt.plot(x,yHat2,'r')
 
           
 
           
 
            
 
             

               16
             
 
            
 
            
plt.scatter(x,y)
 
           
 
           
 
            
 
             

               17
             
 
            
 
            
plt.title('k=0.1')
 
           
 
          
 
         
 
        
 
       
 
      
 
      
 
      
 
      
 
       
 
       
 
       
 
       
 
      
 
     
 
    
 
   
 
  
 
 
 
 查看全部 
 
 
  
 
   

    <matplotlib.text.Text at 0xb35ac50>
   
 
   

    
   
 
   

    
   
 
   

    
   
 
  
 
 
 
 
 
 
 
  
 
   
 可以看到，当ｋ越小时，拟合的效果越好，但是我们拟合的目的在于预测，需要避免过拟合的问题，这时候需要做Bias/Variance Trade-off：
 
   
 Cross Validation可以帮助我们做Bias/Variance Trade-off:
 
   
 一般的Validation就是把数据分为随机的两部分一部分做训练，一部分作验证。还有leave-one-out validation,k-fold Cross Validation等方法。
 
   
 训练的目的是为了让训练误差尽量减小，同时也要注意模型的自由度（参数个数）,避免测试误差很大。
 
   
 训练误差一般用ＭＳＥ来衡量：
 
   
 
    
     
      MSE=1n∑i(yi−y^i)2
     
   
 比如在我们这个简单的例子里，如果要对之后的数据进行预测，就需要通过Validation选择参数ｋ的大小，再对离需要预测的点最近的点做加权线性回归去估计那个点的值。 
   
 
  
 
 
 
 
 
 
 
  
 
   
 三.总结
 
   
 学以致用，统计的方法众多，运用存乎一心，学习的路上从来不孤单。本帖在于学习总结，如有谬误之处，望指正。
 
   
 给时光以生命，给岁月以文明。　　　——刘慈欣
 
  
 
 
 
 
 
 
 
  
 
   
 
    
 
     
 
      
     
 
     
 
      
 
       
 
        
 
         
 
          
 
           
 
           
 
            
 
           
 
          
 
          
 
          
 
          
 
          
 
           
 
           
 
          
 
          
 
           
 
            
 
             

               1
             
 
            
 
            
​