binomial-llm-comp-examp.R

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')
  )
}

Leave a Reply

Your email address will not be published. Required fields are marked *