Este artículo proviene de un motor de traducción automática.

Ejecución de pruebas

Implementación de regresiones logísticas con Newton-Raphson

James McCaffrey

Descargar el ejemplo de código

Regresión logística (LR) es una técnica de aprendizaje automático que puede utilizarse para hacer predicciones sobre datos donde la variable dependiente para predecir toma un valor de 0 ó 1.Ejemplos incluyen la predicción o no un paciente morirá debido a enfermedades del corazón, dentro de un cierto número de años (0 = no mueren, 1 = die), basado en el paciente edad, sexo y nivel de colesterol y predecir si o no un equipo de béisbol va a ganar un juego (0 = perder, 1 = victoria) basándose en factores como equipo de promedio de bateo y jarra limpia ejecución promedio de partida.Regresión logística asume que los datos del problema encaja una ecuación que tiene la forma p = 1.0 / (1.0 + e-z) donde z = b0 + (b1)(x1) + (b2)(x2) +...+ (bn)(xn).Las variables x son los predictores y los valores de b son constantes que deben determinarse.Por ejemplo, supongamos que desea predecir la muerte por enfermedades cardíacas.Las variables predictivas sea x 1 = edad del paciente, x 2 = sexo del paciente (0 = macho, 1 = femenino) y x 3 = nivel de colesterol de paciente.Y supongo que de alguna manera han determinado que b0 =-95.0, b1 = 0.4, b2 = -0,9 y b3 = 11.2.Si hay un paciente varón de 50 año de edad cuyo nivel de colesterol es 6.8, entonces z =-95.0 + (0.4)(50) + (-0.9)(0) + (11.2)(6.8) = 1,16 y así p = 1.0 / (1.0 + exp(-1.16)) = 0.7613.El valor de p vagamente puede interpretarse como una probabilidad, por lo que en este caso había concluir que el paciente tiene una 0.7613 probabilidad de morir en el número especificado de años.

La función 1.0 / (1.0 + e-z) se llama función sigmoide.El dominio de los posibles valores de z es todos los números reales.El resultado de la función es un valor entre 0.0 y 1.0, como se muestra en figura 1.No se puede asumir que los datos que se está trabajando con pueden ser modelados por la función sigmoide, pero muchos conjuntos de datos reales de hecho pueden ser modelados con precisión por la función.

Figura 1 la función sigmoide

Cuando mediante regresión logística, el principal problema es cómo determinar los valores de b (a menudo llamado beta) de la ecuación de LR.En la mayoría de las situaciones, tiene algunos datos históricos con resultados conocidos y utilizar una de varias técnicas para encontrar los valores de beta que mejor se adapte a los datos.Después de que haya determinado los valores de beta, puede utilizarlos para hacer predicciones sobre los nuevos datos.Una de las técnicas más comunes para encontrar los valores de beta para una ecuación de regresión logística se llama iterativamente reweighted mínimos cuadrados (IRLS).IRLS comienza con una estimación de los valores de la beta y luego iterativamente calcula un nuevo y mejor juego de las betas hasta que se cumple alguna condición de parada.Existen varias técnicas que pueden utilizarse para determinar un nuevo y mejor conjunto de valores de beta.Uno de los más comunes se denomina Newton-Raphson (NR), que consiste en encontrar la derivada del cálculo de una función, en este caso la derivada de la función sigmoide.Debido a la estrecha conexión entre IRLS y Newton-Raphson en regresión logística, los dos términos a menudo se usan indistintamente.

Aunque hay muchos recursos disponibles que describen las matemáticas complejas detrás de encontrar regresión logística parámetros beta utilizando el método de Newton-Raphson, existen muy pocos, si alguno, guías de implementación para los programadores.Este artículo explica exactamente cómo logística regresión con obras de Newton-Raphson y cómo implementar una solución utilizando el lenguaje C#.Echa un vistazo a la captura de pantalla en figura 2 para ver donde estoy enfilé.

Figura 2 regresión logística con Newton-Raphson

