Este artigo foi traduzido por máquina.

Execução de teste

Regressão logística de codificação com Newton-Raphson

James McCaffrey

Baixar o código de exemplo

Regressão logística (LR) é uma técnica de aprendizado de máquina que pode ser usada para fazer previsões sobre dados onde a variável dependente a ser previsto tem um valor de 0 ou 1.Exemplos incluem prevendo ou não um paciente vai morrer devido à doença cardíaca dentro de um certo número de anos (0 = não morrer, 1 = die), com base do paciente idade, sexo e nível de colesterol e prevendo ou não uma equipe de beisebol vai ganhar um jogo (0 = perder, 1 = win) baseada em fatores como equipe média de batedura e começando a média de execução de acumulado do arremessador.Regressão logística assume que os dados do problema se encaixa uma equação que tem a forma p = 1.0 / (1.0 + e-z) onde z = b0 + (b1)(x1) + (b2)(x2) +...+ (bn)(xn).As variáveis de x são os preditores e os valores de b são constantes que devem ser determinados.Por exemplo, suponha que você deseja prever morte por doença cardíaca.As variáveis do predictor seja x 1 = a idade dos pacientes, x 2 = paciente sexo (0 = macho, 1 = feminino) e x 3 = nível de colesterol do paciente.E suponha que de alguma forma ter determinado que b0 =-95.0, b1 = 0,4, b2 = -0,9 e b3 = 11,2.Se houver um paciente do sexo masculino 50 ano de idade cujo nível de colesterol é 6,8, então z =-95.0 + (0.4)(50) + (-0.9)(0) + (11.2)(6.8) = 1,16 e então p = 1.0 / (1.0 + exp(-1.16)) = 0.7613.O valor de p pode frouxamente ser interpretada como uma probabilidade, então neste caso você tivesse concluir que o paciente tem uma probabilidade de 0.7613 de morrer dentro do número especificado de anos.

A função 1.0 / (1.0 + e-z) é chamado a função sigmóide.O domínio de valores possíveis para z é todos os números reais.O resultado da função é um valor entre 0,0 e 1,0, como mostrado na Figura 1.Você não pode presumir que os dados que você está trabalhando com podem ser modelados pela função sigmóide, mas muitos conjuntos de dados reais de fato podem ser modelados com precisão pela função.


Figura 1 A função sigmóide

Quando usando regressão logística, o principal problema é como determinar os valores de b (muitas vezes chamado de beta) para a equação de LR.Na maioria das situações, você tem alguns dados históricos com resultados conhecidos e usa uma das várias técnicas para encontrar os valores do beta que melhor se ajustam aos dados.Depois de determinar os valores de beta, você pode usá-los para fazer previsões sobre novos dados.Uma das técnicas mais comuns para encontrar os valores de beta para uma equação de regressão logística é chamada iterativamente reweighted least squares (IRLS).IRLS começa com uma estimativa dos valores de beta e depois iterativamente calcula um novo, melhor conjunto de betas, até que alguma condição de parada.Existem várias técnicas que podem ser usadas para determinar uma nova, melhor conjunto de valores de beta.Um dos mais comuns é chamado de Newton-Raphson (NR), que envolve encontrar a derivada de cálculo de uma função — neste caso, a derivada da função sigmóide.Por causa da estreita conexão entre IRLS e Newton-Raphson na regressão logística, os dois termos são muitas vezes intercambiáveis.

Embora há uma abundância de recursos disponíveis que descrevem a matemática complexa por trás de encontrar regressão logística parâmetros beta usando Newton-Raphson, existem muito poucos, se houver, guias de implementação para programadores.Este artigo explica a regressão logística exatamente como com obras de Newton-Raphson e como implementar uma solução usando a linguagem c#.Dê uma olhada no screenshot no Figura 2 para ver onde quero chegar.


Figura 2 a regressão logística com Newton-Raphson

