validation

Posted on Edited on

Validation for Bayesian model

一个非常不完善的总结/记录 for my painful discovery of Bayesian model validation
可能还有后续

说实话，model validation的总体思路就是split dataset into training and testing data. 然后利用training data在testing data的prediction去判断。有了这个整体思路，就是看需要解决的问题：

  • 怎么划分dataset：training/testing data的定义；在不同研究中结合实际。
  • 如何做training data在testing data的prediction：后验预测这一块（可能是该篇的重点内容）
  • 得到prediction之后怎么判断和testing data的observation：常用的validation判断，median/mean/95% predicted interval等
  • 如果没有达到预期怎么调整（目前我也不知道（笑）） 暂时不涉及

Definition of validation

虽然说到validation肯定都知道是什么，但是还是在这里贴一下（不一定）官方的：

Validation is the process of assessing how well a model performs on held-out data.

Why need validation? measure its predictive accuracy, for its own sake or for purposes of model comparison, selection, or averaging （Vehtari et al. 2017）.

Splitting dataset

首先，关于training/testing data的定义先奉上：

  • Training data: Data used to fit the model
  • Testing data: Data left-out during model fitting, used to evaluate model performance

非常的简单明了，但是，结合实际的时候就要分情况了。但是一般来说，都是80% training, 20% testing.
(这里是用前80%作为training, 后20%作为testing.) 但是！有时候划分不能很精确的按照这个比例，至于如何划分，也是一个值得研究的问题，但这里就不展开研究了 （因为我不会呢还。。。

Posterior predictive sample

Brief eg.

举例代码：

1
2
3
4
5
6
7
8
9
10
11
set.seed(123)
n <- 100
x.i <- 1:n
year.i <- seq(1990, 1990 + n - 1)
y.i <- x.i + rbeta(n, 1, 10) * 3 # 这个按理说是未知的，只给出y.i

data.full <- data.frame(x = x.i, year = year.i, y = y.i)

# 拆分 training 和 full dataset
train.cut <- quantile(data.full$year, 0.8) # 80% training
data.train <- subset(data.full, year <= train.cut)  # training set

这里用的jags建模（模型随意），然后parameters.to.save = c("beta0", "beta1")得到后验抽样：

1
2
3
b0.train <- fit.train$BUGSoutput$sims.list$beta0
b1.train <- fit.train$BUGSoutput$sims.list$beta1
sigma_train <- fit_train$BUGSoutput$sims.list$sigma

Posterior predict

后面就是对应testing data里面的每一条数据，找到对应year值匹配，得到一个matrix: row = length(testing data), col = predicted sample.

1
2
3
4
5
6
7
8
for (yr in data.valid.years) {
  row <- data.full[data.full$year == yr, ]
  x_val <- row$x
  y_obs <- row$y
  
  # Training 模型预测
  mu_draws  <- b0.train + b1.train * x_val
}

重点来了，pred.train在这里是期望值()，也就是latent mean，不是概率密度(posterior predictive density).
如果要得到后验概率likelihood：

1
2
# 把预测值当作分布的参数，去计算某个观测值在该分布下的密度或者概率
yrep <- rnorm(draw_length, mean = mu_draws , sd = sigma_train)

Summary

上面提到了1）预测的均值，是sample里面每一条数据的样本均值，用来获得posterior sample的均值预测区间等。2）分布概率，做 WAIC、LOO-CV、预测验证，需要结合 posterior 估计的 （精度）或（标准差）来计算。
比如这里，如果要和testing data (observation真实数据)对比，也就是预测验证，所以要考虑observation的噪声().

Validation

这里就是看predicted sample在testing data的预测情况，考虑下面两个指标：

  • Median/mean value
  • 95% PI

Error

这里包括Mean error 和 Absolute mean error.

1
2
3
4
5
6
# mean error
pred_mean[j]  <- mean(yrep)
residuals[j]  <- test$y[j] - pred_mean[j]

# abs mean error
mae <- mean(abs(residuals))

Coverage

如果选择95%PI，那么coverage的期望值就是95%（below下限和above上限的都是5%）。

1
2
3
4
pred_lower[j] <- quantile(yrep, 0.025)
pred_upper[j] <- quantile(yrep, 0.975)

coverage95 <- mean(test$y >= pred_lower & test$y <= pred_upper)

End

如果要看training data和full data训练的模型差异，就用full data fit相同的模型，然后用predicted value对比（不需要obs error）。使用同样指标。

在validation过程中，如果预测不达预期（error过大或coverage覆盖太少），首先考虑模型，error过大就是模型选择问题了，可以考虑变量相关（比如控制变量）；如果error很好但是coverage太小（低于60%这样），就要看看有没有考虑observation error！在这之前，可以先画一个PIT图，如果是U型，那就是观测误差的问题。 fit_train$BUGSoutput$sims.list$sigma在取后验样本加上这个东西。

关于pesterior predictive sample的生成，dnorm or rnorm 这个问题询问的是G老师，给出的回答：

rnorm: 后验预测分布, 直接地得到你要的验证指标
dnorm: 基于密度的评分规则（比如 log predictive density、elpd），
那时再用 dnorm(y_obs, mean=mu_draw, sd=sigma_draw) 求密度并取对数、对样本求平均就行

嗯，就是密度函数和随机生成函数的区别。具体的讲解这里讲的挺清楚的。

Reference

  1. Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical bayesian model evaluation using leave-one-out cross-validation and waic Statistics and computing, 27(5):1413–1432.
  2. R语言中dnorm, pnorm, qnorm与rnorm以及随机数。source