2回目後半は線形分離できない2次元データのロジスティック回帰。
黄色○と+を識別したいが、どう見ても直線では識別できなさそう。
n次の空間に射影し、(n + 2)!/(n! * 2!)次元空間で線形分離する。
講義の問題と同様n = 6としてみる
R
ex2_2 <- function()
{
# ----- load data -----
data = read.csv('ex2data2.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))
legend("topleft", legend = c("y = 1", "y = 0"), pch = c(3, 16), col = c("black", "yellow"))
par(new = T)
# ----- variables -----
mappedX = mapFeature(X[, 1], X[, 2])
m = length(mappedX[,1])
n = length(mappedX[1,])
# ----- mapFeatureでバイアス項が足されている
initial_theta = rep(0, n)
# ----- optimize -----
print("start optimization")
lambda = 1
result <- optim(par = initial_theta, fn = costFunctionReg, gr = costGr, method = "CG", X = mappedX, y = y, lambda = lambda)
theta = result$par
plotDB2(theta, X, y)
p = round(sigmoid(mappedX %*% theta))
cat("Train Accuracy : ", mean(p == y) * 100, "\n")
return(result)
}
mapFeature <- function(X1, X2)
{
degree = 6
out = matrix(1, length(X1), 1)
for (i in 1 : degree) {
for (j in 0 : i) {
out = cbind(out, X1 ^ (i-j) * X2 ^ j)
}
}
return(out)
}
costFunctionReg <- function(theta, X, y, lambda)
{
m = length(y)
predicts = X %*% theta
probs = sigmoid(predicts)
costs = -y * log(probs) - (1 - y) * log(1 - probs)
J = sum(costs) / m
# ----- regularize -----
J = J + lambda / (2 * m) * sum(theta[-1] ^ 2) # theta[-1]は1番目以外の要素
return(J)
}
costGr <- function(theta, X, y, lambda)
{
m = length(y)
predicts = X %*% theta
probs = sigmoid(predicts)
grad = (t(X) %*% (probs - y)) / m + lambda / m * theta
grad[1] = (t(X)[1, ] %*% (probs - y)) / m
return(grad)
}
plotDB2 <- function(theta, X, y)
{
num_points = 50
u <- seq(0.9 * min(X[, 1]), 1.1 * max(X[, 1]), length = num_points)
v <- seq(0.9 * min(X[, 2]), 1.1 * max(X[, 2]), length = num_points)
z <- matrix(0, nrow = length(u), ncol = length(v))
for (i in 1 : length(u)) {
for (j in 1 : length(v)) {
z[i, j] = mapFeature(u[i], v[j]) %*% theta
}
}
contour(u, v, z, xlim = c(0.9 * min(X[, 1]), 1.1 * max(X[, 1])), ylim = c(0.9 * min(X[, 2]), 1.1 * max(X[, 2])))
}
optim関数にmethod = "CG"を指定しているのは、デフォルトだと識別があまりよくならなかったため。
結果
n = 2でやってみた
Train Accuracy : 81.35593
こっちでもいいんでないかという気もするレベル。
nを決めるには横軸n、縦軸TrainAccuracyのグラフを描いてみれば いいんだったか。
あるいは縦軸は汎化能力か
コスト関数の値と勾配を別の関数で計算しているせいで、無駄な計算をしている。
matlabの最適化関数では、関数の返り値を多次元にして勾配を返してやるとそれを使って最適化してくれていたのだけど、Rのoptimは同じことはできないらしい。
fn
A function to be minimized (or maximized), with first argument the vector of parameters over which minimization is to take place. It should return a scalar result.
