library(cubature)
demo < - function() {
write.cmd <- function(s) { cat(sprintf('> %s\n', s)); flush.console() }
readln < - function()readline('Hit [ENTER] to continue...\n')
write.cmd('print.config()')
print.config()
readln()
write.cmd('plot.p0.posterior.density(k.data, n.data)')
plot.p0.posterior.density(k.data, n.data)
readln()
write.cmd('p0.posterior.probability(0.2, 0.4, k.data, n.data)')
print(p0.posterior.probability(0.2, 0.4, k.data, n.data))
readln()
write.cmd('plot.p1.posterior.density(k.data, n.data)')
plot.p1.posterior.density(k.data, n.data)
readln()
write.cmd('p1.posterior.probability(0.2, 0.4, k.data, n.data)')
print(p1.posterior.probability(0.2, 0.4, k.data, n.data))
readln()
write.cmd('plot.delta1.posterior.density(k.data, n.data)')
plot.delta1.posterior.density(k.data, n.data)
readln()
write.cmd('delta1.posterior.probability(0, 1, k.data, n.data)')
print(delta1.posterior.probability(0, 1, k.data, n.data))
readln()
write.cmd('plot.p2.posterior.density(k.data, n.data)')
plot.p2.posterior.density(k.data, n.data)
readln()
write.cmd('p2.posterior.probability(0.2, 0.4, k.data, n.data)')
print(p2.posterior.probability(0.2, 0.4, k.data, n.data))
readln()
write.cmd('plot.delta2.posterior.density(k.data, n.data)')
plot.delta2.posterior.density(k.data, n.data)
readln()
write.cmd('delta2.posterior.probability(0, 1, k.data, n.data)')
print(delta2.posterior.probability(0, 1, k.data, n.data))
readln()
write.cmd('plot.p3.posterior.density(k.data, n.data)')
plot.p3.posterior.density(k.data, n.data)
readln()
write.cmd('p3.posterior.probability(0.2, 0.4, k.data, n.data)')
print(p3.posterior.probability(0.2, 0.4, k.data, n.data))
readln()
write.cmd('plot.delta3.posterior.density(k.data, n.data)')
plot.delta3.posterior.density(k.data, n.data)
readln()
write.cmd('delta3.posterior.probability(0, 1, k.data, n.data)')
print(delta3.posterior.probability(0, 1, k.data, n.data))
readln()
invisible(NULL)
}
######################################################################
# Configuration: control parameters, prior parameters, and data
######################################################################
# A very small positive number
eps <- 1e-12
# Relative tolerance for computing integrals
rel.tol <- 1e-3
# each wave is two weeks long
w <- 2
# drift standard deviation
s.sigma <- 0.105 * sqrt(2)
# a large number far out in the right tail of the prior for sigma,
# the drift standard deviation.
inf.sigma <- 4 * s.sigma
# standard deviation for prior on alpha0
s.alpha <- 1.6
# sample sizes and numbers of positive response for time-steps 0 to 3
n.data <- c(100, 50, 63, 45)
k.data <- c( 25, 15, 22, 18)
print.config <- function() {
cat(sprintf('eps <- %g\n', eps))
cat(sprintf('rel.tol <- %g\n', rel.tol))
cat(sprintf('w <- %d\n', w))
cat(sprintf('s.sigma <- %g\n', s.sigma))
cat(sprintf('inf.sigma <- %g\n', inf.sigma))
cat(sprintf('s.alpha <- %g\n', s.alpha))
cat(sprintf('n.data <- c(%s)\n', paste(n.data, collapse=', ')))
cat(sprintf('k.data <- c(%s)\n', paste(k.data, collapse=', ')))
}
######################################################################
# Utilities
######################################################################
logit <- function(p)log(p/(1-p))
# RETURNS: A vector of the midpoints of the N equal-length subintervals of the
# interval from a to b.
#
midpoints <- function(a, b, N) {
a + ((1:N) - 0.5)/N * (b - a)
}
# RETURNS:
# INTEGRAL lo[n] to hi[n] : ... : INTEGRAL lo[1] to hi[1]:
# f(x1, ..., xn) d x1 ... d xn
# where n is the length of lo and hi.
#
multivariate.integrate <- function(f, lo, hi) {
x <- adaptIntegrate(f, lo, hi, tol=rel.tol)
stopifnot(x$returnCode == 0 && x$error < 2 * rel.tol)
x$integral
}
# Runs command and prints out start time and total running time.
#
time.command <- function(command) {
t0 <- Sys.time()
cat(sprintf('Start: %s\n', format(t0))); flush.console()
eval(command)
t1 <- Sys.time()
cat(sprintf('Run time: %s\n', format(t1 - t0)))
}
######################################################################
# X.joint.density() functions giving joint probability density of all
# variables up to a given time step.
######################################################################
# RETURNS: density(p0, k0=k[1] | n0=n[1]).
#
p0.joint.density <- function(p0, k, n) {
# Start with check to avoids some numerical problems.
if (p0 < 0 + eps || p0 > 1 - eps)
return(0)
k0 < - k[1]; n0 <- n[1]
alpha0 <- logit(p0)
dnorm(alpha0, 0, s.alpha) / p0 / (1 - p0) *
dbinom(k0, n0, p0)
}
# RETURNS: density(sigma, p0, k0=k[1], p1, k1=k[2] | n0=n[1], n1=n[2]).
#
p1.joint.density <- function(sigma, p0, p1, k, n) {
# Start with check to avoids some numerical problems.
if (any(c(sigma, p0, p1) < 0 + eps) || any(c(p0, p1) > 1 - eps))
return(0)
k0 < - k[1]; n0 <- n[1]
k1 <- k[2]; n1 <- n[2]
alpha0 <- logit(p0)
alpha1 <- logit(p1)
2 * dnorm(sigma, 0, s.sigma) *
dnorm(alpha0, 0, s.alpha) / p0 / (1 - p0) *
dnorm(alpha1, alpha0, s.alpha) / p1 / (1 - p1) *
dbinom(k0, n0, p0) *
dbinom(k1, n1, p1)
}
# RETURNS: density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3]).
#
p2.joint.density <- function(sigma, p0, p1, p2, k, n) {
# Start with check to avoids some numerical problems.
if (any(c(sigma, p0, p1, p2) < 0 + eps) || any(c(p0, p1, p2) > 1 - eps))
return(0)
k0 < - k[1]; n0 <- n[1]
k1 <- k[2]; n1 <- n[2]
k2 <- k[3]; n2 <- n[3]
alpha0 <- logit(p0)
alpha1 <- logit(p1)
alpha2 <- logit(p2)
2 * dnorm(sigma, 0, s.sigma) *
dnorm(alpha0, 0, s.alpha) / p0 / (1 - p0) *
dnorm(alpha1, alpha0, s.alpha) / p1 / (1 - p1) *
dnorm(alpha2, alpha1, s.alpha) / p2 / (1 - p2) *
dbinom(k0, n0, p0) *
dbinom(k1, n1, p1) *
dbinom(k2, n2, p2)
}
# RETURNS:
# density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3], p3, k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4]).
#
p3.joint.density <- function(sigma, p0, p1, p2, p3, k, n) {
# Start with check to avoids some numerical problems.
ps <- c(p0, p1, p2, p3)
if (any(c(sigma, ps) < 0 + eps) || any(ps > 1 - eps))
return(0)
k0 < - k[1]; n0 <- n[1]
k1 <- k[2]; n1 <- n[2]
k2 <- k[3]; n2 <- n[3]
k3 <- k[4]; n3 <- n[4]
alpha0 <- logit(p0)
alpha1 <- logit(p1)
alpha2 <- logit(p2)
alpha3 <- logit(p3)
2 * dnorm(sigma, 0, s.sigma) *
dnorm(alpha0, 0, s.alpha) / p0 / (1 - p0) *
dnorm(alpha1, alpha0, s.alpha) / p1 / (1 - p1) *
dnorm(alpha2, alpha1, s.alpha) / p2 / (1 - p2) *
dnorm(alpha3, alpha2, s.alpha) / p3 / (1 - p3) *
dbinom(k0, n0, p0) *
dbinom(k1, n1, p1) *
dbinom(k2, n2, p2) *
dbinom(k3, n3, p3)
}
# RETURNS: density(sigma, p0, k0=k[1], delta1, k1=k[2]
# | n0=n[1], n1=n[2])
# where delta1 = p1 - p0.
#
delta1.joint.density <- function(sigma, p0, delta1, k, n) {
p1.joint.density(sigma, p0, p0 + delta1, k, n)
}
# RETURNS: density(sigma, p0, k0=k[1], p1, k1=k[2], delta2, k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3])
# where delta2 = p2 - p0.
#
delta2.joint.density <- function(sigma, p0, p1, delta2, k, n) {
p2.joint.density(sigma, p0, p1, p0 + delta2, k, n)
}
# RETURNS:
# density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3], delta3, k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4])
# where delta3 = p3 - p0.
#
delta3.joint.density <- function(sigma, p0, p1, p2, delta3, k, n) {
p3.joint.density(sigma, p0, p1, p2, p0 + delta3, k, n)
}
######################################################################
# X.data.joint.probability() functions giving joint probability of
# (1) a variable being within a given range, and
# (2) count data up to a given time step having given values.
######################################################################
# RETURNS: Pr(a < p0 < b, k0=k[1] | n0=n[1]),
# computed as
# INTEGRAL a to b: density(p0, k0=k[1] | n0=n[1]) d p0.
#
p0.data.joint.probability <- function(a, b, k, n) {
f <- function(p0) {
p0.joint.density(p0, k, n)
}
multivariate.integrate(f, a, b)
}
# RETURNS: Pr(a < p1 < b, k0=k[1], k1=k[2] | n0=n[1], n1=n[2]),
# computed as
# INTEGRAL a to b: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2] | n0=n[1], n1=n[2])
# d sigma d p0 d p1
# where inf is a large positive number.
#
p1.data.joint.probability <- function(a, b, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
p1.joint.density(sigma, p0, p1, k, n)
}
multivariate.integrate(f, c(0, 0, a), c(inf.sigma, 1, b))
}
# RETURNS: Pr(a < delta1 < b, k0=k[1], k1=k[2] | n0=n[1], n1=n[2]),
# computed as
# INTEGRAL a to b: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], delta1, k1=k[2] | n0=n[1], n1=n[2])
# d sigma d p0 d delta1
# where inf is a large positive number and delta1 = p1 - p0.
#
delta1.data.joint.probability <- function(a, b, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
delta1 <- x[3]
delta1.joint.density(sigma, p0, delta1, k, n)
}
multivariate.integrate(f, c(0, 0, a), c(inf.sigma, 1, b))
}
# RETURNS: Pr(a < p2 < b, k0=k[1], k1=k[2], k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3]),
# computed as
# INTEGRAL a to b: INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3])
# d sigma d p0 d p1 d p2
# where inf is a large positive number.
#
p2.data.joint.probability <- function(a, b, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
p2 <- x[4]
p2.joint.density(sigma, p0, p1, p2, k, n)
}
multivariate.integrate(f, c(0, 0, 0, a), c(inf.sigma, 1, 1, b))
}
# RETURNS: Pr(a < delta2 < b, k0=k[1], k1=k[2], k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3]),
# computed as
# INTEGRAL a to b: INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2], delta2, k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3])
# d sigma d p0 d p1 d delta2
# where inf is a large positive number and delta2 = p2 - p0.
#
delta2.data.joint.probability <- function(a, b, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
delta2 <- x[4]
delta2.joint.density(sigma, p0, p1, delta2, k, n)
}
multivariate.integrate(f, c(0, 0, 0, a), c(inf.sigma, 1, 1, b))
}
# RETURNS: Pr(a < p3 < b, k0=k[1], k1=k[2], k2=k[3], k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4]),
# computed as
# INTEGRAL a to b: INTEGRAL 0 to 1:
# INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3], p3, k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4])
# d sigma d p0 d p1 d p2 d p3
# where inf is a large positive number.
#
p3.data.joint.probability <- function(a, b, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
p2 <- x[4]
p3 <- x[5]
p3.joint.density(sigma, p0, p1, p2, p3, k, n)
}
multivariate.integrate(f, c(0, 0, 0, 0, a), c(inf.sigma, 1, 1, 1, b))
}
# RETURNS: Pr(a < delta3 < b, k0=k[1], k1=k[2], k2=k[3], k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4]),
# computed as
# INTEGRAL a to b: INTEGRAL 0 to 1:
# INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3], delta3, k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4])
# d sigma d p0 d p1 d p2 d delta3
# where inf is a large positive number and delta3 = p3 - p0.
#
delta3.data.joint.probability <- function(a, b, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
p2 <- x[4]
delta3 <- x[5]
delta3.joint.density(sigma, p0, p1, p2, delta3, k, n)
}
multivariate.integrate(f, c(0, 0, 0, 0, a), c(inf.sigma, 1, 1, 1, b))
}
######################################################################
# X.data.joint.density() functions giving joint probability density
# of a variable and the count data up to a given time step.
######################################################################
# RETURNS: density(p1, k0=k[1], k1=k[2] | n0=n[1], n1=n[2]),
# computed as
# INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2] | n0=n[1], n1=n[2])
# d sigma d p0
#
p1.data.joint.density <- function(p1, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1.joint.density(sigma, p0, p1, k, n)
}
multivariate.integrate(f, c(0, 0), c(inf.sigma, 1))
}
# RETURNS: density(delta1, k0=k[1], k1=k[2] | n0=n[1], n1=n[2]),
# computed as
# INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], delta1, k1=k[2] | n0=n[1], n1=n[2])
# d sigma d p0
#
delta1.data.joint.density <- function(delta1, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
delta1.joint.density(sigma, p0, delta1, k, n)
}
multivariate.integrate(f, c(0, 0), c(inf.sigma, 1))
}
# RETURNS: density(p2, k0=k[1], k1=k[2], k2=k[3] | n0=n[1], n1=n[2], n2=n[3]),
# computed as
# INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3])
# d sigma d p0 d p1
#
p2.data.joint.density <- function(p2, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
p2.joint.density(sigma, p0, p1, p2, k, n)
}
multivariate.integrate(f, c(0, 0, 0), c(inf.sigma, 1, 1))
}
# RETURNS: density(delta2, k0=k[1], k1=k[2], k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3]),
# computed as
# INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2], delta2, k2=k[3]
# | n0=n[1], n1=n[2], n2=n[3])
# d sigma d p0 d p1
#
delta2.data.joint.density <- function(delta2, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
delta2.joint.density(sigma, p0, p1, delta2, k, n)
}
multivariate.integrate(f, c(0, 0, 0), c(inf.sigma, 1, 1))
}
# RETURNS: density(p3, k0=k[1], k1=k[2], k2=k[3], k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4]),
# computed as
# INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3], p3, k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4])
# d sigma d p0 d p1 d p2
#
p3.data.joint.density <- function(p3, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
p2 <- x[4]
p3.joint.density(sigma, p0, p1, p2, p3, k, n)
}
multivariate.integrate(f, c(0, 0, 0, 0), c(inf.sigma, 1, 1, 1))
}
# RETURNS: density(delta3, k0=k[1], k1=k[2], k2=k[3], k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4]),
# computed as
# INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to 1: INTEGRAL 0 to inf:
# density(sigma, p0, k0=k[1], p1, k1=k[2], p2, k2=k[3], delta3, k3=k[4]
# | n0=n[1], n1=n[2], n2=n[3], n3=n[4])
# d sigma d p0 d p1 d p2
#
delta3.data.joint.density <- function(delta3, k, n) {
f <- function(x) {
sigma <- x[1]
p0 <- x[2]
p1 <- x[3]
p2 <- x[4]
delta3.joint.density(sigma, p0, p1, p2, delta3, k, n)
}
multivariate.integrate(f, c(0, 0, 0, 0), c(inf.sigma, 1, 1, 1))
}
######################################################################
# X.posterior.density() functions.
######################################################################
# NOTE: vapply(vec, func, 0) creates the vector
# c(func(vec[1]), func(vec[2]), ..., func(vec[n])
# where n is the length of vec.
# RETURNS: A vector of the values
# density(p0 | k0=k[1], n0=n[1])
# for all values p0 in p0.vec.
#
p0.posterior.density <- function(p0.vec, k, n) {
f <- function(p0)p0.joint.density(p0, k, n)
numerator <- vapply(p0.vec, f, 0)
denominator <- p0.data.joint.probability(0, 1, k, n)
numerator / denominator
}
# RETURNS: A vector of the values
# density(p1 | k0=k[1], n0=n[1], k1=k[2], n1=n[2])
# for all values p1 in p1.vec
#
p1.posterior.density <- function(p1.vec, k, n) {
f <- function(p1)p1.data.joint.density(p1, k, n)
numerator <- vapply(p1.vec, f, 0)
denominator <- p1.data.joint.probability(0, 1, k, n)
numerator / denominator
}
# RETURNS: A vector of the values
# density(delta1 | k0=k[1], n0=n[1], k1=k[2], n1=n[2])
# for all values delta1 in delta1.vec
#
delta1.posterior.density <- function(delta1.vec, k, n) {
f <- function(delta1)delta1.data.joint.density(delta1, k, n)
numerator <- vapply(delta1.vec, f, 0)
denominator <- delta1.data.joint.probability(0, 1, k, n)
numerator / denominator
}
# RETURNS: A vector of the values
# density(p2 | k0=k[1], n0=n[1], k1=k[2], n1=n[2], k2=k[3], n2=n[3])
# for all values p2 in p2.vec
#
p2.posterior.density <- function(p2.vec, k, n) {
f <- function(p2)p2.data.joint.density(p2, k, n)
numerator <- vapply(p2.vec, f, 0)
denominator <- p2.data.joint.probability(0, 1, k, n)
numerator / denominator
}
# RETURNS: A vector of the values
# density(delta2 | k0=k[1], n0=n[1], k1=k[2], n1=n[2], k2=k[3], n2=n[3])
# for all values delta2 in delta2.vec
#
delta2.posterior.density <- function(delta2.vec, k, n) {
f <- function(delta2)delta2.data.joint.density(delta2, k, n)
numerator <- vapply(delta2.vec, f, 0)
denominator <- delta2.data.joint.probability(0, 1, k, n)
numerator / denominator
}
# RETURNS: A vector of the values
# density(p3 | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3], k3=k[4], n3=n[4])
# for all values p3 in p3.vec
#
p3.posterior.density <- function(p3.vec, k, n) {
f <- function(p3)p3.data.joint.density(p3, k, n)
numerator <- vapply(p3.vec, f, 0)
denominator <- p3.data.joint.probability(0, 1, k, n)
numerator / denominator
}
# RETURNS: A vector of the values
# density(delta3 | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3], k3=k[4], n3=n[4])
# for all values delta3 in delta3.vec
#
delta3.posterior.density <- function(delta3.vec, k, n) {
f <- function(delta3)delta3.data.joint.density(delta3, k, n)
numerator <- vapply(delta3.vec, f, 0)
denominator <- delta3.data.joint.probability(0, 1, k, n)
numerator / denominator
}
######################################################################
# X.posterior.probability() functions giving the posterior
# probability that a variable is in a given range.
######################################################################
# RETURNS: Pr(a < p0 < b | k0=k[1], n0=n[1]).
#
p0.posterior.probability <- function(a, b, k, n) {
p0.data.joint.probability(a, b, k, n) /
p0.data.joint.probability(0, 1, k, n)
}
# RETURNS: Pr(a < p1 < b | k0=k[1], n0=n[1], k1=k[2], n1=n[2]).
#
p1.posterior.probability <- function(a, b, k, n) {
p1.data.joint.probability(a, b, k, n) /
p1.data.joint.probability(0, 1, k, n)
}
# RETURNS: Pr(a < delta1 < b | k0=k[1], n0=n[1], k1=k[2], n1=n[2]).
#
delta1.posterior.probability <- function(a, b, k, n) {
delta1.data.joint.probability(a, b, k, n) /
delta1.data.joint.probability(-1, 1, k, n)
}
# RETURNS: Pr(a < p2 < b | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3]).
#
p2.posterior.probability <- function(a, b, k, n) {
p2.data.joint.probability(a, b, k, n) /
p2.data.joint.probability(0, 1, k, n)
}
# RETURNS: Pr(a < delta2 < b | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3]).
#
delta2.posterior.probability <- function(a, b, k, n) {
delta2.data.joint.probability(a, b, k, n) /
delta2.data.joint.probability(-1, 1, k, n)
}
# RETURNS: Pr(a < p3 < b | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3], k3=k[4], n3=n[4]).
#
p3.posterior.probability <- function(a, b, k, n) {
p3.data.joint.probability(a, b, k, n) /
p3.data.joint.probability(0, 1, k, n)
}
# RETURNS: Pr(a < delta3 < b | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3], k3=k[4], n3=n[4]).
#
delta3.posterior.probability <- function(a, b, k, n) {
delta3.data.joint.probability(a, b, k, n) /
delta3.data.joint.probability(-1, 1, k, n)
}
######################################################################
# plot.X.posterior.density() functions.
######################################################################
# Plots density(p0 | k0=k[1], n0=n[1]).
#
plot.p0.posterior.density <- function(k, n) {
pvec <- midpoints(0, 1, 500)
time.command(
plot(pvec, p0.posterior.density(pvec, k, n),
type='l', ylab='posterior density', xlab='p0')
)
}
# Plots density(p1 | k0=k[1], n0=n[1], k1=k[2], n1=n[2])
#
plot.p1.posterior.density <- function(k, n) {
pvec <- midpoints(0, 1, 500)
time.command(
plot(pvec, p1.posterior.density(pvec, k, n),
type='l', ylab='posterior density', xlab='p1')
)
}
# Plots density(delta1 | k0=k[1], n0=n[1], k1=k[2], n1=n[2])
#
plot.delta1.posterior.density <- function(k, n) {
delta.vec <- midpoints(-.5, .5, 500)
time.command(
plot(delta.vec, delta1.posterior.density(delta.vec, k, n),
type='l', ylab='posterior density', xlab='delta1')
)
}
# Plots density(p2 | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3])
#
plot.p2.posterior.density <- function(k, n) {
pvec <- midpoints(0, 1, 300)
time.command(
plot(pvec, p2.posterior.density(pvec, k, n),
type='l', ylab='posterior density', xlab='p2')
)
}
# Plots density(delta2 | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3])
#
plot.delta2.posterior.density <- function(k, n) {
delta.vec <- midpoints(-0.5, 0.5, 300)
time.command(
plot(delta.vec, delta2.posterior.density(delta.vec, k, n),
type='l', ylab='posterior density', xlab='delta2')
)
}
# Plots density(p3 | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3], k3=k[4], n3=n[4])
#
plot.p3.posterior.density <- function(k, n) {
pvec <- midpoints(0, 1, 300)
time.command(
plot(pvec, p3.posterior.density(pvec, k, n),
type='l', ylab='posterior density', xlab='p3')
)
}
# Plots density(delta3 | k0=k[1], n0=n[1], k1=k[2], n1=n[2],
# k2=k[3], n2=n[3], k3=k[4], n3=n[4])
#
plot.delta3.posterior.density <- function(k, n) {
delta.vec <- midpoints(-0.5, 0.5, 300)
time.command(
plot(delta.vec, delta3.posterior.density(delta.vec, k, n),
type='l', ylab='posterior density', xlab='delta3')
)
}
