在2015年的二月21日，我的妻子已经33天没有来月经了，她怀孕了，这真是天大的好消息！通常月经的周期是大约一个月，如果你们夫妇打算怀孕，那么月经没来或许是一个好消息。但是33天，这还无法确定这是一个消失的月经周期，或许只是来晚了，那么它是否真的是一个好消息？
为了能获得结论我建立了一个简单的贝叶斯模型，基于这个模型，可以根据你当前距离上一次经期的天数、你历史经期的起点数据来计算在当前经期周期中你怀孕的可能性。在此篇文章中我将阐述我所使用的数据、先验思想、模型假设以及如何使用重点抽样法获取数据并用R语言运算出结果。在最后，我将解释为什么模型的运算结果最终并不重要。另外，我将附上简便的脚本以供读者自行计算.
数据非常幸运的是，在2014年的下半年间我的妻子一直在记录她经期起始日期，否则我只能以仅拥有小量数据而告终。总体上我们拥有8个经期的起始日期数据，但是我采用的数据不是日期而是相邻经期起始日间相隔的天数。 已经有33天。
period_onset <- as.Date(c("2014-07-02", "2014-08-02", "2014-08-29", "2014-09-25",
"2014-10-24", "2014-11-20", "2014-12-22", "2015-01-19"))
days_between_periods <- as.numeric(diff(period_onset))
所以日期发生得相对规律，以28天为一个周期循环。最后一次月经开始日期是在1月19日，所以在2月21日，距离最后一次经期发生日。
模型的建立我要建立一个涵盖生理周期的模型，包括受孕期和不受孕期，这显然需要做大量的简化。我做了一些总体假设如下：
一对情侣受孕与否不受其他因素的影响。
女方拥有固定的经期。
该对想要受孕的夫妻正在积极地尝试受孕。换言之，如Wilcox et al. (2000) 推荐的每周两次到三次受精。
一旦怀孕，期间将不会有经期。
接下来是我所做的具体假设：
假设两个相邻经期间隔的天数（days_between_periods）服从的正态分布，其中均值（mean_period）和标准差（sd_period）未知。
假设在一对生育能力的夫妻（is_fertile为真 ）受孕时一个周期内怀上孕的概率是0.19（更多关于选定该值的由来见参考文献）。不幸的是，并非所有的夫妻都具备生育能力，没有生育能力则怀孕的几率为零。如果生育率被编码为 0-1，那么可生育率可以被简洁的写为 0.19* is_fertile.
在某一些不能受孕的时期（n_non_pregnant_periods）的怀孕失败率则为(1 – 0.19 * is_fertile)^n_non_pregnant_periods
最后，如果你在这一个周期内（从上一次生理期至这一次生理期为一个周期）将不会怀孕；那么最新一次经期距离下一个经期的天数（next_period）将必然会大于最新一次经期距离当前日期的天数（days_since_last_perio）。即，next_period < days_since_last_period的概率为零。这么做看上去很奇怪因为这个事件是显然的，但是我们在模型中将会要用到它。
基本的假设就是这样了。但是为了使其更加实际，需要考虑使用一个似然函数，一个给定了参数和一些数据、计算在给定参数下数据的概率，通常而言是一个与概率成正比例的数值——似然值。因为这个似然值可能极小所以我需要对其取对数，从而避免引起数值问题。当用R语言设计似然函数时，总体上的模式如下：
方程将数据和参数作为选项。
通过预处理，将似然值的初始值设为1.0，相应的对数为0.0。（log_like <- 0.0）
用R语言调用概率密度分布函数（比如dnorm, dbinom and dpois），用该函数计算模型中不同部分的似然值。然后将这些似然值相乘。对应地，将取对数后的似然值log_like相加。
为了让d*函数返回对数似然值，只需添加参数log=TRUE。并且注意似然值0.0对应的取对值为-inf
calc_log_like <- function(days_since_last_period, days_between_periods,
mean_period, sd_period, next_period,
is_fertile, is_pregnant) {
n_non_pregnant_periods <- length(days_between_periods)
log_like <- 0
if(n_non_pregnant_periods > 0) {
log_like <- log_like + sum( dnorm(days_between_periods, mean_period, sd_period, log = TRUE) )
}
log_like <- log_like + log( (1 - 0.19 * is_fertile)^n_non_pregnant_periods )
if(!is_pregnant && next_period < days_since_last_period) {
log_like <- -Inf
}
log_like
}