采用线性回归的方式，Y=bx+a,求出系数
#include 
#include 
int main()
{
//float x[6]={0.0,1.0,2.0,3.0,4.0,5.0};
//float y[6]={1.1,1.9,3.1,3.9,4.4,4.9};
/*float x[6]={0.0,1.0,2.0,3.0,4.0,5.0};
float y[6]={0.0,1.0,2.0,3.0,4.0,5.0};*/
float x[6]={0.0,1.0,2.0,3.0,4.0,5.0};
float y[6]={1.0,1.0,1.0,1.0,1.0,1.0};
float a,b,mxy,xx,yy,x2,x22;
int i;
a=b=mxy=xx=yy=x2=x22=0.0;
for(i=0;i<6;i++)
{
mxy=6.0*x[i]*y[i]+mxy;
xx=1.0*x[i]+xx;
yy=1.0*y[i]+yy;
x2=1.0*x[i]*x[i]*6.0+x2;
x22=1.0*x[i]+x22;
}
b=1.0*(mxy-xx*yy)/(x2-x22*x22);
a=1.0*yy/6.0-b*xx/6.0;
printf("Y=%0.2fx+%0.2f\n",b,a);
system("pause");
return 0;
}
Y=0.00x+1.00

为什么会用纯Python来实现呢？这种陶（hao）冶（wu）情（yong）操（chu）的做法肯定也不是我自己闲着蛋疼搞着玩的。。。  明明有更好的科学计算库不用，非要纯手写正面硬刚，除了老师用挂科作为威胁还能有谁？

下面是一大纯手打的堆理推导... 写的逻辑有些混乱，后续有时间再慢慢整理 觉得麻烦的小伙伴就不用看了，反正我代码里也写的有注释，代码稍后传上github后更新连接

---------------------------------------------------------------------------------------------------------------

 

多元线性回归 最小二乘法实现

 

线性回归包括一元线性回归和多元线性回归，一元的是只有一个x和一个y。多元的是指有多个x和一个y。

 

 

对于平面中的这n个点，可以使用无数条曲线来拟合。要求样本回归函数尽可能好地拟合这组值。综合起来看，这条直线处于样本数据的中心位置最合理。 选择最佳拟合曲线的标准可以确定为：使总的拟合误差（即总残差）达到最小。有以下三个标准可以选择：

        （1）用“残差和最小”确定直线位置是一个途径。但很快发现计算“残差和”存在相互抵消的问题。

        （2）用“残差绝对值和最小”确定直线位置也是一个途径。但绝对值的计算比较麻烦。

        （3）最小二乘法的原则是以“残差平方和最小”确定直线位置。用最小二乘法除了计算比较方便外，得到的估计量还具有优良特性。这种方法对异常值非常敏感。

 

样本x为一个多维向量，x=[x1,x2,…xd]

那么权重W也是与X同维度的一个向量，W=[w1,w2,…wd]

所以Y对x的预测可以写为Y=W ∙ x+b(这里为内积的形式)，但通常将b移到W向量里面去，写为W=[b,w1,w2,…wd]，改变后的W用θ表示，θ=[θ0,θ1,θ2,…θd]

因为W的维度变了所以我们也要给x对应的多加一个维度，此时的x为[x0,x1,…xd]，使x0等于1，让其在矩阵相乘时匹配刚才在W中加入的b

现在Y=W ∙ x+b 就等价为 Y=θ∙x

在机器学习的矩阵运算当中，有个约定俗称的规定，当说某某向量时都指的是一个列向量，所以θ和x都是一个d+1维的列向量[θ0⋮θd]、[x0⋮xd]，同时也都是一个(d+1)x1维的矩阵

这里的维度就对应方程里面的元数，有几维就是几元

