多クラスロジスティック回帰にRを使用する

Question

短い形式：

Rの勾配降下を介してマルチクラスロジスティック回帰分類アルゴリズムを実装する方法は？ 2つ以上のラベルがある場合、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 例、1から10までの手書き桁の20×20画像の400個の画素に起こるそれぞれ含む400の機能、を有する行列であります9：

これら400個のピクセルの値の「コンピュータビジョン」に基づき、手書きの数を予測するロジスティック回帰アルゴリズムを適用するには、二分決定されていないの新たな課題が発生します。係数を最適化することは、このR-bloggers exampleのようにアドホック勾配降下ループでは効率的ではありそうにありません。

R-bloggersにも、2つの説明変数（機能）と2つの結果に基づいてうまく機能する例があります。この例では、 R関数を使用していますが、これはthe way to goと思われます。

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! 

明らかに、これは不完全なようだ、と私にエラーを与える：私はドキュメントを読んでいるにもかかわらず

、私たちは10の可能な結果の中から決定する必要があり、この複雑な例を、設定に問題が発生しますメッセージ：

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

注意：X%*%theta_matrixは問題なく実行されます。したがって、問題は10個の分類器（0〜9）があり、関数costを使用して演算を実現できるように10 y列ベクトルの行列を作成しなければならないということです。上記の私の非動作コードと同じように、y = as.matrix(model.matrix(lm(y ~ as.factor(y))))のような行を含むダミーコードのyベクトルをダミーコードにすることも可能ですが、これもまた「one-versus-all」アイデアをカプセル化しているかどうかはわかりません - OK、おそらくそうではなく、おそらくこれが問題です。

それ以外の場合は、バイナリクラシファイアと非常に平行したコードでR-bloggers postで動作するようです。

この問題の正しい構文は何ですか？

I have tried to work it out one digit against all othersに注意してください。ただし、複雑さの点では意味がないと思います。

Answer 1

に入力するthetaはベクターである必要があります。コスト関数内で行列に変換することができます。

ここで前の質問をご覧ください：How to get optim working with matrix multiplication inside the function to be maximized in R