Программирование логистической регрессии со стохастическим градиентным спуском в R

Я пытаюсь запрограммировать логистическую регрессию со стохастическим нисходящим градиентом в R. Например, я последовал примеру Эндрю Нг по имени: «ex2data1.txt».

Дело в том, что алгоритм работает правильно, но оценка теты не совсем то, что я ожидал. Поэтому я попытался изменить весь алгоритм, чтобы решить эту проблему. Однако для меня это было почти невозможно. Мне не удалось обнаружить ошибку, которая вызывает эту проблему. Таким образом, было бы очень полезно, если бы кто-нибудь мог проверить пример и сказать мне, почему тэты не рассчитываются правильно. Я очень ценю это.

Что касается программирования, я не использую никакие функции, реализованные в R, или вычисление матриц. Я просто использую суммы и вычитания в циклах, потому что я хочу использовать код в Hadoop, и я не могу использовать матричное исчисление или даже функции, которые уже запрограммированы в R, такие как «сумма», «кварт» и т. д.

Стохастический градиентный спуск:

Loop {
   for i = 1 to m, {
     θj := θj + α(y(i) - hθ(x(i)))(xj)(i)
  }
}`

Логистическая регрессия: ссылка на изображение

Мой код:

data1 <- read.table("~/ex2data1.txt", sep = ",")
names(data1) <- c("Exam1", "Exam2", "Admit")

# Sample the data for stochastic gradient decent 

ss<-data1[sample(nrow(data1),size=nrow(data1),replace=FALSE),]

x <- with(ss, matrix(cbind(1, Exam1), nrow = nrow(ss)))
y <- c(ss$Admit)
m <- nrow(x)

# startup parameters

iterations<-1
j<-vector()
alpha<-0.05
theta<-c(0,0)



#My loop

while(iterations<=10){

    coste<-c(0,0)
    suma<-0

    for(i in 1:m){

        # h<-1/(1+exp(-Q*x)

        h<-1/(1+exp((-theta)*x[i,]))

        #Cost(hQ(x),y)=y(i)*log(hQ(x))+(1-y(i))*log(1-hQ(x))

            cost<-((y[i]*log(h))+((1-y[i])*log(1-h)))

        #sum(cost) i=1 to m

            suma<-suma+cost

        #Diferences=(hQ(x(i))-y(i))*x(i)

            difference<-(h-y[i])*x[i,]  

        #sum the differences 

            coste<-coste+difference

        #calculation thetas and upgrade = Qj:= Qj - alpha* sum((h-y[i])*x[i,]*x(i))

            theta[1]<-(theta[1]-alpha*1/m*(coste[1]))
            theta[2]<-(theta[2]-alpha*1/m*(coste[2]))

    }
        #J(Q)=(-1/m)* sum ( y(i)*log(hQ(x))+(1-y(i))*log(1-hQ(x)))

            j[iterations]<-(-1/m)*suma

            iterations=iterations+1

}



#If I compare my thetas with R glm 


Call:  glm(formula = y ~ x[, 2], family = binomial("logit"), data = data1)

Coefficients:

Intercept:-4.71816 

x[, 2]  :0.08091  

Мои теты

Intercept: 0.4624024 
 x[,2]: 1.3650234

person user3488416    schedule 02.04.2014    source источник


Ответы (2)


Я реализовал решение в R для другого набора примеров Ng: ex2data2.txt. Вот мой код:

sigmoid <- function(z) {
return(1/(1 + exp(-z)))
}


mapFeature <- function(X1, X2) {
degree <- 6
out <- rep(1, length(X1))
for (i in 1:degree) {
for (j in 0:i) {
out <- cbind(out, (X1^(i - j)) * (X2^j))
}
}
return(out)
}


## Cost Function
fr <- function(theta, X, y, lambda) {
m <- length(y)
return(1/m * sum(-y * log(sigmoid(X %*% theta)) - (1 - y) *
log(1 - sigmoid(X %*% theta))) + lambda/2/m * sum(theta[-1]^2))
}


## Gradient
grr <- function(theta, X, y, lambda) {
return(1/m * t(X) %*% (sigmoid(X %*% theta) - y) + lambda/m *
c(0, theta[-1]))
}

data <- read.csv("ex2data2.txt", header = F)
X = as.matrix(data[,c(1,2)])
y = data[,3]
X = mapFeature(X[,1],X[,2])
m <- nrow(X)
n <- ncol(X)
initial_theta = rep(0, n)
lambda <- 1
res <- optim(initial_theta, fr, grr, X, y, lambda,
method = "BFGS", control = list(maxit = 100000))
person George Dontas    schedule 02.04.2014
comment
Привет! спасибо за ваш ответ, но я просто использую суммы и вычитания в циклах, потому что я хочу использовать код в hadoop, и я не могу использовать матричное исчисление или даже функции, которые уже запрограммированы в R, такие как сумма, sqrt, optim и т. д. - person user3488416; 03.04.2014
comment
Какова цель функции mapFeature? - person Lecromine; 24.07.2018
comment
Один из способов лучше подогнать данные — создать больше функций из каждой точки данных. В предоставленной функции mapFeature.m мы будем отображать функции во все полиномиальные члены x1 и x2 до шестой степени. В результате этого сопоставления наш вектор из двух признаков (оценок по двум тестам QA) был преобразован в 28-мерный вектор. Классификатор логистической регрессии, обученный на этом векторе признаков более высокого измерения, будет иметь более сложную границу решения и будет казаться нелинейным при рисовании на нашем двумерном графике. - person George Dontas; 24.07.2018

Разве * не должно быть %*% в некоторых случаях? Например. h<-1/(1+exp((-theta) %*% x[i,])) вместо h<-1/(1+exp((-theta)*x[i,]))

person Kai    schedule 18.08.2014