所以这里用大写的X表示有m个样本，X=[⋯x1T⋯⋯⋮⋯⋯xmT⋯]=[x01⋯xd1⋮⋯⋮x0m⋯xdm]，（这里的X是一个矩阵，包含了所有的样本点x1,x2,…xm而每一个样本点都是有相同的多维的， [⋯x1T⋯⋯⋮⋯⋯xmT⋯]* [θ0⋮θd]= y0⋮ym=Y，得到一个m*1维的矩阵Y，这个矩阵中的每一个参数代表给定的θ参数下对每一个样本xi的预测值yi，Y是对应样本的标签值，也是一个列向量。

得到Y（预测值）过后，就可以算损失函数了，L=i=1m（yi-yi）2，y是对应样本的真实值

L=（y1-y1）2+（y2-y2）2+…+（ym-ym）2

 =[ y1-y1,⋯,ym-ym] ∙ [y1-y1⋮ym-ym]

 =(y1,⋯,ym-[y1,⋯,ym]) ∙ [y1-y1⋮ym-ym]

 =(x∙θ-Y)T ∙  (x∙θ-Y)

因为x、θ和Y都是一个m+1维的列向量所以第一项(x∙θ-Y)的结果也是一个列项量，所以加了一个转置

得到损失函数L=(x∙θ-Y)T (x∙θ-Y)，通过改变θ的值来最小化L

对θ求偏导，∂L∂θ=∂(θTxT-YT)(xθ-Y)∂θ=θTxTxθ-θTxTY-YTxθ+YTY∂θ，每一项都对θ求偏导然后加起来

 ∂L∂θ=2xTxθ-xTY-xTY，要使损失函数最小，就是偏导为0的点

所以推出xTxθ=xTY，θ=( xTx)-1xTY，这就要求(xTx)是可逆的

 

最小二乘法与梯度下降法的区别：

实现方法和结果不同：最小二乘法是直接对∆求导找出全局最小，是非迭代法。而梯度下降法是一种迭代法，先给定一个β ，然后向∆下降最快的方向调整β ，在若干次迭代之后找到局部最小。梯度下降法的缺点是到最小点的时候收敛速度变慢，并且对初始点的选择极为敏感，其改进大多是在这两方面下功夫。

 

 

 

 

 

 

 

 

 
Cholesky分解求系数参考：
[1]冯天祥. 多元线性回归最小二乘法及其经济分析[J]. 经济师,2003,11:129.
 
