故上述即为对数似然函数中19%的怀孕概率值的由来，19%亦作为is_pregnant的先验值。现在我有了所有参数的先验，可以建立一个由先验函数的抽样函数了。
sample_from_prior <- function(n) {
prior <- data.frame(mean_period = rnorm(n, 27.7, 2.4),
sd_period = abs(rnorm(n, 0, 2.05)),
is_fertile = rbinom(n, 1, 0.95))
prior$is_pregnant <- rbinom(n, 1, 0.19 * prior$is_fertile)
prior$next_period <- rnorm(n, prior$mean_period, prior$sd_period)
prior$next_period[prior$is_pregnant == 1] <- NA
prior
}
这里使用了一个参数(n)，它输出了一个n行的数据框，每一行是基于先验数值得出的样本数据。输出结果如下：
sample_from_prior(n = 4)
## mean_period sd_period is_fertile is_pregnant next_period
## 1 29 1.24 1 0 30
## 2 29 3.73 1 0 28
## 3 27 1.29 1 1 NA
## 4 27 0.57 0 0 27
现在，我已经收集了贝叶斯统计分析的三大要素：先验信息，似然函数以及数据。为了拟合模型我有很多方法，但是这里有一个非常方便的方法——重要性抽样。我之前曾写文提及过重要性抽样法，这里我们来回顾一下：重要性抽样法是一种蒙特卡洛实验法，它建立起来非常简单并且适用于以下情况：（1）参数空间非常小（2）先验分布与后验分布的形式区别不大。因为我的参数空间比较小，加之我使用了信息量包含得比较丰富的先验数据。因此，认为重点抽样法在此例中是可用的。在重要性抽样法中三个基本的步骤为:
1. 由先验分布产生大样本（这里可以通过sample_from_prior得到）
2. 给定了参数时，对每一个与似然值成比例的先验数据进行赋权。（这里可以通过 calc_log_like 得到）
3. 将权重归一化，从而在先验分布的情况下形成了新的概率分布。最终，根据此概率分布对先验分布的样本进行重新抽样。（这里可以用R函数抽样）
（ 注意存在与该过程不同的多种方法，但是在用来拟合贝叶斯模型时，这是重要性抽样法的常用版本）
因为我已经定义过 sample_from_prior 和 calc_log_like，因此需要定义一个新的方程来做后验抽样：
sample_from_posterior <- function(days_since_last_period, days_between_periods, n_samples) {
prior <- sample_from_prior(n_samples)
log_like <- sapply(1:n_samples, function(i) {
calc_log_like(days_since_last_period, days_between_periods,
prior$mean_period[i], prior$sd_period[i], prior$next_period[i],
prior$is_fertile[i], prior$is_pregnant[i])
})
posterior <- prior[ sample(n_samples, replace = TRUE, prob = exp(log_like)), ]
posterior
}
因此，在2月21日，2015，我的妻子已经没有来月经33天了。这是一个好休息吗？让我们运行这个模型看看结果吧！