O programa de demonstração começa por gerar dois arquivos de dados sintéticos.O primeiro é chamado de arquivo de treinamento e consiste de 80 linhas de dados de idade, sexo, colesterol e morte.O arquivo de formação é usado para calcular os valores de beta de LR.O segundo é chamado de arquivo de teste e possui 20 linhas de dados que são usados para avaliar a precisão da equação de LR com os valores de beta, calculados a partir dos dados de treinamento.O programa de demonstração carrega os valores de x preditor dos dados de treinamento em uma matriz e carrega os variável dependente y valores de dados em um vetor.Observe que a matriz de treinamento X, muitas vezes chamada de uma matriz de projeto, possui uma coluna adicional de primeira que consiste de todos os valores 1.0, e que os valores de prognósticos têm sido convertidos em valores numéricos.Em seguida, o programa de demonstração define três condições de parada do algoritmo IRLS, indicado pelo jumpFactor, epsilon e variáveis maxIterations.O programa de demonstração usa o algoritmo de Newton-Raphson para estimar a b0, b1, b2 e b3 valores de beta que melhor se ajustam aos dados de treinamento.O demo conclui, avaliando como a equação resultante de LR com valores computados beta é precisa dos dados de teste.Neste exemplo, 18 dos 20 valores de Y (90 por cento) foram previu corretamente.

Este artigo pressupõe que você tem avançadas habilidades de programação e, pelo menos, um conhecimento intermediário de terminologia e operações de matriz, mas não assuma que você sabe alguma coisa sobre regressão logística.O código que produziu a captura de tela em Figura 2 é demasiado grande para apresentar na íntegra aqui, mas o código fonte completo está disponível em archive.msdn.microsoft.com/mag201208TestRun.Devido à complexidade do algoritmo IRLS/NR, vou me concentrar principalmente em peças-chave do algoritmo em vez de no próprio código, então você será capaz de modificar o código para satisfazer suas próprias necessidades ou refatorá-lo para outra linguagem de programação, se desejar.

Estrutura geral do programa

Por simplicidade, todo o código que gerou a imagem em Figura 2 está contido em um único aplicativo de console em c#.A estrutura do programa e o método Main, com algumas declarações de WriteLine removidos, são listados em Figura 3.

Figure 3 estrutura do 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

Embora a regressão logística é um tema complexo, o código em Figura 3 não é completamente tão complicado como poderia parecer primeiro porque a maioria dos métodos mostrados são rotinas auxiliares relativamente curto. Os dois principais métodos são ComputeBestBeta e ConstructNewBetaVector.

O método de MakeRawData gera um arquivo de dados de idade-sexo-colesterol-morte quase. A idade é um valor aleatório inteiro entre 35 e 80, sexo é m ou f com probabilidade igual, e o colesterol é um valor real semi-aleatório entre 0,1 e 9,9 que é baseado no valor atual de idade. A variável dependente de morte é calculada usando a equação de regressão logística com valores fixos de beta de b0 =-95.0, b1 = 0,4, b2 = -0,9 e b3 = 11,2. Então MakeRawData gera dados que definitivamente podem ser modelados usando o LR, em vez de gerar dados puramente aleatórios que seriam provavelmente não seguem um modelo de LR.

Um novo vetor de Beta de computação

Coração de regressão logística com Newton-Raphson é uma rotina que calcula um novo, supostamente melhor, conjunto de valores beta do conjunto atual de valores. A matemática é muito profunda, mas felizmente o resultado não é muito complexo. Em forma de pseudo-equation, o processo de atualização é dada por:

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

Eis o novo b [t] ("no tempo t," não matriz indexação) vector de beta. Do lado direito, b [t-1] é o vetor de beta a idade ("no tempo t-1"). A função inv é a inversão de matriz. X maiúsculo é a matriz do projeto — ou seja, os valores das variáveis predictor aumentada com uma coluna à esquerda de 1.0s. Maiúscula X' é a transposta da matriz x design. Maiúsculas y é o vetor de valores de variável dependente (Lembre-se de que cada valor será 0 ou 1). A representa a quantidade p [t-1] vetor previu antigos valores de probabilidade para Y (que será composto de valores entre 0,0 e 1,0). A quantidade de w maiúscula é uma matriz de pesos chamados, o que requer um pouco de explicação.

A equação de atualização do beta e a matriz w são melhor explicadas com um exemplo concreto. Suponha que para manter a simplicidade que o conjunto de treinamento consiste somente as cinco primeiras linhas de dados mostrados na Figura 1. Assim, a matriz de projeto x seria:

    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

O vetor de variável dependente y seria:

0
1
0
1
0

Vamos supor que o vetor antigo de beta de valores a serem atualizados, b [t-1], é:

    1.00
    0.01
    0.01
    0.01

Com esses valores para x e beta, o vector p velho, p [t-1], é:

    0.8226
    0.8428
    0.8242
    0.8512
    0.8085

Observe que, se nós assumimos que os valores de p < 0,5 são interpretados como y = 0 e p valores > = 0,5 são interpretados como y = 1, os valores antigos de beta seriam prever correctamente apenas dois dos cinco casos nos dados de treinamento.

