機械学習手習い: 多項式回帰と正則化
「入門 機械学習」手習い、6日目。「6章 正則化:テキスト回帰」です。
多項式回帰と、過学習を避けるための正則化について学び、最後に正則化を使って書籍の裏表紙の紹介文から人気順を予測します。
# 前準備 > setwd("06-Regularization/") > library('ggplot2')
非線形データの回帰分析
世の中には、直線では関係をうまく表現できないデータがあります。 例えばこんなの。
> set.seed(1) > x <- seq(-10, 10, by = 0.01) > y <- 1 - x ^ 2 + rnorm(length(x), 0, 5) > ggplot(data.frame(X = x, Y = y), aes(x = X, y = Y)) + geom_point() + geom_smooth(method = 'lm', se = FALSE) > ggsave(filename = "01.png")
回帰直線も引いてみましたが、構造をうまくとらえられていなのは明らかです。
このようなデータでも、一般化加法モデルを使うと、非線形な滑らかな表現を得ることができす。
geom_smooth を引数なしで使うと、一般化加法モデルでの回帰線が出力されます。
> ggplot(data.frame(X = x, Y = y), aes(x = X, y = Y)) + geom_point() + geom_smooth(se = FALSE) > ggsave(filename = "02.png")
さて、このような一見、直線には当てはめにくいデータですが、入力を変形させることで、線形回帰の前提を満たすようにすることができたりします。 ↑のデータであれば、xの2乗を使う形に変形することで、x,yの関係を単純な線形問題に変換できます。
> x.squared <- x ^ 2 > ggplot(data.frame(XSquared = x.squared, Y = y), aes(x = XSquared, y = Y)) + geom_point() + geom_smooth(method = 'lm', se = FALSE) > ggsave(filename = "03.png")
変形により、回帰線がどの程度データを表現できるようになったか、R2値を比較して見ます。
> summary(lm(y ~ x))$r.squared [1] 2.972904e-06 > summary(lm(y ~ x.squared))$r.squared [1] 0.9705898
変形前は、ほぼ0%であったのに対して、変形後は、97%を説明できるようになっています。
多項式回帰
データから複雑な形のモデルを構築するアプローチは、多項式回帰と呼ばれます。
簡単な直線では表現できない、正弦波のデータを使って、試してみます。
> set.seed(1) > x <- seq(0, 1, by = 0.01) > y <- sin(2 * pi * x) + rnorm(length(x), 0, 0.1) > df <- data.frame(X = x, Y = y) > ggplot(df, aes(x = X, y = Y)) + geom_point() > ggsave(filename = "04.png")
初めに、線形回帰で分析を行ってみます。
> summary(lm(Y ~ X, data = df)) Call: lm(formula = Y ~ X, data = df) Residuals: Min 1Q Median 3Q Max -1.00376 -0.41253 -0.00409 0.40664 0.85874 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.94111 0.09057 10.39 <2e-16 *** X -1.86189 0.15648 -11.90 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4585 on 99 degrees of freedom Multiple R-squared: 0.5885, Adjusted R-squared: 0.5843 F-statistic: 141.6 on 1 and 99 DF, p-value: < 2.2e-16
R2(Multiple R-squared) は0.58となっていて、60%の変動は表現できているようですが・・・
回帰線を引いてみます。
> ggplot(data.frame(X = x, Y = y), aes(x = X, y = Y)) + geom_point() + geom_smooth(method = 'lm', se = FALSE) > ggsave(filename = "05.png")
右肩下がりの直線になっています。正弦波をさらにもう一周伸ばすと、破綻しますね。 やはり、直線では無理があります。
ということで、まずは、xの2乗、3乗を入力に追加してみます。
> df <- transform(df, X2 = X ^ 2) > df <- transform(df, X3 = X ^ 3) > > summary(lm(Y ~ X + X2 + X3, data = df)) Call: lm(formula = Y ~ X + X2 + X3, data = df) Residuals: Min 1Q Median 3Q Max -0.32331 -0.08538 0.00652 0.08320 0.20239 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.16341 0.04425 -3.693 0.000367 *** X 11.67844 0.38513 30.323 < 2e-16 *** X2 -33.94179 0.89748 -37.819 < 2e-16 *** X3 22.59349 0.58979 38.308 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1153 on 97 degrees of freedom Multiple R-squared: 0.9745, Adjusted R-squared: 0.9737 F-statistic: 1235 on 3 and 97 DF, p-value: < 2.2e-16
R2の値が、58%から97%に向上しました。
このような感じで入力データを追加していけば、精度はさらに向上するらしい。 ただし、単純にxのべき乗を追加すると、新しく追加したべき乗の行と既存の行の相関が高くなってしまい、線形回帰のアルゴリズムが正常に動作しなくなるらしい。 しかし、ここで、xの連続したべき乗のようにふるまうが、それぞれの相関がないようにしたxの変種を追加するようにすれば、この問題を回避できるらしい。
このxの変種は、直行多項式と呼ばれ、 poly 関数を使って生成できます。
# 14個のxのべき乗を使う例 > summary(lm(Y ~ poly(X, degree = 14), data = df)) Call: lm(formula = Y ~ poly(X, degree = 14), data = df) Residuals: Min 1Q Median 3Q Max -0.232557 -0.042933 0.002159 0.051021 0.209959 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.010167 0.009038 1.125 0.2638 poly(X, degree = 14)1 -5.455362 0.090827 -60.063 < 2e-16 *** poly(X, degree = 14)2 -0.039389 0.090827 -0.434 0.6656 poly(X, degree = 14)3 4.418054 0.090827 48.642 < 2e-16 *** poly(X, degree = 14)4 -0.047966 0.090827 -0.528 0.5988 poly(X, degree = 14)5 -0.706451 0.090827 -7.778 1.48e-11 *** poly(X, degree = 14)6 -0.204221 0.090827 -2.248 0.0271 * poly(X, degree = 14)7 -0.051341 0.090827 -0.565 0.5734 poly(X, degree = 14)8 -0.031001 0.090827 -0.341 0.7337 poly(X, degree = 14)9 0.077232 0.090827 0.850 0.3975 poly(X, degree = 14)10 0.048088 0.090827 0.529 0.5979 poly(X, degree = 14)11 0.129990 0.090827 1.431 0.1560 poly(X, degree = 14)12 0.024726 0.090827 0.272 0.7861 poly(X, degree = 14)13 0.023706 0.090827 0.261 0.7947 poly(X, degree = 14)14 0.087906 0.090827 0.968 0.3358 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.09083 on 86 degrees of freedom Multiple R-squared: 0.986, Adjusted R-squared: 0.9837 F-statistic: 431.7 on 14 and 86 DF, p-value: < 2.2e-16
poly を使うことで、複雑な形を捉えることができるようになりますが、問題もあります。
これを確認するため、次数(degree)を増やしたモデルの形を出力してみます。
# 次数1 > poly.fit <- lm(Y ~ poly(X, degree = 1), data = df) > df <- transform(df, PredictedY = predict(poly.fit)) > ggplot(df, aes(x = X, y = PredictedY)) + geom_point() + geom_line() > ggsave(filename = "11.png")
# 次数3 > poly.fit <- lm(Y ~ poly(X, degree = 3), data = df) > df <- transform(df, PredictedY = predict(poly.fit)) > ggplot(df, aes(x = X, y = PredictedY)) + geom_point() + geom_line() > ggsave(filename = "12.png")
# 次数5 > poly.fit <- lm(Y ~ poly(X, degree = 5), data = df) > df <- transform(df, PredictedY = predict(poly.fit)) > ggplot(df, aes(x = X, y = PredictedY)) + geom_point() + geom_line() > ggsave(filename = "13.png")
# 次数25 > poly.fit <- lm(Y ~ poly(X, degree = 25), data = df) > df <- transform(df, PredictedY = predict(poly.fit)) > ggplot(df, aes(x = X, y = PredictedY)) + geom_point() + geom_line() > ggsave(filename = "14.png")
次数3～5では問題なさそうですが、次数25になると歪みはじめていて、元の波形を正しく表現できなくなっています。 この問題は、過学習と呼ばれます。
過学習を防ぐ
交差検定
交差検定は、データを2つに分割し、その一部を使って解析したモデルを、残りに当てはめて検証し、モデルの妥当性を検証する手法です。
交差検定を使って、次数の妥当性を検証してみます。 まずは、正弦波データを再作成。
> set.seed(1) > x <- seq(0, 1, by = 0.01) > y <- sin(2 * pi * x) + rnorm(length(x), 0, 0.1)
今回はデータを半分に分割して。片方を訓練データ、もう一方をテストデータとします。
今回のデータは単純に前半/後半で分割すると傾向が偏るので、sampleを使ってランダム抽出を行います。
> n <- length(x) > # 1～nのベクトルから、0.5 * n 個のインデックスをランダム抽出する > indices <- sort(sample(1:n, round(0.5 * n))) # インデックスのデータを訓練データとする > training.x <- x[indices] > training.y <- y[indices] # 残りはテストデータ > test.x <- x[-indices] > test.y <- y[-indices] > > training.df <- data.frame(X = training.x, Y = training.y) > test.df <- data.frame(X = test.x, Y = test.y)
精度は平均二乗誤差の平方根(RMSE)で確認します。RMSEを算出する関数を定義。
rmse <- function(y, h) { return(sqrt(mean((y - h) ^ 2))) }
準備ができたので、1～12まで次数でモデルを作成し、検証してみます。
> performance <- data.frame() > for (d in 1:12) { poly.fit <- lm(Y ~ poly(X, degree = d), data = training.df) performance <- rbind(performance, data.frame(Degree = d, Data = 'Training', RMSE = rmse(training.y, predict(poly.fit)))) performance <- rbind(performance, data.frame(Degree = d, Data = 'Test', RMSE = rmse(test.y, predict(poly.fit, newdata = test.df)))) } > ggplot(performance, aes(x = Degree, y = RMSE, linetype = Data)) + geom_point() + geom_line() > ggsave(filename = "20.png")
次数が3あたりから精度が高くなってきますが、11,12くらいになると逆に悪化してきています。 このことから、真ん中あたりの次数が妥当と判断できます。
正則化
モデルの精度と複雑さの妥協点を求めることで、過学習を防ぐアプローチもあります。
glmnet 関数を使うと、回帰に適用可能なすべての正則化セットが返されます。
Lambdaの数値が大きいほど、正則化の効果が強く(=複雑さが低く)なります。
各Lambdaごとに、精度を算出することで、妥協点を探ります。
> library('glmnet') > glmnet.fit <- with(training.df, glmnet(poly(X, degree = 10), Y)) > lambdas <- glmnet.fit$lambda > performance <- data.frame() # Lambdaごとに精度(RMSE)を計算 > for (lambda in lambdas) { performance <- rbind(performance, data.frame(Lambda = lambda, RMSE = rmse(test.y, with(test.df, predict(glmnet.fit, poly(X, degree = 10), s = lambda))))) } # 結果を図示 > ggplot(performance, aes(x = Lambda, y = RMSE)) + geom_point() + geom_line() > ggsave(filename = "21.png")
0.5あたりで、RMSEが最小となり、精度が一番高くなっているのがわかります。
テキスト回帰
多項式回帰と正則化の利用例として、オライリーが出版した書籍の売り上げトップ100冊の相対的な人気度を、 裏表紙の紹介文に含まれる単語から推定してみます。
単語を入力値とするようなテキスト回帰では、データの行より列の方が多く、正則化なしの線形回帰では常に過学習したモデルが生成されるらしい。 このため、何らかの正則化をする必要があるらしい。
まずはデータを読み込んで、文書単語行列を作ります。
> ranks <- read.csv(file.path('data', 'oreilly.csv'), stringsAsFactors = FALSE) > library('tm') > documents <- data.frame(Text = ranks$Long.Desc.) > row.names(documents) <- 1:nrow(documents) > > corpus <- Corpus(DataframeSource(documents)) > corpus <- tm_map(corpus, content_transformer(tolower)) > corpus <- tm_map(corpus, stripWhitespace) > corpus <- tm_map(corpus, removeWords, stopwords('english')) > dtm <- DocumentTermMatrix(corpus) # DTMをマトリックスに変換 > x <- as.matrix(dtm) # 順位を、100を最高位、1が最小位となるよう、逆順に変換。(見やすさのため) > y <- rev(1:100)
順位を計算して、結果を評価。
> performance <- data.frame() # 0.1, 0.25, 0.5, 1, 2, 5 の Lambda ごとに計測を実行 > for (lambda in c(0.1, 0.25, 0.5, 1, 2, 5)) { for (i in 1:50) { # 50回試行 # データを2つに分割 indices <- sample(1:100, 80) training.x <- x[indices, ] training.y <- y[indices] test.x <- x[-indices, ] test.y <- y[-indices] # 単語一覧から順位を推定する正則化セットを取得 glm.fit <- glmnet(training.x, training.y) # 指定したLambdaでの推定値を算出 predicted.y <- predict(glm.fit, test.x, s = lambda) # RMSE rmse <- sqrt(mean((predicted.y - test.y) ^ 2)) # 結果をデータフレームに格納 performance <- rbind(performance, data.frame(Lambda = lambda, Iteration = i, RMSE = rmse)) } } > head(performance) Lambda Iteration RMSE 1 0.1 1 34.85195 2 0.1 2 25.09364 3 0.1 3 26.98583 4 0.1 4 29.08673 5 0.1 5 32.75941 6 0.1 6 28.37690
図にしてみます。
> ggplot(performance, aes(x = Lambda, y = RMSE)) + stat_summary(fun.data = 'mean_cl_boot', geom = 'errorbar') + stat_summary(fun.data = 'mean_cl_boot', geom = 'point') > ggsave(filename = "30.png")
図をみれば明らかなように、Lambdaの値が大きくなるとモデルはよくなっているが、これはちょうどモデルが切片に縮退しているときに起こっているもので、テキストデータは全く使っていないことを示すらしい。端的に言えば、今回の回帰で発見された手がかりはないらしい。失敗らしい。
ロジスティック回帰
順位予測は失敗でしたが、書籍が50位に入るかどうかの予測ならできるかもしれません。 ロジスティック回帰は、ある要素が2つのカテゴリの一つに属する確率を予測する、回帰の一種です。
まずは、ロジスティック回帰をデータセット全体に適用し、どうなるか見てみます。
# 順位を0 or 1に分類 > y <- rep(c(1, 0), each = 50) # glmnetのfamilyを'binomial' にするとロジスティック回帰になる > regularized.fit <- glmnet(x, y, family = 'binomial')
モデルの予測がどうなるか見てみます。
> predict(regularized.fit, newx = x, s = 0.001) 1 1 4.573570 2 5.341867 3 4.668052 4 5.310955 5 4.869373 6 5.333862 7 4.487558 8 4.463377 9 5.352454 10 4.793254 11 5.352101 ..
正の値/負の値が含まれるので、ifelseを使って、0/1に変換すれば、期待する予測値が得られます。
> ifelse(predict(regularized.fit, newx = x, s = 0.001) > 0, 1, 0) 1 1 1 2 1 3 1 4 1 5 1
これを使って、書籍が50位以上に属するか推定してみます。
> set.seed(1) > performance <- data.frame() > for (i in 1:250) { indices <- sample(1:100, 80) training.x <- x[indices, ] training.y <- y[indices] test.x <- x[-indices, ] test.y <- y[-indices] for (lambda in c(0.0001, 0.001, 0.0025, 0.005, 0.01, 0.025, 0.5, 0.1)) { glm.fit <- glmnet(training.x, training.y, family = 'binomial') predicted.y <- ifelse(predict(glm.fit, test.x, s = lambda) > 0, 1, 0) error.rate <- mean(predicted.y != test.y) performance <- rbind(performance, data.frame( Lambda = lambda, Iteration = i, ErrorRate = error.rate)) } } > ggplot(performance, aes(x = Lambda, y = ErrorRate)) + stat_summary(fun.data = 'mean_cl_boot', geom = 'errorbar') + stat_summary(fun.data = 'mean_cl_boot', geom = 'point') + scale_x_log10() > ggsave(filename = "31.png")
lambaの値が小さいとき、偶然よりも良い予測になりました。
