逻辑斯谛回归模型

时间:2021-01-20 23:51:34

    逻辑斯谛回归模型是研究因变量为二分类或多分类观察结果与影响因素之间的关系的一种概率型非线性回归模型。逻辑斯谛回归系数通过最大似然估计得到。Logistic函数如下:

    逻辑斯谛回归模型     

    式中x

    逻辑斯谛回归模型     

    这里逻辑斯谛回归模型是输入变量的n个特征,然后按照Logistic函数形式求出。

    假设有n个独立变量的向量逻辑斯谛回归模型,设条件概率逻辑斯谛回归模型x条件下y发生的概率(假设y=1y发生)。则Logistic函数表示为:

    逻辑斯谛回归模型逻辑斯谛回归模型     

    同理,在x条件下y不发生的概率为:

逻辑斯谛回归模型

    Logistic回归都是围绕Logistic函数展开的,如何求解逻辑斯谛回归模型Logistic回归模型的主要问题,采用的最大似然估计来求解这组参数。

    假设有m个观测样本,观测值分别为逻辑斯谛回归模型,设逻辑斯谛回归模型为给定条件下逻辑斯谛回归模型的概率,同理逻辑斯谛回归模型的概率为逻辑斯谛回归模型,得到一个观测值得概率为:

    逻辑斯谛回归模型逻辑斯谛回归模型     

    因为各观测样本间相互独立,于是得到似然函数:

    逻辑斯谛回归模型逻辑斯谛回归模型    

    对似然函数取对数:

    逻辑斯谛回归模型逻辑斯谛回归模型    

    现要求向量逻辑斯谛回归模型使的逻辑斯谛回归模型最大,其中:

    逻辑斯谛回归模型逻辑斯谛回归模型    

    要求的最大似然估计,我们需要确定似然函数存在局部极大值。因此,对似然函数求偏导后得:

    逻辑斯谛回归模型     

    由多元函数极值的必要条件可知,若多元函数在一点取得极值,且一阶偏导存在,则该点处所有一阶偏导为0。由此,可以得出n+1个方程,如下:

    逻辑斯谛回归模型     

    由此方程解出的逻辑斯谛回归模型 不一定是似然函数的极值,需要通过Hessian矩阵来判断得出的解是否为似然函数的极值。

    Hessian矩阵是一个多元函数的二阶偏导构成的方阵,描述了函数的局部曲率。对一个多元函数逻辑斯谛回归模型 ,如果他的二阶偏导都存在,那么Hessian矩阵如下:

    逻辑斯谛回归模型     

    通过Hessian矩阵,我们可以判断一点M处极值的三种情况:

  1. 如果逻辑斯谛回归模型 是正定矩阵,则临界点M处是一个局部极小值;
  2. 如果逻辑斯谛回归模型 是负定矩阵,则临界点M处是一个局部极大值;
  3. 如果逻辑斯谛回归模型 是不定矩阵,则临界点M处不是极值。