A matriz de pesos w é um m x matriz m, onde m é o número de linhas de X. Todos os valores da matriz w são 0,0, exceto para os valores de m que são na diagonal principal. Cada um desses valores é igual ao valor de p correspondente, multiplicado por 1-p. Assim, para este exemplo, W seria Tamanho 5 x 5. A célula superior esquerda em [0,0] seria (0.8226)(1-0.8226) = 0.1459. A célula no [1,1] seria (0.8428)(1-0.8428) = 0,1325 e assim por diante. A quantidade (p)(1-p) representa a derivada do cálculo da função sigmóide.

Na prática, a matriz w não é calculada explicitamente porque seu tamanho pode ser enorme. Se você tivesse 1.000 linhas de dados de treinamento, matriz w tenho 1.000.000 células. Observe que a equação de atualização do beta tem um termo W [t-1] X, o que significa o produto matricial de W [t-1] e X. Porque a maioria dos valores de W [t-1] são zero, a maioria dos termos de multiplicação de matriz são também zero. Isto permite W [t-1] vezes X, para ser calculado diretamente a partir de p [t-1] e X, sem explicitamente construindo W. Várias das referências matemáticas que descrevem IRLS com o algoritmo NR para LR usam o símbolo X ~ (X til) para o produto W [t-1] e X. Consulte o método ComputeXtilde no download de código para obter detalhes de implementação.

Método ConstructNewBetaVector aceita como parâmetros de entrada antigo vetor beta, a matriz de projeto X, o vetor de variável dependente y e o vetor de probabilidade a velha. O método calcula e retorna o vetor de beta atualizada.

O método é implementado da seguinte forma:

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;

Com a coleção de matriz e vetor auxiliar métodos mostrado na Figura 3, um novo vetor de beta de computação é bastante simples. Observe que o método realiza inversão de matriz. Este é um processo que pode dar errado em muitos aspectos e é uma fraqueza significativa de Newton-Raphson.

Continuando o exemplo, a matriz Xt (a transposta de X) seria:

     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 A (X ~) iria ser computada a partir de p de vetor e 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

Intermediário matriz B, que representa o produto da Xt e X ~ (que, por sua vez, é XtW[t-1]X) seria:

     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

Intermediário matriz c é a inversa da matriz b e seria:

     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

Intermediário matriz d é o produto da matriz c e matriz X-transpor e poderia ser calculado 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

Intermediate vector YP é a diferença entre vetores y e p [t-1] e seria:

    -0.8
     0.2
    -0.8
     0.1
    -0.8

Vetor intermediário e é o produto da matriz d e vector YP e mantém os incrementos a ser adicionado para o vetor antigo da beta. Vector e seria:

     4.1
    -0.4
    -2.8
     2.3

O vetor de novo, final beta é obtido adicionando-se que os valores no vetor intermediário e para os valores antigos do beta e neste exemplo seria:

     5.1
    -0.3
    -2.8
     2.4

Com os novos valores de beta, os novos valores para o vetor p seria:

    0.0240
    0.9627
    0.0168
    0.9192
    0.0154

Se esses valores de p são interpretados como Y = 0 quando p < 0,5, então, depois de apenas uma iteração de Newton-Raphson, os valores de beta prevêem correctamente todos os cinco casos nos dados de teste.

Saber quando parar

A técnica de Newton-Raphson para regressão logística iterativamente melhora os valores do vetor beta até que alguma condição de parada. É surpreendentemente difícil saber quando parar a iteração. ComputeBestBeta método lida com essa tarefa. O código é apresentado na Figura 4.

Figura 4 o melhor vetor de Beta de computação

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

Uma das maiores armadilhas da regressão logística é iterar muitas vezes, resultando em um conjunto de valores de beta que ajustar os dados de treinamento perfeitamente, mas não cabem quaisquer outros conjuntos de dados também. Isso é chamado de over-fitting. Saber quando parar o processo de formação em regressão logística é a arte e a ciência da parte. Método ComputeBestBeta começa por inicializar o vetor de beta para todos os valores 0.0, os valores de p associado de computação e informática, em seguida, o erro entre os valores de p e os valores de Y. O código neste artigo utiliza o quadrado médio do erro — a média da soma das diferenças ao quadrado entre p e Y. Outras possibilidades de erro de computação incluem desvio absoluto médio e erro de Cruz-entropia.