还可以采用最小二乘法来估计参数：







 
算法设计也可以参考两种系数最终公式设计。
下面的Java代码由网友设计，采用第一种方法计算参数。


  1 /**
  2  * 多元线性回归分析
  3  *
  4  * @param x[m][n]
  5  *            每一列存放m个自变量的观察值
  6  * @param y[n]
  7  *            存放随即变量y的n个观察值
  8  * @param m
  9  *            自变量的个数
 10  * @param n
 11  *            观察数据的组数
 12  * @param a
 13  *            返回回归系数a0,...,am
 14  * @param dt[4]
 15  *            dt[0]偏差平方和q,dt[1] 平均标准偏差s dt[2]返回复相关系数r dt[3]返回回归平方和u
 16  * @param v[m]
 17  *            返回m个自变量的偏相关系数
 18  */
 19 public static void sqt2(double[][] x, double[] y, int m, int n, double[] a,
 20         double[] dt, double[] v) {
 21     int i, j, k, mm;
 22     double q, e, u, p, yy, s, r, pp;
 23     double[] b = new double[(m + 1) * (m + 1)];
 24     mm = m + 1;
 25     b[mm * mm - 1] = n;
 26     for (j = 0; j <= m - 1; j++) {
 27         p = 0.0;
 28         for (i = 0; i <= n - 1; i++)
 29             p = p + x[j][i];
 30         b[m * mm + j] = p;
 31         b[j * mm + m] = p;
 32     }
 33     for (i = 0; i <= m - 1; i++)
 34         for (j = i; j <= m - 1; j++) {
 35             p = 0.0;
 36             for (k = 0; k <= n - 1; k++)
 37                 p = p + x[i][k] * x[j][k];
 38             b[j * mm + i] = p;
 39             b[i * mm + j] = p;
 40         }
 41     a[m] = 0.0;
 42     for (i = 0; i <= n - 1; i++)
 43         a[m] = a[m] + y[i];
 44     for (i = 0; i <= m - 1; i++) {
 45         a[i] = 0.0;
 46         for (j = 0; j <= n - 1; j++)
 47             a[i] = a[i] + x[i][j] * y[j];
 48     }
 49     chlk(b, mm, 1, a);
 50     yy = 0.0;
 51     for (i = 0; i <= n - 1; i++)
 52         yy = yy + y[i] / n;
 53     q = 0.0;
 54     e = 0.0;
 55     u = 0.0;
 56     for (i = 0; i <= n - 1; i++) {
 57         p = a[m];
 58         for (j = 0; j <= m - 1; j++)
 59             p = p + a[j] * x[j][i];
 60         q = q + (y[i] - p) * (y[i] - p);
 61         e = e + (y[i] - yy) * (y[i] - yy);
 62         u = u + (yy - p) * (yy - p);
 63     }
 64     s = Math.sqrt(q / n);
 65     r = Math.sqrt(1.0 - q / e);
 66     for (j = 0; j <= m - 1; j++) {
 67         p = 0.0;
 68         for (i = 0; i <= n - 1; i++) {
 69             pp = a[m];
 70             for (k = 0; k <= m - 1; k++)
 71                 if (k != j)
 72                     pp = pp + a[k] * x[k][i];
 73             p = p + (y[i] - pp) * (y[i] - pp);
 74         }
 75         v[j] = Math.sqrt(1.0 - q / p);
 76     }
 77     dt[0] = q;
 78     dt[1] = s;
 79     dt[2] = r;
 80     dt[3] = u;
 81 }
 82 private static int chlk(double[] a, int n, int m, double[] d) {
 83     int i, j, k, u, v;
 84     if ((a[0] + 1.0 == 1.0) || (a[0] < 0.0)) {
 85         System.out.println("fail\n");
 86         return (-2);
 87     }
 88     a[0] = Math.sqrt(a[0]);
 89     for (j = 1; j <= n - 1; j++)
 90         a[j] = a[j] / a[0];
 91     for (i = 1; i <= n - 1; i++) {
 92         u = i * n + i;
 93         for (j = 1; j <= i; j++) {
 94             v = (j - 1) * n + i;
 95             a[u] = a[u] - a[v] * a[v];
 96         }
 97         if ((a[u] + 1.0 == 1.0) || (a[u] < 0.0)) {
 98             System.out.println("fail\n");
 99             return (-2);
100         }
101         a[u] = Math.sqrt(a[u]);
102         if (i != (n - 1)) {
103             for (j = i + 1; j <= n - 1; j++) {
104                 v = i * n + j;
105                 for (k = 1; k <= i; k++)
106                     a[v] = a[v] - a[(k - 1) * n + i] * a[(k - 1) * n + j];
107                 a[v] = a[v] / a[u];
108             }
109         }
110     }
111     for (j = 0; j <= m - 1; j++) {
112         d[j] = d[j] / a[0];
113         for (i = 1; i <= n - 1; i++) {
114             u = i * n + i;
115             v = i * m + j;
116             for (k = 1; k <= i; k++)
117                 d[v] = d[v] - a[(k - 1) * n + i] * d[(k - 1) * m + j];
118             d[v] = d[v] / a[u];
119         }
120     }
121     for (j = 0; j <= m - 1; j++) {
122         u = (n - 1) * m + j;
123         d[u] = d[u] / a[n * n - 1];
124         for (k = n - 1; k >= 1; k--) {
125             u = (k - 1) * m + j;
126             for (i = k; i <= n - 1; i++) {
127                 v = (k - 1) * n + i;
128                 d[u] = d[u] - a[v] * d[i * m + j];
129             }
130             v = (k - 1) * n + k - 1;
131             d[u] = d[u] / a[v];
132         }
133     }
134     return (2);
135 }
136 /**
137  * @param args
138  */
139 public static void main(String[] args) {
140     // TODO Auto-generated method stub
141     /**
142      * 一元回归
143      */
144     // int i;
145     // double[] dt=new double[6];
146     // double[] a=new double[2];
147     // double[] x={ 0.0,0.1,0.2,0.3,0.4,0.5,
148     // 0.6,0.7,0.8,0.9,1.0};
149     // double[] y={ 2.75,2.84,2.965,3.01,3.20,
150     // 3.25,3.38,3.43,3.55,3.66,3.74};
151     // SPT.SPT1(x,y,11,a,dt);
152     // System.out.println("");
153     // System.out.println("a="+a[1]+" b="+a[0]);
154     // System.out.println("q="+dt[0]+" s="+dt[1]+" p="+dt[2]);
155     // System.out.println(" umax="+dt[3]+" umin="+dt[4]+" u="+dt[5]);
156     /**
157      * 多元回归
158      */
159     int i;
160     double[] a = new double[4];
161     double[] v = new double[3];
162     double[] dt = new double[4];
163     double[][] x = { { 1.1, 1.0, 1.2, 1.1, 0.9 },
164             { 2.0, 2.0, 1.8, 1.9, 2.1 }, { 3.2, 3.2, 3.0, 2.9, 2.9 } };
165     double[] y = { 10.1, 10.2, 10.0, 10.1, 10.0 };
166     SPT.sqt2(x, y, 3, 5, a, dt, v);
167     for (i = 0; i <= 3; i++)
168         System.out.println("a(" + i + ")=" + a[i]);
169     System.out.println("q=" + dt[0] + "  s=" + dt[1] + "  r=" + dt[2]);
170     for (i = 0; i <= 2; i++)
171         System.out.println("v(" + i + ")=" + v[i]);
172     System.out.println("u=" + dt[3]);
173 }

