使用 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 是一个包含 5,000 个示例的矩阵，每个示例包含 400 个特征，它们恰好是 20 x 20 的手写数字图像的 400 像素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)))) 行，但我不知道这封装了“一对一” -all”的想法 - 好的，可能不是，这可能就是问题所在。

否则，它似乎可以在带有二进制分类器的R-bloggers post 上工作，并且与相同的代码极其平行。

那么这个问题的正确语法是什么？

请注意I have tried to work it out one digit against all others，但我认为这在复杂性方面没有意义。

Answer 1

您提供给optim 的theta 必须是向量。您可以将其转换为成本函数中的矩阵。

在此处查看上一个问题：How to get optim working with matrix multiplication inside the function to be maximized in R