[R - ml] 模型的评估

主要的模型评估方式有三种:

1. hold out

以构建模型的训练数据来评估模型,显然不是最佳选择,因为我们的目标是对未来数据的分析。
而模型是以训练集来建立的。
所以要 train, test

然而我们不能重复的进行测试以选择最佳模型,因此常用的方式是
train, validataion, test
用validatation 数据集来帮助寻找最佳的模型,然后用test 完成最后一次测试,作为对未来数据的验证

random_ids = order(runif(1000))
train = credit[random_ids[1:500], ]
validate = credit[random_ids[501:750], ]
test = credit[random_ids[751:1000], ]

如果数据的class 是非均衡的,就可能有问题

分层随机采样

require(caret)
in_train = createDataPartition(credit$default, p = 0.75, list = FALSE) # 不以list保存
train = credit[in_train, ]
test = credit[-in_train, ]

2. Cross-validatation : CV

10-fold cv

  • 100 行数据
  • 20 行用于测试
  • 80 行用于训练模型
  • 例如 对于KNN 模型, 我们想找到最佳k
  • 对于每个k, 我们进行10-fold CV。即把训练数据分10组,取其中的9组(72行)构建模型,
    然后计算剩余的 8 行数据的预测错误,之后我们对10次计算后的错误平均。最终我们选择最小的那个模型。
    0.2, 0.3, 0.4, 0.5 取最小的那个k
  • 利用 剩余80行数据,以及前面选择的 k 来构建模型
  • 利用剩余的20行数据验证我们的模型

caret for CV

require(caret)
folds = createFolds(credit$Good.Loan, k = 10)
str(folds)
credit01_train = credit[ -folds$Fold01, ]
credit01_test = credit[folds$Fold01, ]

信用卡的例子:

library(caret)
library(C50)
library(irr)
set.seed(123)
folds = createFolds(credit$Good.Loan, k = 10)
cv_results = lapply(folds, function(x) {
  credit_train = credit[-x, ]
  credit_test = credit[x, ]
  credit_model = C5.0(Good.Loan ~. , data = credit_train)
  credit_pred = predict(credit_model, credit_test)
  credit_actual = credit_test$Good.Loan
  kappa = kappa2(data.frame(credit_actual, credit_pred))$value
  return(kappa)
})

str(cv_results)
mean(unlist(cv_results))

说明前面 我们关于信用的坏账模型不怎么样。
k-fold 选多大是一个bias 与 vias 的问题。
k 越小 显然 bias(偏差) 更大, variance(方差)小,
而k越大, bias 小, 但是 variance大。

3. boostrap

boostrap 源于统计中利用随机抽样来估计总体的方法。
类似机器学习中我们用多个随机选取的训练与测试及来估计模型的性能统计。
最终利用均值来作为未来新数据下的性能估计。
k-fold 与 bootstrap 的区别在于 bootstrap是有放回的,也就是说某些样本可能重复出现。

bootstrap 错误估计同时考虑过于乐观的训练集与过于悲观的测试集。
error = 0.632 × error_{test} + 0.368 × error_{train}

回归问题

MeanSqureErroe MSE = \frac{1}{N} {\sum{(y_{est}-y_{actual})^2}}
RootMeanSquareError RMSE = \sqrt{MSE}
Relative Squre Erroe RSE = \frac{MSE} {MSE_{baseline}} = \frac{ \sum{(y_{est}-y_{actual})^2}}{\sum{(y_{mean}-y_{actual})^2}}
R^2 模型的解释性:
R^2 = \frac{MSE_{baseline}-MSE}{MSE_{baseline}}
R^2 = 1- RSE

demo

require(car)
?Prestige
Prestige_clean = Prestige[!is.na(Prestige$type), ]
model = lm(prestige~. , data = Prestige_clean)
actual = Prestige_clean$prestige
rmse = (mean(score - actual)^2)^0.5
rmse
mu = mean(actual)
rse = sum(score - actual)^2 /sum(mu - actual)^2
rse
rsquare = 1 - rse
rsquare

MeanSquareLogError MSLE = \frac{1}{N} {\sum(logy_est - logy_actual)^2 }
RootMeanSquareLogError MSLE = \sqrt MSLE

分类问题

如何判断分类器性能?

实际分类值
预测分类
估计的预测分类概率

实际与预测对比

actual_outcom = test_data$out_come
predicted_outcome = predict(model, test_data)

预测分类的概率

R语言不同模型,参数不同。prob, posterior, raw, probability 都有可能

# naive bayes
predict(model, test_data, type = 'raw')
# decision tree
predict(model, test_data, type = 'prob')

Confusion Matrix

accuracy = \frac{TN + TN} {TP + TN + FP + FN}
erroratte =\frac{FP + FN} {TP + TN + FP + FN} = 1 - accuracy
table()