View Code



下面的C++代码由网友提供，采用第二中方法计算系数。

#include<iostream>
#include<fstream>
#include<iomanip>
using namespace std;
void transpose(double **p1,double **p2,int m,int n);
void multipl(double **p1,double **p2,double **p3,int m,int n,int p);
void Inver(double **p1,double **p2,int n);
double SD(double **p1,double **p2,double **p3,double **p4,int m,int n);
double ST(double **p1,int m);
void de_allocate(double **data,int m);
int main() {
int row,col;
char filename[30];
double SDsum,STsum,F,R2;
cout<<"Input original data file: \n";
ifstream infile;  //打开文件
cin>>filename;
infile.open(filename);
if(!infile) {
cout<<"Opening the file failed!\n";
exit(1);
}
infile>>row>>col; //读入文件中的行数和列数
double **matrix=new double*[row]; //为动态二维数组分配内存
double **X=new double*[row];
double **Y=new double*[row];
double **XT=new double*[col];
double **XTX=new double*[col];
double **XTXInv=new double*[col];
double **XTXInvXT=new double*[col];
double **B=new double*[col];
double **YE=new double*[row];
for(int i=0;i<row;i++) {
  matrix[i]=new double[col];
  X[i]=new double[col];
  Y[i]=new double[1];
  Y[i]=new double[1];
  YE[i]=new double[1];
}
for(int i=0;i<col;i++) {
  XT[i]=new double[row];
  XTX[i]=new double[2*col];/为什么必须分配2*col列空间而不是col？在矩阵求逆时，XTX变增广矩阵，列数变为原来2吧倍，跟求逆算法有关。
  XTXInv[i]=new double[col];
  XTXInvXT[i]=new double[row];
  B[i]=new double[1];
}
for(int i=0;i<row;i++)
  for(int j=0;j<col;j++)
    infile>>matrix[i][j];
infile.close();
for(int i=0;i<row;i++) { //提取1X和Y数组列
  X[i][0]=1;
  Y[i][0]=matrix[i][col-1];
  for(int j=0;j<col-1;j++)
    X[i][j+1]=matrix[i][j];
}
transpose(X,XT,row,col);
multipl(XT,X,XTX,col,row,col);
Inver(XTX,XTXInv,col);
multipl(XTXInv,XT,XTXInvXT,col,col,row);
multipl(XTXInvXT,Y,B,col,row,1);
SDsum=SD(Y,X,B,YE,row,col);
STsum=ST(Y,row);
F=((STsum-SDsum)/(col-1))/(SDsum/(row-col));
R2=1/(1+(row-col)/F/(col-1));
cout<<"输出B:\n";  //屏幕输出结果B，SD，ST，F，R2
for(int i=0;i<col;i++)
  cout<<setiosflags(ios::fixed)<<setprecision(4)<<B[i][0]<<' ';
  cout<<endl;
cout<<"SD="<<SDsum<<';'<<"ST="<<STsum<<';'<<"F="<<F<<';'<<"R2="<<R2<<endl;
ofstream outfile; // 结果写入文件
cout<<"Output file'name:\n";
cin>>filename;
outfile.open(filename);
if(!outfile) {
  cout<<"Opening the file failed!\n";
  exit(1);
}
outfile<<"输出B:\n";
for(int i=0;i<col;i++)
  outfile<<B[i][0]<<' ';
  outfile<<endl;
outfile<<setiosflags(ios::fixed)<<setprecision(4)<<"SD="<<SDsum<<';'<<"ST="<<STsum<<';'<<"F="<<F<<';'<<"R2="<<R2<<endl;
outfile<<"Y and YE and Y-YE's value are:\n";
for(int i=0;i<row;i++)
    outfile<<Y[i][0]<<"  "<<YE[i][0]<<"  "<<Y[i][0]-YE[i][0]<<endl;
outfile.close();
de_allocate(matrix,row);
de_allocate(X,row);
de_allocate(Y,row);
de_allocate(XT,col);
de_allocate(XTX,col);
de_allocate(XTXInv,col);
de_allocate(XTXInvXT,col);
de_allocate(B,col);
de_allocate(YE,row);
system("pause");
return(0);
}
void de_allocate(double **data,int m) { //释放内存单元
  for(int i=0;i<m;i++)
    delete []data[i];
  delete []data;
}
double ST(double **p1,int m) { //求总离差平方和ST
  double sum1=0,sum2=0,Yave=0;
  for(int i=0;i<m;i++)
    sum1+=p1[i][0];
  Yave=sum1/m;
  for(int i=0;i<m;i++)
    sum2+=(p1[i][0]-Yave)*(p1[i][0]-Yave);
  return sum2;
}
double SD(double **p1,double **p2,double **p3,double **p4,int m,int n) { //求偏差平方和SD
  double sum1=0,sum2=0;
  for(int i=0;i<m;i++) {
    sum1=0;
    for(int k=0;k<n;k++)
      sum1+=p2[i][k]*p3[k][0];
    p4[i][0]=sum1;
  }
  for(int i=0;i<m;i++)
    sum2+=(p1[i][0]-p4[i][0])*(p1[i][0]-p4[i][0]);
  return sum2;
}
void transpose(double **p1,double **p2,int m,int n) {  //矩阵转置
for(int i=0;i<n;i++)
  for(int j=0;j<m;j++)
      p2[i][j]=p1[j][i];
}
void multipl(double **p1,double **p2,double **p3,int m,int n,int p) { //矩阵相乘
double sum;
for(int i=0;i<m;i++) {
  for(int j=0;j<p;j++) {
    sum=0;
    for(int k=0;k<n;k++)
      sum+=p1[i][k]*p2[k][j];
    p3[i][j]=sum;
  }
}
}
void Inver(double **p1,double **p2,int n) {  //求逆矩阵
//初始化矩阵在右侧加入单位阵
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++) {
      p1[i][j+n]=0;
      p1[i][i+n]=1;
    }
  }