O loop de processamento principal em ComputeBestBeta repetidamente calcula um novo vetor de beta e o termo de erro correspondente. A condição de parada principal no algoritmo é um parâmetro de maxIterations. Lembre-se de que o algoritmo de Newton-Raphson calcula a matriz inversa, que é suscetível à falha. Nessa caso ComputeBestBeta retorna o vetor de beta melhor encontrado (que, infelizmente, poderia ser o vetor inicial de all-zero). Você pode querer lançar uma exceção aqui, em vez disso. Outra alternativa é tentar uma fuga, modificando os novos valores de beta, ajustando-os para a média dos valores anteriores e os novos valores.

A condição de parada próxima verifica para ver se a alteração em todos os novos valores de beta é menor do que algum valor pequeno de parâmetro epsilon, usando o método auxiliar NoChange. Isso indica que a Newton-Raphson tem convergido e, na verdade, há uma boa chance que seu modelo é agora over-fitted. Em vez de retornar o melhor vetor de beta encontrado neste momento, você talvez queira retornar o melhor vetor de beta de uma iteração anterior. A condição de parada seguinte verifica se qualquer um dos valores de beta novo saltaram de um montante maior do que alguns grande valor fornecido pelo parâmetro jumpFactor. Isso indica que Newton-Raphson possivelmente tem girado fora de controle e você deseja lançar uma exceção em vez de nova sintonização o melhor vetor de beta encontrado.

Outra condição de parada em ComputeBestBeta verifica ver se os novos valores de beta realmente produzem um erro maior do que os valores do beta anterior fez. Neste exemplo, se betas novos produzem um erro maior quatro vezes em uma fileira, o algoritmo termina e retorna o vetor de beta melhor encontrado. Você pode querer parametrizar o valor para o número máximo de aumentos consecutivos em erro. Quando é detectado um aumento no erro, o método tenta escapar da situação, alterando os valores do vetor de beta atual para a média dos valores na atual versão beta e os valores no beta calculado.

Não há nenhum único melhor conjunto de condições de parar. Cada problema de regressão logística requer um pouco de experimentação para sintonizar o algoritmo de Newton-Raphson.

Vantagens vs. Desvantagens

Existem várias alternativas de Newton-Raphson para estimar o melhor conjunto de valores de beta em regressão logística. Newton-Raphson é uma técnica de otimização numérica determinística. Você também pode empregar não-determinísticas (probabilística ou seja) técnicas como otimização de enxame de partículas, algoritmos evolutivos e otimização de forrageamento bacteriana.

A principal vantagem de Newton-Raphson comparado a alternativas probabilísticas é que, na maioria das situações, NR é muito mais rápido. Mas Newton-Raphson tem várias desvantagens. Porque NR usa a inversão de matriz, o algoritmo irá falhar quando se deparar com uma matriz singular durante o cálculo. Simples implementações de Newton-Raphson exigem que todos os dados em memória, o que significa que há um limite para o tamanho da formação e teste de conjuntos de dados, que você pode usar. Embora seja raro, alguns conjuntos de dados não podem convergir para um conjunto de melhores valores de beta, usando NR.

Quando eu uso de regressão logística, emprego, frequentemente, um dupla de ataque. Eu primeiro experimento com uma abordagem de Newton-Raphson. Se essa técnica falha, eu caio de volta usando otimização de enxame de partículas para encontrar o melhor conjunto de valores de beta. É importante notar que a regressão logística não é mágica, e nem todos os dados se ajusta a um modelo de regressão logística. Outras técnicas de aprendizado de máquina para o modelo de dados com uma variável dependente binária incluem redes neurais, suporte a máquinas de vector e análise discriminante linear.

A explicação do processo de atualização do vetor de beta do algoritmo Newton-Raphson apresentado no presente artigo e o download de código que acompanha você deve se levantar e correr com regressão logística usando NR. A regressão logística é um tema complexo, fascinante e pode fazer uma adição valiosa ao seu conjunto de habilidades pessoais.

Dr. James McCaffrey trabalha para a Volt Information Sciences Inc., onde gerencia o treinamento técnico para engenheiros de software trabalham em Microsoft Redmond, Wash., campus. Ele trabalhou em vários produtos da Microsoft, incluindo o Internet Explorer e o MSN Busca e McCaffrey é autor do livro “.NET Test Automation Recipes” (Apress, 2006). Ele pode ser contatado em jammc@microsoft.com.

Agradecemos aos seguintes especialistas técnicos pela revisão deste artigo: Dan Liebling e Anne Loomis Thompson