table(sms_raw_test$type, sms_test_pred)
require(gmodels)
CrossTable(sms_test_pred, sms_raw_test$type,
                     prop.chisq = FALSE, prop.t = FALSE,
                     dnn = c('predicted', 'actual'))

caret::ConfusionMatrix()

confusionMatrix(sms_test_pred, sms_raw_test$type, positive = 'spam')

次函数将输出Kappa
kappa statistic
.< 0.2 预测比差
.> 0.8 好
k = [Pr(a) - Pr(e)] / [1 - Pr(e)]

Kappa 统计 是相对预期一致的准确性调整
Pr(e) 是预期值与实际值一致的概率
Pr(e) = Pr(actual postive)Pr(predicted positive) + Pr(actual negative)Pr(predicted negative)

require(vcd)
table(sms_raw_test$type)
table(sms_raw_test$type, sms_test_pred)
kappa(table(sms_raw_test$type, sms_test_pred))
pr.ap = 186 / (186 + 1206)
pr.pp = (157 + 4) / (4 + 157 + 1202 + 29)
pr.an = 1206 / (1206 + 186)
pr.pn = (1202 + 29) / (1202 + 29 + 4 + 157)
pr.e = pr.ap*pr.pp+pr.an*pr.pn
pr.a = (1202 + 157) / (1202 + 29 + 4 + 157)
kappa.c = (pr.a - pr.e) / (1 - pr.e)
kappa.c

sensitivity 和 specificity

sensitivity 度量 postive 样本被正确分类的比率
specificity 度量 negative 样本被正确分类的比率

# manual calculation
sensitivity(sms_test_pred, sms_raw_test$type, positive = 'spam')
specificity(sms_test_pred, sms_raw_test$type, negative = 'ham')

0~1 越接近 1 越好
同时考虑的方式 trade-off

precision 和 recall

正确性: precision =\frac{TP}{TP + FP}
完备性: recall =\frac{TP}{TP + FN}
实际上 recall = sensitivity
posPredValue 由caret提供函数计算

posPredValue(sms_test_pred, sms_raw_test$type, positive = 'spam')

F - measure

F - measure 是以调和均值将precision和recall 以一个指标来衡量。
F - measure = \frac{ 2 * precision * recall} {precision + recall} = \frac{2TP} {2TP + FN + FP}

ROC 曲线

receiver operating characteristic curve

AUC

Area Under Curve 被定义为ROC曲线下的面积,显然这个面积的数值不会大于1

require(ROCR)
sms_test_pred = predict(sms_classifier, sms_test, type = 'raw')
pred = prediction(sms_test_pred[, 2], sms_raw_test$type)
perf = performance(pred, measure = 'tpr', x.measure = 'fpr')
plot(perf, main = "ROC curve", col = 'blue', lwd = 2)
abline(a = 0, b = 1, lwd = 2, lty = 2) # 两步要一起操作

AUC 可以通过以下方式获得:

perf.auc = performance(pred, measure = 'auc')
# performance 是s4 的一个对象
str(perf.auc)
unlist(perf.auc@y.values)
library(ROCR)
library(e1071)
inTrain = createDataPartition(y = iris$Species, p = 0.75, list = FALSE)
iristrain = iris[inTrain, ]
iristest = iris[-inTrain, ]

nb_model = naiveBayes(Species~. , data = iristrain)
nb_prediction = predict(nb_model, iristest[, -5], type = 'raw')
score = nb_prediction[, c('virginica')]
actual_class = iristest$Species == 'virginica'
# Mix some noise to the score
# Make the score less precise for illustration
score = (score + runif(length(score))) / 2
pred = prediction(score, actual_class)
perf = performance(pred, 'prec', 'rec')
plot(perf)

library(ROCR)
library(e1071)
nb_model = naiveBayes(Species~. , data = iristrain)
nb_prediction = predict(nb_model, iristest[, -5], type = 'raw')
score = nb_prediction[, c('virginica')]
actual_class = iristest$Species == 'virginica'
# Mix some noise to the score
# Make the score less precise for illustration
score = (score + runif(length(score))) / 2
pred = prediction(score, actual_class)
perf = performance(pred, 'prec', 'rec')
auc = performance(pred, 'auc')
auc = unlist(slot(auc, 'y.values'))
plot(perf)
legend(0.6, 0.3, c(c(paste('AUC is', auc),'\n'), border = 'white', cex = 1.0,
                   box.col = 'white'))

ROC - cost

根据FP, FN 的cost来选择最佳的cutoff点
cost.fp, cost.fn 惩罚因子

# Plot the cost curve to find the best cutoff
# Assign the cost for False + ve and False - ve
perf = performance(pred, 'cost', cost.fp = 4, cost.fn = 1) 
plot(perf)

Bias vs Variance

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,294评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,493评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,790评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,595评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,718评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,906评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,053评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,797评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,250评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,570评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,711评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,388评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,018评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,796评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,023评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,461评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,595评论 2 350

推荐阅读更多精彩内容