次の方法で共有


この記事は機械翻訳されたものです。

テストの実行

ニュートン ラフソン法を使用したロジスティック回帰のコードを作成する (機械翻訳)

James McCaffrey

コード サンプルのダウンロード

James McCaffreyロジスティック回帰 (LR) 上のデータを予測する予測される従属変数が 0 または 1 の値をかかります使用できます機械学習手法です。 例予測患者心臓病による特定数年内で死ぬかどうか (0 死ぬことはない 1 = ダイ =) は、患者の年齢、性別、コレステ ロール レベルおよび野球チームは試合に勝つでしょうかどうかの予測に基づいて (0 = 1 を失う = win) 平均打撃成績・投手平均開始チームなどの要因に基づいて。 ロジスティック回帰を想定しています問題データ フォーム 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。 コレステロールのレベルがある 6.8 は、z 50 歳の男性患者の場合-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。 シグモイド関数では、作業しているデータをモデル化できますが、多くの現実のデータ セット実際関数によって正確にモデリングすることができると仮定することはできません。

The Sigmoid Function
図 1 シグモイド関数

ロジスティック回帰を使用する場合、主な問題は LR 方程式の b ("ベータ"とも呼ばれます) の値を確認する方法です。 ほとんどの状況でいくつかの履歴データと既知の結果があるし、データに最適なベータの値を検索するいくつかのテクニックのいずれかを使用します。 ベータ版の値を決定した後は、新しいデータで予測を行うできます。 ロジスティック回帰方程式のベータ値を見つけるための最も一般的な手法の 1 つは、繰り返し reweighted 最小二乗 (IRLS) と呼ばれます。 IRLS は、ベータ値の推定値で開始し、繰り返しを計算、新しい一部停止条件が満たされるまでより良いのベータ版を設定します。 新しいを決定するより良いベータの値使用できますいくつかの方法があります。 1 つの最も一般的なニュートン ・ ラプソン法 (NR) は、関数の微積分の誘導体を見つけること含むと呼ばれる — この場合はシグモイド関数の微分。 ロジスティック回帰のため閉じる接続 IRLS とニュートン ・ ラプソン法の間、2 つの用語は同義です。

ロジスティック回帰ニュートン ・ ラプソン法を用いたベータ パラメーターを見つけることの背後にある複雑な数学を記述する利用可能なリソースの多くがありますが、プログラマのための非常に少ない場合は、実装ガイドが存在します。 正確にどのようにロジスティック回帰ニュートン-ラフソン作品と c# 言語を使用してソリューションを実装する方法について説明します。 スクリーン ショットを見て図 2 にどこに向かっているを参照してください。

Logistic Regression with Newton-Raphson
図 2 ニュートン ・ ラプソン法とロジスティック回帰

デモ プログラムは、2 つの合成データ ファイルを生成することによって開始します。 最初は、トレーニング ファイルと呼ばれ、年齢、性別、コレステロールおよび死データの 80 の行で構成されています。 トレーニング ファイルは、LR のベータ値を計算するために使用されます。 2 番目は、テスト ファイルと呼ばれ、20 行の訓練データから計算ベータ値を持つ LR 方程式の精度を評価するために使用するデータを保持しています。 デモ プログラムは、行列のトレーニング データの X 予測値の読み込み、ベクトルにデータの Y 従属変数の値を読み込みます。 デザイン マトリックスと呼ばれる X トレーニング マトリックス、すべての 1.0 の値で構成される追加の最初の列が、その予測値はすべて数値に変換されています。 次に、デモ プログラムは、3 つの停止条件変数 maxIterations, イプシロンと jumpFactor によって示される、IRLS のアルゴリズムを設定します。 デモ プログラムは、b0、b1、トレーニング データに最適な b2 と b3 のベータ値を推定するニュートン ・ ラプソン法のアルゴリズムを使用します。 どのように正確な結果の LR 方程式計算ベータ値を持つテスト データには評価することによって、デモを終了します。 この例では、18 のうち 20 Y 値 (90%) が適正に予測されました。

この資料では、プログラミングのスキルとは、少なくとも中間の知識の行列演算や専門用語、高度がロジスティック回帰について何を知っているものではありませんを想定しています。 スクリーン ショットを生成するコード図 2 ここでは、完全に提示するには大きすぎるが、完全なソース コードでは archive.msdn.microsoft.com/mag201208TestRun。 IRLS/NR のアルゴリズムの複雑さのため、コードにだからあなた独自のニーズを満たす、またはご希望の場合は別のプログラミング言語にリファクタリングするコードを変更することができますではなく、アルゴリズムのキーの部分に主を取り上げます。

全体的なプログラムの構造