//对于对角元素为0的进行换行操作
  for(int i=0;i<n;i++)
   {
    while(p1[i][i]==0)
     {
      for(int j=i+1;j<n;j++)
        {
         if (p1[j][i]!=0)
         { double temp=0;
          for(int r=i;r<2*n;r++)
            {temp=p1[j][r];p1[j][r]=p1[i][r];p1[i][r]=temp;}
         }
         break;
        }
     }
    //if (p1[i][i]==0) return 0;
   }
 //行变换为上三角矩阵
  double k=0;
  for(int i=0;i<n;i++) {
    for(int j=i+1;j<n;j++) {
    k=(-1)*p1[j][i]/p1[i][i];
    for(int r=i;r<2*n;r++)
      p1[j][r]+=k*p1[i][r];
    }
  }
  //行变换为下三角矩阵
  //double k=0;
  for(int i=n-1;i>=0;i--) {
    for(int j=i-1;j>=0;j--) {
      k=(-1)*p1[j][i]/p1[i][i];
      for(int r=0;r<2*n;r++)
        p1[j][r]+=k*p1[i][r];
      }
  }
  //化为单位阵
  for(int i=n-1;i>=0;i--) {
    k=p1[i][i];
    for(int j=0;j<2*n;j++)
      p1[i][j]/=k;
  }
  //拆分出逆矩阵
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++)
      p2[i][j]=p1[i][n+j];
  }
}

 


