機械学習ということであるが，例のごとく「機械的に」train:test = 8:2 に分割し，1 回だけの分析である。
以下の 2 つの結論だけを述べている。
「調整済みの決定係数が 0.5777 であった。特別良くもないが悪くもない結果だった。」
「テストデータの RMSE を評価する。 8.143424」
機械学習は，もっと統計学的見地からの分析を踏まえてやっておかないと，足下をすくわれる。
今回のデータセットのサンプルサイズは 414 である。これを 8:2 に分割すると，テストデータはサンプルサイズが 83 しかなくなる。
わずか 83 例でモデルの評価ができるのか？？
しかも，このデータセットの従属変数には外れ値がある。全体からすれば 1/414 = 0.2% の外れ値であるが，この外れ値が 83 個のテストデータに含まれるような場合には，一挙に 1/83 = 1.2% を占めることになる。
そのため，評価に多かれ少なかれ影響を及ぼすことになる。
元のプログラムを Julia に書き換えると 30 倍速くなる。いいかえれば，同じ時間で R の30倍の回数のシミュレーションを行える。
30万回のシミュレーションを行うと，最初の図に示すように RMSE のヒストグラムを描くと右に瘤のある分布になっている。
この瘤の原因が 1 個の外れ値である。外れ値がテストデータに含まれるとこの瘤に含まれるのである。RMSE が 10.3 以上になったのが 約 20% である。
トレインデータで得られたモデルが有用であるとしても，テストデータであまり好ましくない評価を得られる確率が 20 % もあるのは甚だ不適切な話だ。
今回のデータは，トレインデータで得られたモデルもさして有用でないので，論外ではあるが。
この 1 個の外れ値を除いて同じく30万回のシミュレーションを行うと，RMSE は綺麗な分布を示す。
この後の議論は何段階かあるのだが，結論は以下のようになろう。
数少ないデータを細切れにしてモデルを作成して，信頼性を評価しようとしても，全体としてモデルの信頼性はさらに低くなるだけである。
逆に，十分な数のデータがあるときには，そこからどのようにデータを抽出・分割してもだいたい同じような結果が得られる。それが確率というものである。
機械的に機械学習を行ってはいけない。
#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
某 qiita の r タグにある機械学習のプログラム
library(rsample)
df <- read.csv("Real_estate_valuation.csv")
#coln <- colnames(df) #列名を保存しておく
#colnames(df) <- c("No","X1","X2","X3","X4","X5","X6","Y") #新しい名前を付けておく
set.seed(123)
sim = function() {
df_split <- initial_split(df, prop = 0.8) #train:test = 8:2 の割合で分ける
df_train <- training(df_split)
df_test <- testing(df_split)
df_train.lm <- lm(Y~X1+X2+X3+X4+X5+X6, df_train)
summary(df_train.lm)
predict2 = function(c, df){
for (i in names(c)){
if (i == "(Intercept)") { tmp <- c[i] }
else { tmp <- tmp + c[i]*df[i] }
}
return(tmp[,1])
}
RMSE = function(m, o) {
tmp <- sqrt(mean((m-o)^2))
return(tmp)
}
c <- coef(df_train.lm)
df_test.y <- predict2(c, df_test)
df_test.rmse <- RMSE(df_test.y, df_test$Y)
df_test.rmse
}
翻訳するときに書いたメモ
自前の predict2 は無駄である。既存の predict() を使えばよい。
==========#
using CSV, DataFrames, Query, GLM, Statistics, Plots, Printf
function sim(df, prop = 0.8)
n = size(df, 1)
m = round(Int, n * prop)
df_split = shuffle(collect(1:n));
df_train = df[df_split[1:m], :];
df_test = df[df_split[m+1:n], :];
df_train_lm = lm(@formula(Y~X1+X2+X3+X4+X5+X6), df_train);
df_test_y2 = predict(df_train_lm, df_test);
sqrt(mean((df_test_y2 .- df_test.Y) .^2)) # RMSE
end
function sim0(df, trial=10000)
n = size(df, 1)
a = lm(@formula(Y ~ X1+X2+X3+X4+X5+X6), df)
pred = predict(a);
rmse0 = sqrt(mean((pred .- df.Y) .^ 2))
R2 = r2(a)
R2adj = adjr2(a)
@printf("n = %d, mean(Y) = %.5g, RMSE = %.5g, R2 = %.5f, R2.adj = %.5f\n",
n, mean(df.Y), rmse0, R2, R2adj)
rmse = zeros(trial);
for i in 1:trial
rmse[i] = sim(df)
end
mean_rmse = mean(rmse)
median_rmse = median(rmse)
@printf("%d trials, Mean(RMSE) = %.5g, Median(RMSE) = %.5g)\n",
trial, mean_rmse, median_rmse)
pyplot()
histogram(rmse, bins=50, label="", xlabel="RMSE", ylabel="Frequency")
scatter!([rmse0], [0], label="\$RMSE_0\$")
scatter!([mean_rmse], [0], label="\$Mean(RMSE)\$")
scatter!([median_rmse], [0], label="\$Median(RMSE)\$")
#savefig(@sprintf "sim%d.png" n)
end
using Random
Random.seed!(123);
df = CSV.read("Real_estate_valuation.csv", DataFrame, drop=[1]);
# names(df) # "X1", "X2", "X3", "X4", "X5", "X6", "Y"
trial = 300000
# 外れ値を含んだままのデータの分析
sim0(df, trial)
n = 414, mean(Y) = 37.98, RMSE = 8.7825, R2 = 0.58237, R2.adj = 0.57621
300000 trials, Mean(RMSE) = 8.8087, Median(RMSE) = 8.3622)
# 外れ値を除いたデータの分析
df2 = copy(df);
delete!(df2, 271); # 外れ値除去
sim0(df2, trial)
n = 413, mean(Y) = 37.788, RMSE = 7.9608, R2 = 0.62675, R2.adj = 0.62123
300000 trials, Mean(RMSE) = 8.0752, Median(RMSE) = 8.0537)