El programa de demostración comienza por generar dos archivos de datos sintéticos.El primero se llama el fichero de capacitación y consta de 80 líneas de datos de edad, sexo, colesterol y muerte.El fichero de capacitación se utiliza para calcular los valores de la beta de LR.El segundo se llama el archivo de prueba y contiene 20 líneas de datos que se utilizan para evaluar la precisión de la ecuación de LR con los valores de beta calculados a partir de los datos de entrenamiento.El programa de demostración carga los valores de predictor de x de los datos de entrenamiento en una matriz y carga los valores de la variable dependiente de y de los datos en un vector.Observe que la matriz de capacitación de X, a menudo llamada una matriz de diseño, tiene una primera columna adicional que consta de todos los valores de 1.0, y que los valores de la Quiniela han sido convertidos a valores numéricos.A continuación, el programa de demostración establece tres condiciones de parada para el algoritmo IRLS, indicado por las variables maxIterations, jumpFactor y epsilon.El programa de demostración utiliza el algoritmo de Newton-Raphson para estimar la b0, b1, b2 y b3 valores de beta que mejor se adapten a los datos de entrenamiento.La demostración concluye evaluando la ecuación resultante de la LR con los valores de la beta calculada precisa cómo está en los datos de prueba.En este ejemplo, se predijeron correctamente 18 de 20 Y valores (90%).

Este artículo se supone que ha avanzado conocimientos de programación y al menos un conocimiento intermedio de la terminología y las operaciones de la matriz, pero no asume que sabes algo de regresión logística.El código que se produjo la captura de pantalla en figura 2 es demasiado grande para presentar en su totalidad aquí, pero el código fuente completo está disponible en archive.msdn.microsoft.com/mag201208TestRun.Debido a la complejidad del algoritmo IRLS/NR, me centraré principalmente en piezas clave del algoritmo en lugar de en el propio código, así usted podrá modificar el código para satisfacer sus propias necesidades o refactorizar a otro lenguaje de programación si lo desea.

Estructura general del programa

Para mayor simplicidad, todo el código que genera la captura de pantalla en figura 2 está contenida en una sola aplicación de consola de C#.La estructura del programa y el método Main, con algunas declaraciones de WriteLine eliminado, se enumeran en figura 3.

Figura 3 estructura de programa

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

Aunque la regresión logística es un tema complejo, el código en figura 3 no es muy complicado como podría parecer en primer lugar porque la mayoría de los métodos mostrados son rutinas auxiliares relativamente corto. Los dos métodos principales son ComputeBestBeta y ConstructNewBetaVector.

El método MakeRawData genera un archivo de datos de edad-sexo-colesterol-muerte cuasialeatorios. La edad es un valor entero aleatorio entre 35 y 80, sexo es m o f con igual probabilidad y colesterol es un semi-random valor real entre 0,1 y 9.9 que se basa en el valor de edad actual. La variable dependiente de muerte se calcula utilizando una ecuación de regresión logística con valores fijos de beta de b0 =-95.0, b1 = 0.4, b2 = -0,9 y b3 = 11.2. Tan MakeRawData genera datos que definitivamente pueden ser modelados usando LR, en lugar de generar datos puramente aleatorios que no es probable que se siguen un modelo LR.

Calcular un nuevo Vector de Beta

El corazón de la regresión logística con Newton-Raphson es una rutina que calcula una nueva, supuestamente mejor, conjunto de valores de beta del conjunto actual de valores. La matemática es muy profunda, pero afortunadamente el resultado neto no es demasiado complejo. En forma de pseudo-equation, el proceso de actualización está dada por:

b [t] = b [t-1] + inv (X'W [t-1] X) X'(Y-p[t-1])

Aquí está la nueva b [t] ("en el tiempo t," no array indexación) vector de beta. En el lado derecho, b [t-1] es el vector de beta antiguo ("en el tiempo t-1"). La función inv es inversión de la matriz. Mayúsculas x es la matriz de diseño — es decir, los valores de las variables predictivas aumentada con una columna principal de 1.0s. Mayúsculas X' es la transpuesta de la matriz de diseño de X. Mayúscula y es el vector de los valores de la variable dependiente (recordar cada valor será 0 o 1). El representa de [t-1] cantidad p el vector predijo de antiguos valores de probabilidad y (que consistirá en valores entre 0.0 y 1.0). La cantidad w mayúscula es una matriz de pesos llamados, que requiere un poco de explicación.

La ecuación de actualización beta y la matriz w se explican mejor con un ejemplo concreto. Supongamos que por simplicidad que el conjunto de entrenamiento consta únicamente de las primeras cinco líneas de datos se muestra en la figura 1. Por lo tanto la matriz de diseño x sería:

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

El vector de variable dependiente y sería:

0

1

0

1

0

Supongamos que el viejo vector beta de valores actualizarse, b [t-1], es:

1.00

0.01

0.01

0.01

Con estos valores para x y beta, el viejo vector de p, p [t-1], es:

0.8226

0.8428

0.8242

0.8512

0.8085

Observe que si suponemos que los valores de p < 0,5 se interpretan como y = 0 y p valores > = 0,5 se interpretan como y = 1, los viejos valores beta predeciría correctamente sólo dos de cada cinco casos en los datos de entrenamiento.

La matriz de pesos w es una m x matriz m donde m es el número de filas de X. Todos los valores de la matriz w son 0.0 excepto los valores de m que son de la diagonal principal. Cada uno de estos valores es igual al valor de p correspondiente multiplicado por 1-p. Para que este ejemplo, W sería tamaño 5 x 5. Sería la celda superior izquierda [0,0] (0.8226)(1-0.8226) = 0.1459. Sería la celda a [1,1] (0.8428)(1-0.8428) = 0.1325 y así sucesivamente. La cantidad (p)(1-p) representa la derivada del cálculo de la función sigmoide.

En la práctica, la matriz w no se calcula explícitamente porque su tamaño podría ser enorme. Si tuvieras 1.000 filas de datos de entrenamiento, matriz w tendría 1.000.000 células. Observe que la ecuación de actualización beta tiene un término W X [t-1], lo que significa la matriz producto de W [t-1] y X. Porque la mayoría de los valores de W [t-1] son cero, la mayoría de los términos de la multiplicación de matriz también son cero. Esto permite tiempos de [t-1] W X a calcularse directamente desde p [t-1] y X, sin construir explícitamente w el. Varias de las referencias de matemáticas que describen IRLS con el algoritmo NR para LR utilizan el símbolo X ~ (X tilde) para el producto de W [t-1] y X. Vea el método ComputeXtilde en la descarga de código para detalles de implementación.

Método ConstructNewBetaVector acepta como parámetros de entrada el vector viejo de beta, la matriz de diseño de X, el vector de variable dependiente y y el vector de probabilidad viejo. El método calcula y devuelve el vector beta actualizada.

El método se implementa como sigue:

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;

Con la colección de matrix y vector métodos auxiliares que se muestra en la figura 3, calcular un nuevo vector de beta es bastante sencillo. Tenga en cuenta que el método realiza inversión de la matriz. Este es un proceso que puede ir mal en muchos sentidos y es una debilidad importante de Newton-Raphson.

Continuando con el ejemplo, matriz Xt (la transpuesta de X) sería:

 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

Matriz (X ~) se calcula de vector p y la matriz x por método auxiliar ComputeXtilde como:

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

Matriz intermedia B, que representa el producto de Xt y X ~ (que a su vez es XtW[t-1]X) sería:

 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

Matriz intermedia c es la inversa de la matriz b y sería:

 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

Matriz intermedia d es el producto de la matriz c y matriz transpuesta de x y se calcula como:

-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

Vector intermedio YP es la diferencia entre vectores y y p [t-1] y sería:

-0.8

 0.2

-0.8

 0.1

-0.8

Vector intermedio e es el producto de la matriz d y vector YP y sostiene los incrementos para agregarse al vieja vector beta. Vector e sería:

 4.1

-0.4

-2.8

 2.3

El vector de nuevo, final beta se calcula sumando que los valores en el vector intermedio e a los viejos valores de beta y en este ejemplo sería:

 5.1

-0.3

-2.8

 2.4

Con los nuevos valores de la beta, los nuevos valores para el vector de p sería:

0.0240

0.9627

0.0168

0.9192

0.0154

Si estos valores de p son interpretados como Y = 0 cuando p < 0.5, entonces, después de sólo una iteración de Newton-Raphson, los valores de beta predicen correctamente los cinco casos en los datos de prueba.

Saber cuándo parar

La técnica de Newton-Raphson para regresión logística iterativamente mejora los valores del vector beta hasta que se cumple alguna condición de parada. Es sorprendentemente difícil saber cuándo detener la iteración. Método ComputeBestBeta maneja esta tarea. El código se presenta en figura 4.

Figura 4 computación el mejor Vector de Beta

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

Uno de los mayores escollos de la regresión logística es iterar demasiadas veces, resultando en un conjunto de valores de beta que encajan perfectamente los datos del entrenamiento pero no encajan bien a otros conjuntos de datos. Esto se llama over-fitting. Saber cuándo parar el proceso de formación en la regresión logística es arte y ciencia de parte. ComputeBestBeta método comienza por inicializar el vector beta para todos los valores de 0.0, calcular los valores de p asociado y luego calcular el error entre los valores de p y los valores de Y. El código de este artículo utiliza el error cuadrado medio — el promedio de las sumas de las diferencias al cuadrado entre p y Y. Otras posibilidades de error informático incluyen desviación absoluta promedio y error de Cruz-entropía.

El bucle principal de procesamiento en ComputeBestBeta repetidamente calcula un nuevo vector de beta y el término de error correspondiente. La condición de parada principal en el algoritmo es un parámetro maxIterations. Recordar que el algoritmo de Newton-Raphson calcula una matriz inversa, que es susceptible a errores. En este caso ComputeBestBeta devuelve encontró el mejor vector de beta (que, desgraciadamente, podría ser el vector de all-cero inicial). Puede que desee iniciar una excepción aquí en su lugar. Otra alternativa es intentar un escape modificando los valores de beta nuevo mediante el ajuste de la media de los valores anteriores y los nuevos valores.

La siguiente condición de parada, se comprueba si el cambio en los nuevos valores de beta es menor que cierto valor pequeño de epsilon de parámetro, utilizando el método auxiliar NoChange. Esto indica que ha convergido Newton-Raphson y, de hecho, existe una buena posibilidad que su modelo es ahora over-fitted. En lugar de devolver el mejor vector de beta se encuentra en este momento, quizás desee devolver el mejor vector de beta de una iteración anterior. La condición de parada siguiente comprueba si alguno de los nuevos valores de la beta han saltado por una cantidad mayor de algún gran valor dado por el parámetro jumpFactor. Esto indica que el método Newton-Raphson posiblemente ha girado fuera de control y desea una excepción en lugar de reoptimización el mejor vector de beta se encuentra.

Otra condición de parada en ComputeBestBeta comprueba ver si los nuevos valores de beta producen un error mayor que los valores de beta anterior hizo. En este ejemplo, si nuevas betas producen un error mayor cuatro veces en una fila, el algoritmo termina y vuelve a encontrar el mejor vector de beta. Quizás desee parametrizar el valor para el número máximo de aumentos consecutivos en error. Cuando se detecta un aumento en el error, el método intenta salir de la situación cambiando los valores en el vector de beta actual a la media de los valores de la beta actual y los valores de la beta recién calculada.

No hay ningún solo mejor conjunto de condiciones de detención. Cada problema de regresión logística requiere un poco de experimentación para sintonizar el algoritmo de Newton-Raphson.

Ventajas vs. Desventajas

Existen varias alternativas a Newton-Raphson para calcular el mejor conjunto de valores de beta en la regresión logística. Newton-Raphson es una técnica de optimización numérica determinista. También puede emplear técnicas no deterministas (es decir, probabilístico) como optimización de enjambre de partículas, algoritmos evolutivos y optimización de forrajeo bacteriana.

La principal ventaja de Newton-Raphson en comparación con alternativas probabilísticos es que en la mayoría de las situaciones NR es mucho más rápido. Pero Newton-Raphson tiene varias desventajas. Porque NR utiliza la inversión de la matriz, el algoritmo fallará cuando se enfrentan a una matriz singular durante el cálculo. Implementaciones simples de Newton-Raphson requieren todos los datos en memoria, lo que significa que hay un límite para el tamaño de la formación y conjuntos de datos que puede utilizar pruebas. Aunque es raro, algunos conjuntos de datos no pueden converger a un conjunto de valores de beta mejores en todo con NR.

Cuando utilizo la regresión logística, a menudo utilizo un ataque doble. Yo experimente primero con un enfoque de Newton-Raphson. Si esta técnica falla, caigo volver a usar optimización de enjambre de partículas para encontrar el mejor conjunto de valores de beta. Es importante notar que la regresión logística no es mágico, y no todos los datos se ajusta a un modelo de regresión logística. Otras técnicas de aprendizaje máquina para modelar datos con una variable dependiente binaria incluyen redes neuronales, máquinas de vector y análisis discriminante lineal de apoyo.

La explicación del proceso de actualización del vector de beta del algoritmo de Newton-Raphson presentado en este artículo y la descarga de código que lo acompaña le debe aprender a manejar con regresión logística utilizando NR. Regresión logística es un tema fascinante, complejo y puede hacer una valiosa adición a su conjunto de habilidades personales.

Dr.James McCaffrey trabaja para Volt Information Sciences Inc., donde administra la capacitación técnica para ingenieros de software que trabajan en el Microsoft Redmond, Washington, campus. Ha trabajado en varios productos de Microsoft Internet Explorer y MSN Search. McCaffrey es el autor de "Recetas de automatización de prueba .net" (Apress, 2006). Él puede ser contactado en jammc@microsoft.com.

Gracias a los siguientes expertos técnicos por su ayuda en la revisión de este artículo: Dan Liebling y Anne Loomis Thompson