对于中的n+1个方程,要求Hessian矩阵,先要求似然函数的二阶偏导,即:

    逻辑斯谛回归模型     

    则似然函数的Hessian矩阵为

    逻辑斯谛回归模型     

    设有矩阵X、A:

    逻辑斯谛回归模型     

    逻辑斯谛回归模型     

    则似然函数的Hessian矩阵可表示为:

    逻辑斯谛回归模型     

    显然,矩阵A是负定的,则可以证明H也是负定的,说明似然函数存在局部极大值。因此,可以使用牛顿迭代法(Newton's Method)来求逻辑斯谛回归模型

    对一元函数,使用牛顿迭代法来求零点。假设要求逻辑斯谛回归模型 的解,首先选取一个点逻辑斯谛回归模型 作为迭代起点,通过下面的式子进行迭代,直到达到指定精度为止。

    逻辑斯谛回归模型     

    由此,有时起始点选择很关键,如果函数只存在一个零点,那么这个起始点选取就无关重要。对已Logistic回归问题,Hessian矩阵对于任意数据都是负定的,所以说极值点只有一个,初始点的选取无关紧要。

    因此,对于上述Logistic回归的似然函数, 令:

    逻辑斯谛回归模型     

    则由可以得到如下的迭代式子:

    逻辑斯谛回归模型     

    由于Hessian矩阵是负定的,将矩阵A提取一个负号,得:

    逻辑斯谛回归模型     

    然后Hessian矩阵变为

    逻辑斯谛回归模型     

    这样,Hessian矩阵就是对称正定的了。那么牛顿迭代式变为:

    逻辑斯谛回归模型     

    现在,关键是如何快速并有效的计算逻辑斯谛回归模型,即解方程组逻辑斯谛回归模型 。由于逻辑斯谛回归模型 是对称正定的,可以使用Cholesky矩阵分解法来解。

    若逻辑斯谛回归模型 对称正定,则存在一个对角元为正数的下三角矩阵逻辑斯谛回归模型,使得逻辑斯谛回归模型 成立。对于逻辑斯谛回归模型,可以通过以下步骤求解:

  1. 逻辑斯谛回归模型 的Cholesky分解,得到逻辑斯谛回归模型
  2. 求解逻辑斯谛回归模型 ,得到逻辑斯谛回归模型
  3. 求解逻辑斯谛回归模型 ,得到逻辑斯谛回归模型

现在的关键问题是对逻辑斯谛回归模型 进行Cholesky分解。假设:

    逻辑斯谛回归模型     

    通过逻辑斯谛回归模型比较两边的关系,首先由

    逻辑斯谛回归模型     

再由

    逻辑斯谛回归模型     

    这样,得到了矩阵逻辑斯谛回归模型 的第一列元素。假设,已经算出了逻辑斯谛回归模型 的前逻辑斯谛回归模型 列元素,通过

    逻辑斯谛回归模型     

    可以得出

    逻辑斯谛回归模型     

    进一步由

    逻辑斯谛回归模型     

    最终:

    逻辑斯谛回归模型     

    这样便通过逻辑斯谛回归模型 的前逻辑斯谛回归模型 列求出了第逻辑斯谛回归模型 列,一直递推下去即可求出逻辑斯谛回归模型 。这种方法称为平方根法。

    利用上述方法需要进行开方,这有可能损失精度和增加运算量,为了避免开方,将Cholesky分解进行改进,即:

    逻辑斯谛回归模型     

    其中:逻辑斯谛回归模型 是单位下三角矩阵,逻辑斯谛回归模型 为对角均为正数的对角矩阵。把这一分解叫逻辑斯谛回归模型分解。设:

    逻辑斯谛回归模型     

    则对于逻辑斯谛回归模型,求解步骤变为:

  1. 逻辑斯谛回归模型逻辑斯谛回归模型分解,得到逻辑斯谛回归模型
  2. 求解逻辑斯谛回归模型 ,得到逻辑斯谛回归模型
  3. 求解逻辑斯谛回归模型 ,得到逻辑斯谛回归模型

对比两边元素,可以得到:

    逻辑斯谛回归模型     

    由此可以确定逻辑斯谛回归模型逻辑斯谛回归模型 的公式如下:

    逻辑斯谛回归模型     

    逻辑斯谛回归模型     

    逻辑斯谛回归模型     

牛顿迭代法

  1. using System;
  2. using System.Collections.Generic;
  3. using System.Text;
  4. using Model.MatrixDecomposition;
  5. using Model.Matrix;
  6. using System.IO;
  7.  
  8. namespace Model.NewtonRaphson
  9. {
  10.     internal class NewtonRaphsonIterate
  11.     {
  12.         Cholesky_LDL_Decomposition decompositionMatrixByLDL = new Cholesky_LDL_Decomposition();
  13.  
  14.         private double[,] matrix_Hessian = null;
  15.         private double[,] matrix_X = null;
  16.         private double[,] matrix_A = null;
  17.  
  18.  
  19.         private double[,] vector_HU = null;
  20.         private double[,] vector_U = null;
  21.         private double[,] vector_Y = null;
  22.         private double[,] vector_Omega = null;
  23.         private double[,] vector_PI = null;
  24.         private double[,] old_VectorOmega = null;
  25.  
  26.         #region 属性
  27.         public double[,] MatrixL
  28.         {
  29.             get
  30.             {
  31.                 return decompositionMatrixByLDL.MatrixL;
  32.             }
  33.         }
  34.  
  35.  
  36.         public double[,] MatrixD
  37.         {
  38.             get
  39.             {
  40.                 return decompositionMatrixByLDL.MatrixD;
  41.             }
  42.         }
  43.  
  44.         public double[,] Hessian
  45.         {
  46.             get
  47.             {
  48.                 return this.matrix_Hessian;
  49.             }
  50.         }
  51.         #endregion
  52.  
  53.         /// <summary>
  54.         /// 执行牛顿迭代法
  55.         /// </summary>
  56.         /// <param name="matrix_X">输入矩阵</param>
  57.         /// <param name="vector_Y">输出向量</param>
  58.         /// <param name="error_Threashold">迭代停止阈值</param>
  59.         /// <param name="omega">计算完成的参数</param>
  60.         internal void Run(double[,] matrix_X, double[,] vector_Y, double error_Threashold, ref double[,] omega)
  61.         {
  62.             double error = 1;
  63.             this.matrix_X = matrix_X;
  64.             this.vector_Y = vector_Y;
  65.             this.vector_Omega = omega;
  66.             InitializeParameter(matrix_X.GetLength(0));
  67.             int i = 0;
  68.             while (error > error_Threashold)
  69.             {
  70.                 i++;
  71.                 error = 0;
  72.                 old_VectorOmega = (double[,])vector_Omega.Clone();
  73.                 GetMatrixAAndVectorPI();
  74.                 matrix_Hessian = MatrixOperation.MultipleMatrix(
  75.                     MatrixOperation.MatrixMultiDiagMatrix(
  76.                     MatrixOperation.TransportMatrix(matrix_X), matrix_A),
  77.                     matrix_X);
  78.                 GetMatrixU();
  79.                 decompositionMatrixByLDL.Cholesky((double[,])matrix_Hessian.Clone(), vector_U, ref vector_HU);
  80.                 vector_Omega = MatrixOperation.AddMatrix(vector_Omega, vector_HU);
  81.                 GetIterationError(ref error);
  82.                 //TODO:迭代计算
  83.             }
  84.             omega = (double[,])vector_Omega.Clone();
  85.         }
  86.  
  87.         private void InitializeParameter(int rowNumber)
  88.         {
  89.             matrix_A = new double[rowNumber, 1];
  90.             vector_PI = new double[rowNumber, 1];
  91.         }
  92.  
  93.         /// <summary>
  94.         /// 获取H=X^TAX的A以及PI(Xi)
  95.         /// </summary>
  96.         private void GetMatrixAAndVectorPI()
  97.         {
  98.             for (int i = 0; i < matrix_X.GetLength(0); i++)
  99.             {
  100.                 double temp_A = 0;
  101.                 //matrix_A[i, 0] += vector_Omega[0, 0];
  102.                 for (int j = 0; j < matrix_X.GetLength(1); j++)
  103.                 {
  104.                     temp_A += matrix_X[i, j] * vector_Omega[j, 0];
  105.                 }//for2
  106.                 matrix_A[i, 0] = (1.0) / (1 + Math.Exp(-temp_A));
  107.                 vector_PI[i, 0] = matrix_A[i, 0];
  108.                 matrix_A[i, 0] *= (1 - matrix_A[i, 0]);
  109.  
  110.             }//for1
  111.         }
  112.  
  113.         //计算HU中的U
  114.         //U=X^T(Y-PI)
  115.         private void GetMatrixU()
  116.         {
  117.             vector_U = MatrixOperation.MultipleMatrix(MatrixOperation.TransportMatrix(matrix_X),
  118.                 MatrixOperation.SubtracteMatrix(vector_Y, vector_PI));
  119.         }
  120.  
  121.         /// <summary>
  122.         /// 计算每次迭代完成后的误差
  123.         /// </summary>
  124.         /// <param name="error">误差</param>
  125.         private void GetIterationError(ref double error)
  126.         {
  127.             for (int i = 0; i < vector_Omega.GetLength(0); i++)
  128.             {
  129.                 error += Math.Abs(vector_Omega[i, 0] - old_VectorOmega[i, 0]);
  130.             }
  131.         }
  132.     }
  133. }

Cholesky分解

  1. using System;
  2. using System.Collections.Generic;
  3. using System.Text;
  4. using Model.Matrix;
  5.  
  6. namespace Model.MatrixDecomposition
  7. {
  8.     internal class Cholesky_LDL_Decomposition
  9.     {
  10.         private double[,] matrix_L = null;
  11.         private double[,] matrix_D = null;
  12.  
  13.         public double[,] MatrixL
  14.         {
  15.             get
  16.             {
  17.                 return this.matrix_L;
  18.             }
  19.         }
  20.  
  21.  
  22.         public double[,] MatrixD
  23.         {
  24.             get
  25.             {
  26.                 return this.matrix_D;
  27.             }
  28.         }
  29.         private int matrix_Dimension = 0;
  30.  
  31.         #region Decomposition Matrix
  32.         /// <summary>
  33.         /// 方程AX=B
  34.         /// 利用Cholesky LDL^T分解法,求方程的解
  35.         /// </summary>
  36.         /// <param name="m_A">系数矩阵A</param>
  37.         /// <param name="v_B">列向量,B</param>
  38.         /// <param name="v_X">求得方程的解</param>
  39.         public void Cholesky(double[,] m_A, double[,] v_B, ref double[,] v_X)
  40.         {
  41.             this.matrix_Dimension = m_A.GetLength(0);
  42.             matrix_L = new double[matrix_Dimension, matrix_Dimension];
  43.             matrix_D = new double[matrix_Dimension, matrix_Dimension];
  44.             //为了计算方便,将L的严格下三角存储在A的对应位置上,
  45.             //而将D的对角元素存储A的对应对角位置上
  46.             //double[,] q = (double[,])m_A.Clone();
  47.             for (int i = 0; i < matrix_Dimension; i++)
  48.             {
  49.                 for (int j = 0; j < i; j++)
  50.                 {
  51.                     m_A[i, i] -= m_A[j, j] * m_A[i, j] * m_A[i, j];
  52.                 }//for
  53.                 for (int k = i + 1; k < matrix_Dimension; k++)
  54.                 {
  55.                     for (int n = 0; n < i; n++)
  56.                     {
  57.                         m_A[k, i] -= m_A[k, n] * m_A[n, n] * m_A[i, n];
  58.                     }//for
  59.                     m_A[k, i] /= m_A[i, i];
  60.                 }//for
  61.             }//for
  62.  
  63.             this.GetLDMatrix(m_A);
  64.             this.SolveEquation(v_B,ref v_X);
  65.             //double[,] resut = MatrixOperation.MultipleMatrix(MatrixOperation.MultipleMatrix(MatrixL, MatrixD), MatrixOperation.TransportMatrix(MatrixL));
  66.         }
  67.  
  68.         /// <summary>
  69.         /// 将L和D矩阵分别赋值
  70.         /// </summary>
  71.         /// <param name="m_A">经过Cholesky分解后的矩阵</param>
  72.         private void GetLDMatrix(double[,] m_A)
  73.         {
  74.             for (int i = 0; i < matrix_Dimension; i++)
  75.             {
  76.                 matrix_L[i, i] = 1;
  77.                 matrix_D[i, i] = m_A[i, i];
  78.                 for (int j = 0; j < i; j++)
  79.                 {
  80.                     matrix_L[i, j] = m_A[i, j];
  81.                 }
  82.             }
  83.         }
  84.         #endregion End Decomposition Matrix
  85.  
  86.         #region Solve Equation
  87.         /// <summary>
  88.         /// 求解LY=B => Y
  89.         /// DL^TX=Y => X
  90.         /// </summary>
  91.         /// <param name="v_B">方程AX=B的列向量B</param>
  92.         /// <param name="v_X">求解结果</param>
  93.         private void SolveEquation(double[,] v_B, ref double[,] v_X)
  94.         {
  95.             v_X =(double[,])v_B.Clone();
  96.             //求解LY=B=>Y
  97.             for (int i = 0; i < matrix_Dimension; i++)
  98.             {
  99.                 for (int j = 0; j < i; j++)
  100.                 {
  101.                     v_X[i,0] -= v_X[j,0] * matrix_L[i, j];
  102.                 }
  103.                 v_X[i,0] /= matrix_L[i, i];
  104.             }
  105.  
  106.             //计算DL^T
  107.             double[,] dMultiLT = MatrixOperation.MultipleMatrix(matrix_D,
  108.                 MatrixOperation.TransportMatrix(matrix_L));
  109.             //求解DL^TX=Y => X
  110.             for (int i = matrix_Dimension - 1; i >= 0; i--)
  111.             {
  112.                 for (int j = i + 1; j < matrix_Dimension; j++)
  113.                 {
  114.                     v_X[i,0] -= v_X[j,0] * dMultiLT[i, j];
  115.                 }
  116.                 v_X[i,0] /= dMultiLT[i, i];
  117.             }
  118.         }//end Method
  119.         #endregion
  120.     }
  121. }

MatrixOperation是一个有关矩阵加、减、乘以及特殊矩阵求逆的一个类。