二元线性回归最小二乘法拟合空间直线。网友提供


#include "stdafx.h"
using namespace std;
using std::cout;
using std::cin;
using std::endl;
#include<math.h>
#include <windows.h>
/*
这是一个控制台程序，任何一个空间直线方程都能以如下的方式表示
x=az+b
y=cz+d
z=z
即
(x-b)/a=(y-d)/c=z/1
*/
#include <vector>
using std::vector;
/* the fitting vaule of this line are: a=2 b=3 c=3 d=0
double z[10]={0.5,0.7,1,1.2,1.5,1.8,2,2.5,2.8,3};
double y[10]={1.5,2.1,3,3.6,4.5,5.4,6,7.5,8.4,9};
double x[10]={4,4.4,5,5.4,6,6.6,7,8,8.6,9};
*/
double x[]={1,1.5,2,2.5,3,3.5,4,4.5,5};
double y[]={-8.1,-7.2,-6.2,-5.5,-4.8,-3.8,-3,-2.2,-1.3};
double z[]={-12,-11.8,-10.7,-9.5,-8.2,-7,-6,-4.5,-3.5};
int _tmain(int argc, _TCHAR* argv[])
{
     
    vector<double>position_z;
    vector<double>position_y;
    vector<double>position_x;
    if ((sizeof(z)/sizeof(double))!=(sizeof(x)/sizeof(double))||(sizeof(z)/sizeof(double))!=(sizeof(y)/sizeof(double))||(sizeof(x)/sizeof(double))!=(sizeof(y)/sizeof(double)))
    {
        ::MessageBox(NULL,"请检查输入数组的长度是否相等","方程无解!",MB_OK);
    }
    //int ss=sizeof(z)/sizeof(double);
    for (int i=0;i<(sizeof(z)/sizeof(double));++i)
    {
        position_z.push_back(z[i]);
        position_y.push_back(y[i]);
        position_x.push_back(x[i]);
    }
         
     
        
    //double m = position_z.size();
    //caculate o1,p1,c1,o2,p2,c2
    double o1 = 0.0;
    double p1 = 0.0;
    double x1 = 0.0;
    double o2 = 0.0;
    double p2 = 0.0;
    double x2 = 0.0;
    double y1 = 0.0;
    double y2 = 0.0;
    for (vector<double>::size_type t=0;t< position_z.size();++t)
    {
        o1 += position_z[t]*position_z[t];
        p1 += position_z[t];
        x1 += position_x[t]*position_z[t];
        o2 += position_z[t];
        p2  = position_z.size();
        x2 += position_x[t] ;
        
        y2 += position_y[t];
        y1 += position_y[t]*position_z[t];
    }
    //caculate a b
    double a = 0.0 ,b = 0.0, c = 0.0, d = 0.0 ;
    if ((o1*p2-o2*p1)!= 0)
    {
        a = (x1*p2 - x2*p1) /(o1*p2-o2*p1);
        if (a<1e-10)
        {
            a=0;
        }
        b = (x2*o1 - x1*o2) /(o1*p2-o2*p1);
        if (b<1e-10)
        {
            b=0;
        }
        c = (y1*p2 - y2*p1) /(o1*p2-o2*p1);
        if (c<1e-10)
        {
            c=0;
        }
        d = (y2*o1 - y1*o2) /(o1*p2-o2*p1);
        if (d<1e-10)
        {
            d=0;
        }
    }
    else
    {
        ::MessageBox(NULL,"OK","方程无解!",MB_OK);
    }
     
   cout<<a<<"  " <<b<<"  "<<c<<"  "<<d<<endl;
    system ("pause");
    return 0;
}


　　


来自为知笔记(Wiz)



转载于:https://www.cnblogs.com/star91/p/4750102.html

线性回归是确定两种及两种以上变量的相互依赖关系。在数据分析中，线性回归是最简单且最有效的分析方法。举个简单的例子，某商品的利润在售价为2元、5元、10元时分别为4元、10元、20元，我们很容易得出商品的利润与售价的关系符合直线：.在上面这个简单的一元线性回归方程中，我们称“2”为回归系数，即斜率为其回归系数，回归系数表示商品的售价(x)每变动一个单位，其利润(y)与之对应的变动关系。

1.线性回归

                                                    



如上图所示，上图为一个简单的一元线性回归示意图，线性回归表示这些离散的点总体上“最逼近”哪条直线。同理，在二元线性回归中，离散的点在空间中最逼近哪个平面。其它多元线性回归亦是如此。线性回归在数据量化分析，数据预测方面有着重要应用。

2.最小二乘法

同样是上图，那这些离散的点总体上“最逼近”哪条直线呢，这需要一个量化。对于同一x，实际数据为y，回归方程的推断数据为ax，则误差d=y-ax。最小二乘法定义为，当，D最小时，回归方程的拟合度最高，根据求出a即可，对最小二乘法的定义方程进行微分求极值即可得出a。

对于多元变量，设

                                        ,    ,    

其回归方程为Y=Xa。同理，现在我们要求矩阵a的值，使得

                                                               

取最小值即可。对上式进行微分求解(具体求解步骤可参考线性代数相关知识)，可得

                                                                              

如果该矩阵满秩，则

                                                                           

上述便为最小二乘法的基本思想。

3.python实现

一元线性回归测试：

from numpy.linalg import inv  # 矩阵求逆
from numpy import dot  # 矩阵点乘
from numpy import mat  # 二维矩阵

X = mat([1, 2, 3]).reshape(3, 1)  # x为1,2,3
Y = mat([5, 10, 15]).reshape(3, 1)  # y为5,10,15
a = dot(dot(inv(dot(X.T, X)), X.T), Y)  # 最小二乘法公式
print(a)

首先，这里推荐一个公式的代码实现技巧，代码从内向外写，即X的转置与X点乘，再求逆， 再与X的转置点乘，再与Y点乘，一层层加括号即可。

由上述代码不难看出，y与x的关系为y=5x，下面我们看测试结果是否为5：



测试结果正确。

实战：

这里准备了一组数据，数据示例截图如下：



 

上述数据为某一商品的销售量与售价、服务投资和其它投资的对应关系，即我们要求出：

                                                     ，    ,      

中，矩阵a。

代码如下：

import numpy as np
import pandas as pd
from numpy.linalg import inv  # 矩阵求逆
from numpy import dot  # 矩阵点乘


dataset = pd.read_csv('C:\\Users\\57105\\Desktop\\data.csv')  # 读入数据
X = dataset.iloc[:, 2: 5]  # x为所有行，2到4列
Y = dataset.iloc[:, 1]  # y为所有行，第1列
a = dot(dot(inv(np.dot(X.T, X)), X.T), Y)  # 最小二乘法求解公式
print(a)

测试结果：



即，对于单条数据，即，该回归方程反应了商品销售量与售价、服务投资和其它投资的关系。

上述便为线性回归的最小二乘法求解方式，关于线性回归的另一种求解方式—梯度下降法可见笔者相关文章。