簡略化のためのすべてのコードが生成のスクリーン ショットで図 2 1 つの c# コンソール アプリケーションに含まれます。 プログラムの構造と削除、いくつかの WriteLine ステートメントの Main メソッドは表示されます図 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 のメソッドのほとんどは比較的短いヘルパー ルーチンですので最初に見えるかもしれませんかなり複雑ではありません。 2 つの主要メソッドは、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 を使用してモデル化できるデータを生成します。

新しいベータ ベクトルを計算

ロジスティック回帰ニュートン ・ ラプソン法での中心を計算する新しいルーチンは、おそらく良く、ベータの値の現在のセットから値します。 数学は非常に深いですが、幸いにも、ネットの結果はあまり複雑でないです。 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 は設計 — つまり、予測変数の値が 1.0年の先頭の列を拡張現実感します。 大文字の X' X デザイン行列の転置行列です。 大文字 Y がベクトルの従属変数の値 (各値は、0 または 1 になりますを思い出してください)。 数量 p [t 1] を表すベクトルの古い確率 (0.0 と 1.0 の間の値で構成されます) Y の予測値。 大文字 W 数量は説明のビットを必要とするいわゆる重み行列です。

ベータ版更新方程式と W の行列は具体的な例と最もよく説明されます。 表示されるデータの最初の 5 行のみトレーニングのセットで構成されますシンプルさのためと仮定図 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

古いベータ ベクトルの値を更新するは、b [t-1] であることと仮定しましょう。

1.00
0.01
0.01
0.01

これらの値を X と p [t-1] 古い p ベクトル ベータ版です:

0.8226
0.8428
0.8242
0.8512
0.8085

私たちの p 値を推定することに気づく < 0.5 y として解釈される 0 と p 値 = > = 0.5 y として解釈されます = 1、古いベータ値が正しくトレーニング データ内の唯一の 2 のうち 5 つのケースを予測します。

ウェイト マトリックス W x m マトリックスの行 X の数は m です。 W のマトリックス内のすべての値は、メインの対角線上にあるそれらの m 値を除き 0.0 です。 これらの各値に対応するを 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 行列計算されていません。 1,000 行のトレーニング データがあった場合は、マトリックス W 1,000,000 セル必要があります。 注意ベータ版更新方程式 W [t 1] と X のマトリックス製品を意味する用語 W [t 1] X を持っています。 ほとんどの W [t-1] の値がゼロのため、行列乗算の用語のほとんどはまたゼロです。 これは、W [t 1] 回できます直接 p [t 1] から、X を明示的に構築することがなく計算する X w. いくつかの LR の NR のアルゴリズムに IRLS を記述する数学記号 X を使用する ~ (チルダ) X W [t 1] と X の製品。 メソッド ComputeXtilde は、コード ダウンロードの実装の詳細を参照してください。

メソッド ConstructNewBetaVector は古いベータ ベクトルの X デザイン行列、Y 従属変数のベクトル、古い確率ベクトル入力パラメーターとして受け入れます。 メソッドを計算し、更新されたベータのベクターを返します。

メソッドを実装しましょう。

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、新しいベータ ベクトル計算の非常に簡単です。 メソッドがマトリックスの反転を実行することに注意してください。 これは多くの方法で間違って行くことができ重要な弱さニュートン ・ ラプソン法のプロセスです。

例に続けて、行列 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 ベクトルおよび行列 X ヘルパー メソッド ComputeXtilde として計算されます。

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 は、古いベータ版ベクターに追加する増分値を保持しています。 ベクトル E はようになります。

4.1
-0.4
-2.8
 2.3

新しい、最終的なベータ ベクトル中間ベクトル E のベータ版は、古い値をこの例の値になるを追加することによって得られます。

5.1
-0.3
-2.8
 2.4

ベータ版の新しい値では、p ベクトルの新しい値になります。

0.0240
0.9627
0.0168
0.9192
0.0154

これらの p 値は、Y として解釈される場合 = 0 と p < 0.5、そして後にのみ 1 つのイテレーション ニュートン ・ ラプソン法のベータ値は正しくテスト データ内のすべての 5 例を予測します。

停止するときに知ること

ロジスティック回帰は、ニュートン-ラフソン法は繰り返し一部停止条件が満たされるまでベータ ベクトルの値を向上させます。 イテレーションを停止するときに知っているは驚くほど難しい。 ComputeBestBeta メソッドはこのタスクを処理します。 そのコードについては、図 4を参照してください。

図 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;
}

トレーニング データを完璧にフィットが、他のデータ セットに適合していないベータ値のセットを生成ロジスティック回帰の最大の落とし穴の 1 つ何度も繰り返し処理です。 これは仮説と呼ばれます。 ロジスティック回帰の訓練過程を停止するときに知ることは、部分芸術と部分科学です。 メソッド ComputeBestBeta ベータ ベクトルをすべて 0.0 の値を初期化、関連付けられている p 値の計算、p 値と Y の値との誤差を計算によって開始されます。 この資料のコードは平均 2乗誤差を使用します — p と Y の間の平方差の合計の平均。 計算誤差のための他の可能性は、平均絶対偏差とクロス エントロピー エラーがあります。

