Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
A script that implements a Bayesian model calculating the probability that a couple is fertile and is going to be pregnant.
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment