A script that implements a Bayesian model calculating the probability that a couple is fertile and is going to be pregnant.

probability_of_pregnancy.R

# A Bayesian model that calculates a probability that a couple is fertile 
# and pregnant. Please use this for fun only, not for any serious purpose 
# like *actually* trying to figure out whether you are pregnant.

# Enter your own period onsets here:
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"))
# If you have no dates you can just set days_between_periods to c() instead like:
# days_between_periods <- c()
days_between_periods <- as.numeric(diff(period_onset))

days_since_last_period <- 33

# The prior probability the couple is fertile
prior_is_fertile <- 0.95

# The probability of becoming pregnant any given cycle
# You can adjust this according to the age of the woman in this table
#   19-26 years 27-34 years 35-39 years 
#     0.25      0.19      0.16 
prior_is_pregnant <- 0.19


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 - prior_is_pregnant * is_fertile)^n_non_pregnant_periods )
  if(!is_pregnant && next_period < days_since_last_period) {
    log_like <- -Inf
  }
  log_like
}

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, prior_is_fertile))
  prior$is_pregnant <- rbinom(n, 1, prior_is_pregnant * prior$is_fertile)
  prior$next_period <- rnorm(n, prior$mean_period, prior$sd_period)
  prior$next_period[prior$is_pregnant == 1] <- NA
  prior
}

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
}

posterior <- sample_from_posterior(days_since_last_period, days_between_periods, n_samples = 100000)

# Probability the couple is fertile
mean(posterior$is_fertile)
# Probability of pregnancy this period cycle
mean(posterior$is_pregnant)