ComputeBestBeta のメインの処理ループには繰り返し新しいベータ ベクトルと対応するエラー用語を計算します。 アルゴリズムの一次停止状態は、maxIterations パラメーターです。 ニュートン ・ ラプソン法アルゴリズム失敗に影響を受けやすい、行列の逆行列を計算することを思い出してください。 (これは、残念ながら、初期の全零ベクトルかもしれない) 最高のベータ ベクトルを発見このケース ComputeBestBeta リターンで。 代わりにここでは例外をスローする可能性があります。 別の方法としては、それらの以前の値の平均を調整することによって、新しいベータ値と新しい値を変更することによってエスケープを試みることです。

次の停止条件をすべての新しいベータの値の変更がいくつか小さな NoChange のヘルパー メソッドを使用して、パラメーター epsilon の値より小さいかどうかをチェックします。 これは、ニュートン ・ ラプソン法の収束し、実際には、あなたのモデルを今すぐ over-fitted 良いチャンスがあることを示します。 この時点で発見した最高のベータ ベクトルを返す代わりに、以前のイテレーションから最高のベータ ベクトルを返すしたいと思うかもしれない。 次の停止条件の新しいベータ値パラメーター jumpFactor によって与えられたいくつか大きい値よりも大きい量急増しているかどうかを確認します。 これは、ニュートン-ラフソンおそらくスピンアウトしています見つけた最高のベータ版のベクトルを終えてのではなく、例外をスローすることを示します。

新しいベータ値が実際には以前のベータ値より大きいエラーを生成するかどうかを見ていた ComputeBestBeta で別の停止条件をチェックします。 この例より大きいエラー 4 回連続で新しいベータが生成する場合アルゴリズム終了し、最高のベータ ベクトルを検索して返します。 エラーの連続増加の最大数の値をパラメーター化します。 エラーの増加が検出されると、メソッドは現在のベータ版の値と新しく計算されたベータの値の平均を現在ベータ ベクトルの値を変更することによって、状況を脱出しようとします。

停止条件の 1 つの最高のセットはありません。 すべてのロジスティック回帰問題ニュートン ・ ラプソン法アルゴリズムをチューニングする実験のビットが必要です。

利点対。 短所

最高のロジスティック回帰のベータ値セットを推定するためのニュートン ・ ラプソン法にいくつかの選択肢があります。 ニュートン-ラフソン決定的数値最適化手法です。 粒子群最適化、進化的アルゴリズムと細菌の採餌の最適化などの非決定的の (確率的意味) の手法を使用することもできます。

確率的に比べてニュートン ・ ラプソン法の主な利点は、ほとんどの状況で NR がはるかに高速であることです。 しかしニュートン ・ ラプソン法はいくつかの欠点があります。 NR は逆行列を使用するため、アルゴリズムは特異行列計算中に出会うときに失敗します。 ニュートン ・ ラプソン法の単純な実装は、訓練のサイズに制限は意味する、メモリ内にすべてのデータを必要し、テスト データ セットを使用できます。 まれですが、いくつかのデータセットはすべてで NR を使用して最高のベータ値のセットに収束可能性があります。

ロジスティック回帰を使用するは、2 つの方面からの攻撃を使用します。 私は最初、ニュートン-ラフソン アプローチで実験します。 その技術が失敗した場合、私は戻って粒子群最適化を使用して最高のベータ値セットを見つけるために落ちる。 ロジスティック回帰は、魔法ではない、すべてのデータに適合ロジスティック回帰モデルに注意する重要です。 その他機械学習手法のモデル データをバイナリに依存する変数には、サポート ベクトル マシンと線形判別分析ニューラル ネットワークなどがあります。

ベータ ベクトルの更新プロセスのこの記事と付属のコード ダウンロードに提示ニュートン ・ ラプソン法アルゴリズムの説明 NR を使用してロジスティック回帰を起動および実行を取得する必要があります。 ロジスティック回帰は、魅力的な複雑なトピックです、あなたの個人的な技術セットに貴重な追加を行うことができます。

Dr.James McCaffrey ボルト情報科学株式会社は、ここで彼はマイクロソフトのレドモンド、ワシントン州でキャンパスで働くソフトウェア エンジニアの技術トレーニングを管理するための作品します。 これまでに、Internet Explorer、MSN サーチなどの複数のマイクロソフト製品にも携わってきました。 McCaffrey は、「.NET テスト オートメーション レシピ」(Apress、2006年) の著者です。 彼は到達することができます jammc@microsoft.com

この記事のレビュー、次技術専門家のおかげで。ダン Liebling、アン ・ ルーミス ・ トンプソン