生活随笔
收集整理的這篇文章主要介紹了
逻辑斯谛回归模型
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
????邏輯斯諦回歸模型是研究因變量為二分類或多分類觀察結果與影響因素之間的關系的一種概率型非線性回歸模型。邏輯斯諦回歸系數通過最大似然估計得到。Logistic函數如下：
???? ????
????式中x為
???? ????
????這里是輸入變量的n個特征，然后按照Logistic函數形式求出。
????假設有n個獨立變量的向量，設條件概率在x條件下y發生的概率（假設y=1為y發生）。則Logistic函數表示為：
???? ????
????同理，在x條件下y不發生的概率為：
????Logistic回歸都是圍繞Logistic函數展開的，如何求解是Logistic回歸模型的主要問題，采用的最大似然估計來求解這組參數。
????假設有m個觀測樣本，觀測值分別為，設為給定條件下的概率，同理的概率為，得到一個觀測值得概率為：
???? ????
????因為各觀測樣本間相互獨立，于是得到似然函數：
????????
????對似然函數取對數：
????????
????現要求向量使的最大，其中：
????????
????要求的最大似然估計，我們需要確定似然函數存在局部極大值。因此，對似然函數求偏導后得：
???? ????
????由多元函數極值的必要條件可知，若多元函數在一點取得極值，且一階偏導存在，則該點處所有一階偏導為0。由此，可以得出n+1個方程，如下：
???? ????
????由此方程解出的 不一定是似然函數的極值，需要通過Hessian矩陣來判斷得出的解是否為似然函數的極值。
????Hessian矩陣是一個多元函數的二階偏導構成的方陣，描述了函數的局部曲率。對一個多元函數 ，如果他的二階偏導都存在，那么Hessian矩陣如下：
???? ????
????通過Hessian矩陣，我們可以判斷一點M處極值的三種情況：
如果 是正定矩陣，則臨界點M處是一個局部極小值； 如果 是負定矩陣，則臨界點M處是一個局部極大值； 如果 是不定矩陣，則臨界點M處不是極值。 對于中的n+1個方程，要求Hessian矩陣，先要求似然函數的二階偏導，即：
???? ????
????則似然函數的Hessian矩陣為
???? ????
????設有矩陣X、A：
???? ????
???? ????
????則似然函數的Hessian矩陣可表示為：
???? ????
????顯然，矩陣A是負定的，則可以證明H也是負定的，說明似然函數存在局部極大值。因此，可以使用牛頓迭代法（Newton's Method）來求 。
????對一元函數，使用牛頓迭代法來求零點。假設要求 的解，首先選取一個點 作為迭代起點，通過下面的式子進行迭代，直到達到指定精度為止。
???? ????
????由此，有時起始點選擇很關鍵，如果函數只存在一個零點，那么這個起始點選取就無關重要。對已Logistic回歸問題，Hessian矩陣對于任意數據都是負定的，所以說極值點只有一個，初始點的選取無關緊要。
????因此，對于上述Logistic回歸的似然函數， 令：
???? ????
????則由可以得到如下的迭代式子：
???? ????
????由于Hessian矩陣是負定的，將矩陣A提取一個負號，得：
???? ????
????然后Hessian矩陣變為
???? ????
????這樣，Hessian矩陣就是對稱正定的了。那么牛頓迭代式變為：
???? ????
????現在，關鍵是如何快速并有效的計算，即解方程組 。由于 是對稱正定的，可以使用Cholesky矩陣分解法來解。
????若 對稱正定，則存在一個對角元為正數的下三角矩陣，使得 成立。對于，可以通過以下步驟求解：
求 的Cholesky分解，得到 求解 ，得到 求解 ，得到 現在的關鍵問題是對 進行Cholesky分解。假設：
???? ????
????通過比較兩邊的關系，首先由
???? ????
再由
???? ????
????這樣，得到了矩陣 的第一列元素。假設，已經算出了 的前 列元素，通過
???? ????
????可以得出
???? ????
????進一步由
???? ????
????最終：
???? ????
????這樣便通過 的前 列求出了第 列，一直遞推下去即可求出 。這種方法稱為平方根法。
????利用上述方法需要進行開方，這有可能損失精度和增加運算量，為了避免開方，將Cholesky分解進行改進，即：
???? ????
????其中： 是單位下三角矩陣， 為對角均為正數的對角矩陣。把這一分解叫分解。設：
???? ????
????則對于，求解步驟變為：
求 的分解，得到 求解 ，得到 求解 ，得到 對比兩邊元素，可以得到：
???? ????
????由此可以確定 和 的公式如下：
???? ????
???? ????
???? ????
牛頓迭代法
using System; using System.Collections.Generic; using System.Text; using Model.MatrixDecomposition; using Model.Matrix; using System.IO; namespace Model.NewtonRaphson { ????internal class NewtonRaphsonIterate ????{ ????????Cholesky_LDL_Decomposition decompositionMatrixByLDL = new Cholesky_LDL_Decomposition(); ????????private double[,] matrix_Hessian = null; ????????private double[,] matrix_X = null; ????????private double[,] matrix_A = null; ????????private double[,] vector_HU = null; ????????private double[,] vector_U = null; ????????private double[,] vector_Y = null; ????????private double[,] vector_Omega = null; ????????private double[,] vector_PI = null; ????????private double[,] old_VectorOmega = null; ????????#region 屬性 ????????public double[,] MatrixL ????????{ ????????????get ????????????{ ????????????????return decompositionMatrixByLDL.MatrixL; ????????????} ????????} ????????public double[,] MatrixD ????????{ ????????????get ????????????{ ????????????????return decompositionMatrixByLDL.MatrixD; ????????????} ????????} ????????public double[,] Hessian ????????{ ????????????get ????????????{ ????????????????return this.matrix_Hessian; ????????????} ????????} ????????#endregion ????????/// <summary> ????????/// 執行牛頓迭代法 ????????/// </summary> ????????/// <param name="matrix_X">輸入矩陣</param> ????????/// <param name="vector_Y">輸出向量</param> ????????/// <param name="error_Threashold">迭代停止閾值</param> ????????/// <param name="omega">計算完成的參數</param> ????????internal void Run(double[,] matrix_X, double[,] vector_Y, double error_Threashold, ref double[,] omega) ????????{ ????????????double error = 1; ????????????this.matrix_X = matrix_X; ????????????this.vector_Y = vector_Y; ????????????this.vector_Omega = omega; ????????????InitializeParameter(matrix_X.GetLength(0)); ????????????int i = 0; ????????????while (error > error_Threashold) ????????????{ ????????????????i++; ????????????????error = 0; ????????????????old_VectorOmega = (double[,])vector_Omega.Clone(); ????????????????GetMatrixAAndVectorPI(); ????????????????matrix_Hessian = MatrixOperation.MultipleMatrix( ????????????????????MatrixOperation.MatrixMultiDiagMatrix( ????????????????????MatrixOperation.TransportMatrix(matrix_X), matrix_A), ????????????????????matrix_X); ????????????????GetMatrixU(); ????????????????decompositionMatrixByLDL.Cholesky((double[,])matrix_Hessian.Clone(), vector_U, ref vector_HU); ????????????????vector_Omega = MatrixOperation.AddMatrix(vector_Omega, vector_HU); ????????????????GetIterationError(ref error); ????????????????//TODO:迭代計算 ????????????} ????????????omega = (double[,])vector_Omega.Clone(); ????????} ????????private void InitializeParameter(int rowNumber) ????????{ ????????????matrix_A = new double[rowNumber, 1]; ????????????vector_PI = new double[rowNumber, 1]; ????????} ????????/// <summary> ????????/// 獲取H=X^TAX的A以及PI(Xi) ????????/// </summary> ????????private void GetMatrixAAndVectorPI() ????????{ ????????????for (int i = 0; i < matrix_X.GetLength(0); i++) ????????????{ ????????????????double temp_A = 0; ????????????????//matrix_A[i, 0] += vector_Omega[0, 0]; ????????????????for (int j = 0; j < matrix_X.GetLength(1); j++) ????????????????{ ????????????????????temp_A += matrix_X[i, j] * vector_Omega[j, 0]; ????????????????}//for2 ????????????????matrix_A[i, 0] = (1.0) / (1 + Math.Exp(-temp_A)); ????????????????vector_PI[i, 0] = matrix_A[i, 0]; ????????????????matrix_A[i, 0] *= (1 - matrix_A[i, 0]); ????????????}//for1 ????????} ????????//計算HU中的U ????????//U=X^T(Y-PI) ????????private void GetMatrixU() ????????{ ????????????vector_U = MatrixOperation.MultipleMatrix(MatrixOperation.TransportMatrix(matrix_X), ????????????????MatrixOperation.SubtracteMatrix(vector_Y, vector_PI)); ????????} ????????/// <summary> ????????/// 計算每次迭代完成后的誤差 ????????/// </summary> ????????/// <param name="error">誤差</param> ????????private void GetIterationError(ref double error) ????????{ ????????????for (int i = 0; i < vector_Omega.GetLength(0); i++) ????????????{ ????????????????error += Math.Abs(vector_Omega[i, 0] - old_VectorOmega[i, 0]); ????????????} ????????} ????} } Cholesky分解
using System; using System.Collections.Generic; using System.Text; using Model.Matrix; namespace Model.MatrixDecomposition { ????internal class Cholesky_LDL_Decomposition ????{ ????????private double[,] matrix_L = null; ????????private double[,] matrix_D = null; ????????public double[,] MatrixL ????????{ ????????????get ????????????{ ????????????????return this.matrix_L; ????????????} ????????} ????????public double[,] MatrixD ????????{ ????????????get ????????????{ ????????????????return this.matrix_D; ????????????} ????????} ????????private int matrix_Dimension = 0; ????????#region Decomposition Matrix ????????/// <summary> ????????/// 方程AX=B ????????/// 利用Cholesky LDL^T分解法，求方程的解 ????????/// </summary> ????????/// <param name="m_A">系數矩陣A</param> ????????/// <param name="v_B">列向量，B</param> ????????/// <param name="v_X">求得方程的解</param> ????????public void Cholesky(double[,] m_A, double[,] v_B, ref double[,] v_X) ????????{ ????????????this.matrix_Dimension = m_A.GetLength(0); ????????????matrix_L = new double[matrix_Dimension, matrix_Dimension]; ????????????matrix_D = new double[matrix_Dimension, matrix_Dimension]; ????????????//為了計算方便，將L的嚴格下三角存儲在A的對應位置上， ????????????//而將D的對角元素存儲A的對應對角位置上 ????????????//double[,] q = (double[,])m_A.Clone(); ????????????for (int i = 0; i < matrix_Dimension; i++) ????????????{ ????????????????for (int j = 0; j < i; j++) ????????????????{ ????????????????????m_A[i, i] -= m_A[j, j] * m_A[i, j] * m_A[i, j]; ????????????????}//for ????????????????for (int k = i + 1; k < matrix_Dimension; k++) ????????????????{ ????????????????????for (int n = 0; n < i; n++) ????????????????????{ ????????????????????????m_A[k, i] -= m_A[k, n] * m_A[n, n] * m_A[i, n]; ????????????????????}//for ????????????????????m_A[k, i] /= m_A[i, i]; ????????????????}//for ????????????}//for ????????????this.GetLDMatrix(m_A); ????????????this.SolveEquation(v_B,ref v_X); ????????????//double[,] resut = MatrixOperation.MultipleMatrix(MatrixOperation.MultipleMatrix(MatrixL, MatrixD), MatrixOperation.TransportMatrix(MatrixL)); ????????} ????????/// <summary> ????????/// 將L和D矩陣分別賦值 ????????/// </summary> ????????/// <param name="m_A">經過Cholesky分解后的矩陣</param> ????????private void GetLDMatrix(double[,] m_A) ????????{ ????????????for (int i = 0; i < matrix_Dimension; i++) ????????????{ ????????????????matrix_L[i, i] = 1; ????????????????matrix_D[i, i] = m_A[i, i]; ????????????????for (int j = 0; j < i; j++) ????????????????{ ????????????????????matrix_L[i, j] = m_A[i, j]; ????????????????} ????????????} ????????} ????????#endregion End Decomposition Matrix ????????#region Solve Equation ????????/// <summary> ????????/// 求解LY=B => Y ????????/// DL^TX=Y => X ????????/// </summary> ????????/// <param name="v_B">方程AX=B的列向量B</param> ????????/// <param name="v_X">求解結果</param> ????????private void SolveEquation(double[,] v_B, ref double[,] v_X) ????????{ ????????????v_X =(double[,])v_B.Clone(); ????????????//求解LY=B=>Y ????????????for (int i = 0; i < matrix_Dimension; i++) ????????????{ ????????????????for (int j = 0; j < i; j++) ????????????????{ ????????????????????v_X[i,0] -= v_X[j,0] * matrix_L[i, j]; ????????????????} ????????????????v_X[i,0] /= matrix_L[i, i]; ????????????} ????????????//計算DL^T ????????????double[,] dMultiLT = MatrixOperation.MultipleMatrix(matrix_D, ????????????????MatrixOperation.TransportMatrix(matrix_L)); ????????????//求解DL^TX=Y => X ????????????for (int i = matrix_Dimension - 1; i >= 0; i--) ????????????{ ????????????????for (int j = i + 1; j < matrix_Dimension; j++) ????????????????{ ????????????????????v_X[i,0] -= v_X[j,0] * dMultiLT[i, j]; ????????????????} ????????????????v_X[i,0] /= dMultiLT[i, i]; ????????????} ????????}//end Method ????????#endregion ????} } MatrixOperation是一個有關矩陣加、減、乘以及特殊矩陣求逆的一個類。
轉載于:https://www.cnblogs.com/reddatepz/p/4496362.html
總結
以上是生活随笔為你收集整理的逻辑斯谛回归模型的全部內容，希望文章能夠幫你解決所遇到的問題。
如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。
