使用 Newton-Raphson 对逻辑回归进行编码
.jpg)
回归 (LR) 是一种机器学习技术，可用于对数据进行预测变量要预测哪里需要值为 0 或 1。例子包括预测是否一个病人会内若干年心脏病死 （0 = 不死，1 = 模），基于病人的年龄、 性别和胆固醇水平，和预测是否一个棒球队会赢一场比赛 （0 = 失去，1 = 赢） 基于团队击球平均和启动投手盈余分析运行平均水平等因素。回归假定问题数据拟合的公式，已形成 p = 1.0 / (1.0 + e-z) 凡 z = b0 + (b1)(x1) + (b2)(x2) +。..+ （bn)(xn)。X 变量是预测和 b 值都必须确定的常量。例如，假设您想要预测死于心脏病。让预测变量 x 1 = 病人的年龄、 x 2 = 病人性别 （0 = 男，1 = 女） 和 x 3 = 病人的胆固醇水平。同时假设您以某种方式确定，b0 =-95.0，b1 = 0.4，b2 =-0.9、 b3 = 11.2。如果有 50 岁的男性病人，其胆固醇水平是 6.8，然后 z =-95.0 + (0.4)(50) + (-0.9)(0) + (11.2)(6.8) = 1.16 和如此 p = 1.0 / (1.0 + exp(-1.16)) = 0.7613。P 值松散可以被解释为一个概率，所以在这种情况下您已经结束病人有 0.7613 死亡概率之内指定的年数。
该函数 1.0 / (1.0 + e-z) 被称为乙状结肠功能。Z 的可能值的域是所有实数。函数的结果是 0.0 和 1.0 之间的值，如中所示图 1。您不能假定您正在使用的数据可以通过乙状结肠功能，建模，但很多真实的数据集可其实准确地模拟功能。
.png)
图 1 的乙状结肠功能
当使用逻辑回归时，面临的首要问题是如何确定的 LR 方程的 b （通常称为 beta） 值。在大多数情况下，您有与已知结果的一些历史数据，并使用几种方法之一来查找最适合数据的测试版的值。您已经确定的 β 值后，可以使用它们来对新数据做出的预测。调用迭代加权的最小二乘法 （IRLS） 寻找的 β 值逻辑回归方程的最常见的技术之一。IRLS 开头的 β 值的估计，然后反复的计算一个新的、 更好地设置的 beta 直到满足停车的某些条件。有几种方法，可用于确定新的、 更好的 β 值设置。其中一个最常见的叫牛顿 (NR)，其中涉及查找函数的微积分导数 — — 在这种情况下乙状结肠函数的导数。IRLS 和牛顿之间在回归中的密切联系，因为这两个术语通常可以互换。
虽然有很多资源可用描述找到回归测试参数使用牛顿的背后复杂的数学，程序员很少的情况下，如果有，实施指南存在。这篇文章解释究竟如何回归与牛顿的作品，以及如何实现一个使用 C# 语言的解决方案。看一看在截图图 2 来看看我驶向哪里。
.png)
图 2 牛顿与回归
演示计划开始通过生成两个综合数据文件。第一次被称为培训文件，由 80 行的年龄、 性别、 胆固醇和死亡数据组成。培训文件用来计算 LR beta 值。第二个被称为测试文件，持有 20 行数据，用来从培训数据计算的 β 值评估 LR 公式的准确性。演示计划将培训数据的 X 预估值加载到一个矩阵，并加载到一个向量的 Y 变量值的数据。请注意 X 培训矩阵，通常被称为设计矩阵，有更多的第一列包含 1.0 的所有值，并预估值都已转换为数字值。下一步，演示计划设置的 IRLS 算法，表示由变量 maxIterations、 电磁辐射及 jumpFactor 三个停车条件。演示计划使用牛顿算法估计 b0、 b1、 b2 和 b3 beta 值，最适合的训练数据。评估测试的数据如何准确的计算的 beta 值生成 LR 方程最后演示。在此示例中，正确预计 18 出 20 Y 值 (90%)。
本文假定您具有先进的编程技巧和至少一个中间知识的矩阵运算和术语，但并不假设你知道什么逻辑回归。所产生的截图中的代码图 2 远太大，目前全部在这里，但完整的源代码是可在 archive.msdn.microsoft.com/mag201208TestRun。由于 IRLS/NR 算法的复杂性，我会重点主要关键部件的算法的而不是对自身的代码，那么你将能修改代码，以满足自己的需要或重构它为另一种编程语言，如果你想。
程序的整体结构
为简单起见，所有代码的都生成中的截图图 2 包含在单个 C# 控制台应用程序中。在列出的程序结构和 Main 方法中的，与某些 WriteLine 语句删除， 图 3。
图 3 程序结构
using System;
using System.IO;
namespace LogisticRegressionNewtonRaphson
{
class LogisticRegressionNRProgram
{
static void Main(string[] args)
{
try
{
string trainFile = "trainFile.txt";
string testFile = "testFile.txt";
MakeRawDataFile(80, 3, trainFile);
MakeRawDataFile(20, 4, testFile);
Console.WriteLine("First 5 lines of training data file are: \n");
DisplayRawData(trainFile, 5);
double[][] xTrainMatrix = LoadRawDataIntoDesignMatrix(trainFile);
double[] yTrainVector = LoadRawDataIntoYVector(trainFile);
double[][] xTestMatrix = LoadRawDataIntoDesignMatrix(testFile);
double[] yTestVector = LoadRawDataIntoYVector(testFile);
int maxIterations = 25;
double epsilon = 0.01;
double jumpFactor = 1000.0;
double[] beta = ComputeBestBeta(xTrainMatrix, yTrainVector,
maxIterations, epsilon, jumpFactor);
Console.WriteLine("Newton-Raphson complete");
Console.WriteLine("The beta vector is: ");
Console.WriteLine(VectorAsString(beta, int.MaxValue, 4, 10));
double acc = PredictiveAccuracy(xTestMatrix, yTestVector, beta);
Console.WriteLine("The predictive accuracy of the model on the test
data is " + acc.ToString("F2") + "%\n");
}
catch (Exception ex)
{
Console.WriteLine("Fatal: " + ex.Message);
Console.ReadLine();
}
} // Main
static void MakeRawDataFile(int numLines, int seed, string fileName)
static void DisplayRawData(string fileName, int numLines)
static double[][] LoadRawDataIntoDesignMatrix(string rawDataFile)
static double[] LoadRawDataIntoYVector(string rawDataFile)
static double PredictiveAccuracy(double[][] xMatrix,
double[] yVector, double[] bVector)
static double[] ComputeBestBeta(double[][] xMatrix, double[] yVector,
int maxNewtonIterations, double epsilon, double jumpFactor)
static double[] ConstructNewBetaVector(double[] oldBetaVector,
double[][] xMatrix,
double[] yVector, double[] oldProbVector)
static double[][] ComputeXtilde(double[] pVector, double[][] xMatrix)
static bool NoChange(double[] oldBvector, double[] newBvector, double epsilon)
static bool OutOfControl(double[] oldBvector, double[] newBvector,
double jumpFactor)
static double[] ConstructProbVector(double[][] xMatrix, double[] bVector)
// Matrix and vector routines:
static double MeanSquaredError(double[] pVector, double[] yVector)
static double[][] MatrixCreate(int rows, int cols)
static double[] VectorCreate(int rows)
static string MatrixAsString(double[][] matrix, int numRows,
int digits, int width)
static double[][] MatrixDuplicate(double[][] matrix)
static double[] VectorAddition(double[] vectorA, double[] vectorB)
static double[] VectorSubtraction(double[] vectorA, double[] vectorB)
static string VectorAsString(double[] vector, int count, int digits, int width)
static double[] VectorDuplicate(double[] vector)
static double[][] MatrixTranspose(double[][] matrix) // assumes
matrix is not null
static double[][] MatrixProduct(double[][] matrixA, double[][] matrixB)
static double[] MatrixVectorProduct(double[][] matrix, double[] vector)
static double[][] MatrixInverse(double[][] matrix)
static double[] HelperSolve(double[][] luMatrix, double[] b) // helper
static double[][] MatrixDecompose(double[][] matrix, out int[] perm,
out int tog)
} // Program
} // ns
虽然回归是一个复杂的主题中的代码图 3 并不是因为它可能会首次出现因为所示的方法中的大多数都是相对较短的 helper 例程的那么复杂。 这两种主要方法是 ComputeBestBeta 和 ConstructNewBetaVector。
MakeRawData 方法将生成一个随机的年龄性别胆固醇死亡数据的文件。 年龄是 35 和 80 之间的随机整数值、 性别平等的概率，是 M 或 F 和胆固醇是随机真实之间的一个值为 0.1 至 9.9 所基于的当前时代价值。 死亡因变量使用具有固定的 β 值的 b0 的回归方程计算 =-95.0，b1 = 0.4，b2 =-0.9、 b3 = 11.2。 所以 MakeRawData 生成数据，肯定使用 LR，而不产生纯粹是随机数据，因此可能不会跟随 LR 模型来建模。
计算新的 beta 版矢量
回归与牛顿的核心是一个计算新的例程，想必更好，设置的 β 值从当前集合中的值。 数学是很深，但幸运的是最终的结果不是太复杂了。 在 pseudo-equation 表中，则获更新过程：
b [t] = b [t-1] + inv (X'W [t 1] X) X'(Y-p[t-1])
B [t] 这里是新的 （"在时间 t，"不数组索引） 测试向量。 在右手边，b [t-1] 是旧 （"在时间 t-1"） 的测试向量。 Inv 功能是矩阵反演。 大写 X 是设计矩阵 — — 预测变量的值，即增加用领先列。 大写 X' 是的 X 设计矩阵转置。 大写 Y 是变量值的矢量 （还记得每个值将为 0 或 1）。 数量 p [t-1] 代表向量的旧预测的概率值 y （这将包括 0.0 和 1.0 之间的值）。 大写字母 W 数量是一个所谓的权重矩阵，其中需要解释的一点。
Beta 更新方程和 W 矩阵与一个具体的例子是最好的解释。 假设为简单起见，训练集仅包含数据中所示的首批五名行图 1。 所以设计矩阵 X 就是：
1.00 48.00 1.00 4.40
1.00 60.00 0.00 7.89
1.00 51.00 0.00 3.48
1.00 66.00 0.00 8.41
1.00 40.00 1.00 3.05
Y 变量向量将会是：
0
1
0
1
0
让我们假设要更新的值的旧 beta 向量，b [t-1] 是：
1.00
0.01
0.01
0.01
这些值与 X 和测试版，旧的 p 向量，p [t-1] 是：
0.8226
0.8428
0.8242
0.8512
0.8085
请注意，如果我们假设 p 值 < 0.5 被解释为 y = 0 和 p 值 > = 0.5 被解释为 y = 1，旧的 β 值会正确预测训练数据中的只有两个五例。
权重矩阵 W 是 m x m 矩阵，其中 m 是 X 的行数。但这些是主要的对角线的 m 值 0.0 W 矩阵中的所有值都。每个这些值等于相应的 p 值乘以 1-p。因此，对于此示例，W 将大小 5 x 5。[0，0] 的左上角单元格会 （0.8226)(1-0.8226) = 0.1459。1，[1] 处的单元格会 (0.8428)(1-0.8428) = 0.1325，等等。数量 （p)(1-p) 代表乙状结肠函数的微积分导数。
在实践中，W 矩阵不计算是显式因为它的大小可能会巨大。如果你有 1000 行的培训数据，矩阵 W 会有 1,000,000 的单元格。请注意任期 [t-1] X，这意味着 W [t-1] 和 X 矩阵产品的试用版更新方程。因为大多数的 W [t-1] 中的值均为零，大部分的矩阵乘法条款也都为零。这使得 W [t-1] 次 X 没有显式构造计算直接从 p [t-1] 和 X，w ·几个 LR NR 算法描述 IRLS 的数学引用使用符号 X ~ (X 颚化符) W [t-1] 和 X 的产品。代码下载 ComputeXtilde 方法实现的详细信息，请参阅。
方法 ConstructNewBetaVector 作为输入参数接受 X 设计矩阵、 Y 变量向量和旧的概率向量的旧 beta 向量。该方法计算，并返回更新的测试向量。
此方法实现如下所示：
double[][] Xt = MatrixTranspose(xMatrix); // X'
double[][] A = ComputeXtilde(oldProbVector, xMatrix); // WX
double[][] B = MatrixProduct(Xt, A); // X'WX
double[][] C = MatrixInverse(B); // inv(X'WX)
if (C == null) // inverse failed!
return null;
double[][] D = MatrixProduct(C, Xt); // inv(X'WX)X'
double[] YP = VectorSubtraction(yVector, oldProbVector); // Y-p
double[] E = MatrixVectorProduct(D, YP); // inv(X'WX)X'(y-p)
double[] result = VectorAddition(oldBetaVector, E);
return result;
矩阵与向量中所示的帮助器方法的集合与图 3，计算新的 beta 版矢量是相当简单。 注意该方法执行矩阵反演。 这是一个过程，可以去错在很多方面是牛顿的一项重大缺陷。
继续该示例，将矩阵 Xt (X 的转置）：
1.00 1.00 1.00 1.00 1.00
48.00 60.00 51.00 66.00 40.00
1.00 0.00 0.00 0.00 1.00
4.40 7.89 3.48 8.41 3.05
矩阵 A (X ~） 将从向量 p 和帮助器方法 ComputeXtilde X 矩阵作为计算：
0.15 7.00 0.15 0.64
0.13 7.95 0.00 1.05
0.14 7.39 0.00 0.50
0.13 8.36 0.00 1.07
0.15 6.19 0.15 0.47
中间矩阵 B，代表该产品的 Xt 和 X ~ （这又是 XtW[t-1]X) 将会：
0.70 36.90 0.30 3.73
36.90 1989.62 13.20 208.46
0.30 13.20 0.30 1.11
3.73 208.46 1.11 23.23
中间矩阵 C 是 B 矩阵的逆矩阵，将会：
602.81 -14.43 -110.41 38.05
-14.43 0.36 2.48 -1.02
-110.41 2.48 26.43 -5.77
38.05 -1.02 -5.77 3.36
中间矩阵 D 是矩阵 C 和 X 转置矩阵的产品，并会计算为：
-33.00 37.01 -0.90 -29.80 31.10
0.77 -0.96 0.30 0.66 -0.72
9.52 -7.32 -4.17 4.54 -2.51
-1.86 3.39 -2.24 -0.98 1.76
中间矢量 YP Y 向量和 p [t-1] 之间的区别是，和将会：
-0.8
0.2
-0.8
0.1
-0.8
中间载体 E D 矩阵的产品和矢量 YP，包含要添加到的旧的 beta 向量的增量。 矢量电子将会是：
4.1
-0.4
-2.8
2.3
新的、 最后的测试向量是由添加的值在中间向量 E 到旧值的测试版，并在此示例中，将会获得的：
5.1
-0.3
-2.8
2.4
Beta 版的新值，p 向量的新值将是：
0.0240
0.9627
0.0168
0.9192
0.0154
如果这些 p 值被解释为 Y = 0 时 p < 牛顿的 0.5，然后后只有一, 次迭代，β 值正确预测所有 5 个案件中的测试数据。
知道何时停止
回归的牛顿技术反复提高 beta 向量的值，直到满足一些停车条件。很难令人惊讶，知道什么时候停止迭代。方法 ComputeBestBeta 处理这项任务。该代码显示在图 4。
计算最佳 Beta 矢量图 4
static double[] ComputeBestBeta(double[][] xMatrix, double[] yVector,
int maxIterations, double epsilon, double jumpFactor)
{
int xRows = xMatrix.Length;
int xCols = xMatrix[0].Length;
if (xRows != yVector.Length)
throw new Exception("xxx (error)");
// Initialize beta values
double[] bVector = new double[xCols];
for (int i = 0; i < xCols; ++i) { bVector[i] = 0.0; }
double[] bestBvector = VectorDuplicate(bVector);
double[] pVector = ConstructProbVector(xMatrix, bVector);
double mse = MeanSquaredError(pVector, yVector);
int timesWorse = 0;
for (int i = 0; i < maxIterations; ++i)
{
double[] newBvector = ConstructNewBetaVector(bVector, xMatrix,
yVector, pVector);
if (newBvector == null)
return bestBvector;
if (NoChange(bVector, newBvector, epsilon) == true)
return bestBvector;
if (OutOfControl(bVector, newBvector, jumpFactor) == true)
return bestBvector;
pVector = ConstructProbVector(xMatrix, newBvector);
double newMSE = MeanSquaredError(pVector, yVector); // Smaller is better
if (newMSE > mse) // new MSE is worse
{
++timesWorse; // Update got-worse counter
if (timesWorse >= 4)
return bestBvector;
bVector = VectorDuplicate(newBvector); // Attempt escape
for (int k = 0; k < bVector.Length; ++k)
bVector[k] = (bVector[k] + newBvector[k]) / 2.0;
mse = newMSE;
}
else // Improved
{
bVector = VectorDuplicate(newBvector);
bestBvector = VectorDuplicate(bVector);
mse = newMSE;
timesWorse = 0;
}
} // End main iteration loop
return bestBvector;
}
回归的最大缺陷之一迭代次数过多，导致一套完全适合的训练数据，但不适合任何其他数据集的 β 值。 这被称为过去式。 部分艺术和一部分科学，知道什么时候停止在回归中的训练过程。 方法 ComputeBestBeta 开始初始化为 0.0 的所有值的 beta 向量、 计算所有相关联的 p 值，然后计算的 p 值和 Y 值的误差。 在本文中的代码使用平均平方的误差 — — 平方和的区别 p Y 款项的平均水平。 计算错误的其他可能性包括平均绝对偏差和交叉熵错误。
在 ComputeBestBeta 中的主要加工循环反复计算新的 beta 版矢量和相应错误术语。 MaxIterations 参数是在算法中的主停止条件。 记得牛顿算法计算矩阵的逆，，容易发生故障。 在此案例 ComputeBestBeta 返回最佳测试向量找到 （这不幸的是，可能是最初的全为零矢量）。 您可能想要在这里而引发异常。 另一种方法是通过修改通过调整他们以前的值的平均新的 beta 值和新值来试图逃走。
下一步的停止条件检查，看看所有新的 beta 值的变化是否小于参数 epsilon，使用 NoChange 的帮助器方法的一些小的值。 这表明牛顿已收敛，事实上，有一个好的机会，您的模型现在 over-fitted。 而不是返回此时发现的最佳测试向量，可能需要从更早的迭代返回最佳测试向量。 下一步的停止条件检查是否任何新的 beta 值已经涨了金额大于给定的参数 jumpFactor 一些大型值。 这表明牛顿可能失控，并且您想要引发异常而不是我找到的最佳测试向量旗。
在 ComputeBestBeta 中的另一个停止条件检查，看看是否新的 beta 值实际上产生更大的错误，比以前的 β 值吗。 在此示例中，如果新 beta 产生更大的错误，四次在行中，该算法终止，并返回最佳测试向量发现。 您可能想要参数的值为错误连续上涨的最大数目。 当检测到错误增加时，该方法将尝试通过当前测试向量中的值更改为当前 beta 版中的值和在测试版中新计算值的平均数摆脱这种情况。
没有单一的最佳套停止条件。 每个逻辑回归问题需要有点实验进行调优的牛顿算法。
优势 vs。 缺点
估计回归中的 β 值的最佳的一组有好几种选择到牛顿。 牛顿是确定性的数值优化技术。 您也可以使用非确定性函数 （指概率） 的技术，如粒子群优化，进化算法和细菌觅食优化。
牛顿相比概率替代方案的主要优点是在大多数情况下 NR 要快得多。 但牛顿也有几个缺点。 由于天然橡胶使用矩阵求逆，算法将失败时在计算过程中遇到了一个奇异的矩阵。 牛顿的简单的实现需要所有的数据在内存中，这意味着有培训的大小限制，并测试您可以使用的数据集。 虽然它很罕见，但有些数据集可能没有会聚到一组最佳 beta 值在所有使用天然橡胶。
当我使用逻辑回归时，我常常采用一项双管齐下的攻击。 我第一次用牛顿方法进行试验。 如果这种技术失败，我回到使用粒子群优化查找最佳值集的 beta 版。 它是重要的是要注意回归并不神奇，并不是所有数据都拟合回归模型。 以二进制文件相关变量的模型数据的其他机器学习技术包括神经网络、 支持向量机和线性鉴别分析。
在这篇文章和附带的代码下载中的牛顿算法的 beta 矢量更新过程的解释应该跟你启动并运行回归使用天然橡胶。 逻辑回归是一个引人入胜的、 复杂的主题，可以使你个人的技能集宝贵的补充。
**博士。**詹姆斯 · 麦卡弗里为伏信息科学 Inc.，凡他管理软件工程师工作在微软雷德蒙德，华盛顿，校园的技术培训工作。 他已为多种 Microsoft 产品效过力，包括 Internet Explorer 和 MSN Search。 McCaffrey 是《.NET 软件测试自动化之道》（Apress，2006）的作者。 他可以在达到 jammc@microsoft.com。
衷心感谢以下技术专家对本文的审阅：丹利布灵和安妮 · 梅斯 · 汤普森