在R中实现逻辑回归公式

Question

在R中实现逻辑回归公式

6

我正在尝试使用随机梯度下降法在R中构建自己的逻辑回归函数，但目前的结果使权重无限增长，因此永远不会停止：

# Logistic regression
# Takes training example vector, output vector, learn rate scalar, and convergence delta limit scalar
my_logr <- function(training_examples,training_outputs,learn_rate,conv_lim) {
  # Initialize gradient vector
  gradient <- as.vector(rep(0,NCOL(training_examples)))
  # Difference between weights
  del_weights <- as.matrix(1)
  # Weights
  weights <- as.matrix(runif(NCOL(training_examples)))
  weights_old <- as.matrix(rep(0,NCOL(training_examples)))

  # Compute gradient
  while(norm(del_weights) > conv_lim) {

    for (k in 1:NROW(training_examples)) {
      gradient <- gradient + 1/NROW(training_examples)*
        ((t(training_outputs[k]*training_examples[k,]
            /(1+exp(training_outputs[k]*t(weights)%*%as.numeric(training_examples[k,]))))))
    }

    # Update weights
    weights <- weights_old - learn_rate*gradient
    del_weights <- as.matrix(weights_old - weights)
    weights_old <- weights

    print(weights)
  }
    return(weights)
}

该函数可通过以下代码进行测试：

data(iris) # Iris data already present in R    
# Dataset for part a (first 50 vs. last 100)
iris_a <- iris
iris_a$Species <- as.integer(iris_a$Species)
# Convert list to binary class
for (i in 1:NROW(iris_a$Species)) {if (iris_a$Species[i] != "1") {iris_a$Species[i] <- -1}}    
random_sample <- sample(1:NROW(iris),50)

weights_a <- my_logr(iris_a[random_sample,1:4],iris_a$Species[random_sample],1,.1)

我对我的算法进行了双重检查，与 Abu-Mostafa的算法如下：

初始化权重向量
对于每个时间步骤计算梯度：
gradient <- -1/N * sum_{1 to N} (training_answer_n * training_Vector_n / (1 + exp(training_answer_n * dot(weight,training_vector_n))))
weight_new <- weight - learn_rate*gradient
重复执行，直到权重差足够小

这里有什么遗漏吗？

- bright-star

我是否错过了权重的归一化项？这可能是一个交叉验证问题吗？ - bright-star

1

从数学角度来看，权重向量上的无约束大小并不产生唯一解。当我将这两行代码添加到分类器函数中时，它在两个步骤内收敛：weights <- weights/norm(weights)...weights <- weights_old - learn_rate*gradient weights <- weights / norm(weights) - bright-star

以下的回答有帮助吗？ - Simon O'Hanlon

2个回答

1

有几个问题。首先，您可以更好地利用R的矢量化方法。其次，我不是随机梯度下降的专家，但您在问题下方给出的算法与您在函数中计算梯度的方式不符。请仔细检查此代码，但它似乎会收敛，我认为它遵循Abu-Mostfafa的方法。我了解您想要计算这个的梯度;

gradient <- -1/N * sum(training_outputs * training_examples / (1 + exp(training_outputs * dot( weights ,training_outputs) ) ) )

这部分算法应该是这样写的...

while(norm(del_weights) > conv_lim) {  
gradient <- -1 / NROW(iris_a) * sum( training_outputs * training_examples / ( 1 + exp( training_outputs * as.matrix(training_examples) %*% weights ) ) )

# Update weights
weights <- weights_old - learn_rate*gradient
del_weights <- as.matrix(weights_old - weights)
weights_old <- weights
print(weights)

}

您可以更轻松地使用以下方法从物种变量创建二元分类：

iris_a$Species <- as.numeric( iris_a$Species )
iris_a$Species[ iris_a$Species != 1 ] <- -1

我无法告诉你返回的结果是否合理，但该代码应遵循第二步。仔细检查每一步，并记住 R 是向量化的，因此您可以在向量上进行元素操作而不需要循环。例如：

x <- 1:5
y <- 1:5
x*y
#[1]  1  4  9 16 25

- Simon O'Hanlon

那很有帮助（我也错误地屏蔽了其他二进制虹膜集，这里没有显示），但现在我又得到了另一个奇怪的结果，即每个分类器都以相同的方式训练：（http://i.imgur.com/qvamVQW.png） - bright-star

现在回归函数返回的所有权重似乎都是相同的。你确定那些乘法是正确的吗？ - bright-star

1

@TrevorAlexander 你好，那些乘法不正确，你发现了。我看到你发布了一个可行的解决方案，但我想回去整理向量化代码并使其工作，因为这是很好的练习！ - Simon O'Hanlon

我同意。我觉得最初接受C语言培训在进行线性代数编程时束缚了我。 - bright-star

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- bright-star · Accepted Answer

从数学角度来看，权重向量上的非约束性大小不会产生唯一的解。当我将这两行代码添加到分类器函数中时，它在两步内收敛：

# Normalize
weights <- weights/norm(weights)

...

# Update weights
weights <- weights_old - learn_rate*gradient
weights <- weights / norm(weights)

我无法理解 @SimonO101 的工作，而且我不打算在实际工作中使用此代码(因为有像 glm 这样的内置函数)，因此只需要编写我能理解的循环即可。

整个函数如下:

# Logistic regression
# Takes training example vector, output vector, learn rate scalar, and convergence delta limit scalar
my_logr <- function(training_examples,training_outputs,learn_rate,conv_lim) {
  # Initialize gradient vector
  gradient <- as.vector(rep(0,NCOL(training_examples)))
  # Difference between weights
  del_weights <- as.matrix(1)
  # Weights
  weights <- as.matrix(runif(NCOL(training_examples)))
  weights_old <- as.matrix(rep(0,NCOL(training_examples)))

  # Normalize
  weights <- weights/norm(weights)

  # Compute gradient
  while(norm(del_weights) > conv_lim) {

    for (k in 1:NCOL(training_examples)) {
      gradient <- gradient - 1/NROW(training_examples)*
        ((t(training_outputs[k]*training_examples[k,]
            /(1+exp(training_outputs[k]*t(weights)%*%as.numeric(training_examples[k,]))))))
    }
#     gradient <- -1/NROW(training_examples) * sum(training_outputs * training_examples / (1 + exp(training_outputs * weights%*%training_outputs) ) )

    # Update weights
    weights <- weights_old - learn_rate*gradient
    weights <- weights / norm(weights)
    del_weights <- as.matrix(weights_old - weights)
    weights_old <- weights

    print(weights)
  }
    return(weights)
}