2回目のお題はロジスティック回帰
簡単な関数いくつか作って、後はＲの最適化関数を使うだけ
ex2 <- function()
{
# ----- load data -----
data = read.csv('ex2data1.txt', header = F)
X <- data[, 1:2]
y <- data[, 3]
# ----- plot -----
plot(0, 0, type = "n", xlim = c(0.9 * min(X[, 1]), 1.1 * max(X[, 1])), ylim = c(0.9 * min(X[, 2]), 1.1 * max(X[, 2])),xlab = "x1", ylab = "x2")
points(X, col = ifelse(y == 1, "black", "yellow"), pch = ifelse(y == 1, 3, 16))
# ----- variables -----
X_mat <- cbind(X[, 1],X[, 2])
m <- length(X[,1])
n <- length(X[1,])
# ----- add bias -----
initial_theta <- rep(0, n + 1)
X_mat <- cbind(rep(1, m), X_mat)
# ----- optimize -----
result <- optim(par = initial_theta, costFunction, X = X_mat, y = y)
theta <- result$par
plotDecisionBoundary(theta, X, y)
}
sigmoid <- function(x)
{
return(1 / (1 + exp(-x)))
}
costFunction <- function(theta, X, y)
{
m <- length(y)
predicts <- X %*% theta
probs <- sigmoid(predicts)
costs <- -y * log(probs) - (1 - y) * log(1 - probs)
J <- sum(costs) / m
return(J)
}
plotDecisionBoundary <- function(theta, X, y)
{
# ----- plot -----
plot(0, 0, type = "n", xlim = c(0.9 * min(X[, 1]), 1.1 * max(X[, 1])), ylim = c(0.9 * min(X[, 2]), 1.1 * max(X[, 2])),xlab = "x1", ylab = "x2")
points(X, col = ifelse(y == 1, "black", "yellow"), pch = ifelse(y == 1, 3, 16))
par(new=T)
x1 <- c(min(X[, 1]), max(X, 1))
x2 <- -1 / theta[3] *( theta[2] * x1 +theta[1])
plot(x1, x2, ,type = 'l', xlim = c(0.9 * min(X[, 1]), 1.1 * max(X[, 1])), ylim = c(0.9 * min(X[, 2]), 1.1 * max(X[, 2])))
}
----- 追記 -----
ちなみにglm関数を使うと
ex2_glm <- function()
{
# ----- load data -----
data = read.csv('ex2data1.txt', header = F)
X <- data[, 1:2]
y <- data[, 3]
# ----- GLM -----
fit <- glm(formula = y ~ V1 + V2 + 1, family = binomial, data = X)
theta <- fit$coef
plotDecisionBoundary(theta, X, y)
print(summary(fit))
return(fit)
}
