Использование R для мульти-класса логистической регрессии

Question

Короткого формата:

Как реализовать логистические алгоритмы регрессионного классификации мульти-класса с помощью градиентного спуска в R? Можно ли использовать optim(), если имеется более двух ярлыков?

Код MatLab является:

function [J, grad] = cost(theta, X, y, lambda) 
    m = length(y); 
    J = 0; 
    grad = zeros(size(theta)); 
    h_theta = sigmoid(X * theta); 
    J = (-1/m)*sum(y.*log(h_theta) + (1-y).*log(1-h_theta)) +... 
    (lambda/(2*m))*sum(theta(2:length(theta)).^2); 
    trans = X'; 
    grad(1) = (1/m)*(trans(1,:))*(h_theta - y); 
    grad(2:size(theta, 1)) = 1/m * (trans(2:size(trans,1),:)*(h_theta - y) +... 
    lambda * theta(2:size(theta,1),:)); 
    grad = grad(:); 
end 

и ...

function [all_theta] = oneVsAll(X, y, num_labels, lambda) 
    m = size(X, 1); 
    n = size(X, 2); 
    all_theta = zeros(num_labels, n + 1); 
    initial_theta = zeros(n+1, 1); 
    X = [ones(m, 1) X]; 
    options = optimset('GradObj', 'on', 'MaxIter', 50); 
     for c = 1:num_labels, 
    [theta] = ... 
     fmincg (@(t)(cost(t, X, (y == c), lambda)), ... 
       initial_theta, options); 
    all_theta(c,:) = theta'; 
end 

---

Длинный формат:

Хотя, вероятно, не нужно следовать вопрос, набор данных может быть скачан here и загружен и помещен в каталог R, загружен как:

library(R.matlab) 
data <- readMat('data.mat') 
str(data) 
List of 2 
$ X: num [1:5000, 1:400] 0 0 0 0 0 0 0 0 0 0 ... 
$ y: num [1:5000, 1] 10 10 10 10 10 10 10 10 10 10 ... 

Так X представляет собой матрицу с 5000 примеров, каждая из которых содержит 400 особенности, которые происходят, чтобы быть 400 пикселей в 20 х 20 изображения рукописной цифры от 1 до 10, как, например, этой 9:

Применение логистической регрессии алгоритма для прогнозирования рукописного числа, основанное на «компьютерное зрении» значений в этих 400 пикселях представляет дополнительный вызов не является бинарным решением. Оптимизация коэффициентов вряд ли будет эффективной с ad hoc градиентной петли спуска, как в этом R-bloggers example.

Существует также хорошо разработанный пример также в R-bloggers на основе двух объясняющих переменных (признаков) и дихотомического результата. В примере используется функция optim() R, которая, как представляется, the way to go.

Даже если я прочитал документацию, у меня возникли проблемы с настройкой это более сложный пример, где мы должны решить среди 10 возможных исходов:

library(R.matlab) 
    data <- readMat('data.mat') 

    X = data$X     # These are the values for the pixels in all 5000 examples. 
    y = data$y     # These are the actual correct labels for each example. 
    y = replace(y, y == 10, 0) # Replacing 10 with 0 for simplicity. 

    # Defining the sigmoid function for logistic regression. 
     sigmoid = function(z){ 
      1/(1 + exp(-z)) 
     } 

    X = cbind(rep(1, nrow(X)), X) # Adding an intercept or bias term (column of 1's). 

    # Defining the regularized cost function parametrized by the coefficients. 

     cost = function(theta){ 
      hypothesis = sigmoid(X%*%theta) 
      # In "J" below we will need to have 10 columns of y: 
      y = as.matrix(model.matrix(lm(y ~ as.factor(y)))) 
      m = nrow(y) 
      lambda = 0.1 
      # The regularized cost function is: 
      J = (1/m) * sum(-y * log(hypothesis) - (1 - y) * log(1 - hypothesis)) + 
    (lambda/(2 * m)) * sum(theta[2:nrow(theta), 1]^2) 
      J 
     } 

    no.pixels_plus1 = ncol(X)  # These are the columns of X plus the intercept. 
    no.digits = length(unique(y)) # These are the number of labels (10). 
    # coef matrix rows = no. of labels; cols = no. pixels plus intercept: 
    theta_matrix = t(matrix(rep(0, no.digits*no.pixels_plus1), nrow = no.digits)) 
    cost(theta_matrix) # The initial cost: 
    # [1] 0.6931472 
    theta_optim = optim(par = theta_matrix, fn = cost) # This is the PROBLEM step! 

Очевидно, что это, кажется неполным, и дает мне ошибку сообщение:

Error in X %*% theta : non-conformable arguments 

Обратите внимание, что X%*%theta_matrix осуществляется без каких-либо проблем. Поэтому проблема заключается в том, что у меня есть 10 классификаторов (от 0 до 9) и что я вынужден создать матрицу с векторами столбцов 10 y, чтобы сделать операции выполнимыми с помощью функции cost. Возможно, что решение проходит через фиктивный код y с некоторой строкой вроде: y = as.matrix(model.matrix(lm(y ~ as.factor(y)))), как и в моем нерабочем коде выше, но опять же, я не знаю, что это инкапсулирует идею «один против всех» - Хорошо, наверное, нет, и, вероятно, это проблема.

В противном случае он работает на R-bloggers post с бинарным классификатором и чрезвычайно параллелен идентичному коду.

Итак, какой правильный синтаксис для этой проблемы?

Обратите внимание, что I have tried to work it out one digit against all others, но я не думаю, что это имеет смысл с точки зрения сложности.

Answer 1

theta вы кормите до optim должен быть вектором. Вы можете превратить его в матрицу внутри функции стоимости.

Смотреть предыдущий вопрос здесь: How to get optim working with matrix multiplication inside the function to be maximized in R