Calculate multinomial probabilities

Description

Calculate multinomial probabilities

Usage

calc_multinom_probs(sar, logt, logc, d, n, nt)

Arguments

Value

A vector of probabilities

Examples

sspace_3_5 <- sspace_multinom(3, 5)
calc_multinom_probs(sspace_3_5, sample_unit_simplexn(3, 10),
  apply(matrix(sspace_3_5, ncol = 3, byrow = TRUE), 1, log_multinom_coef, sumx = 5), 3, 5, 10)

Calculate probability for given parameters

Description

Given a set of candidate parameter vectors, the enumerated sample space, and a logical vector with the same number of elements of the sample space, compute the probability for each element of the sample space and take the sum.

Usage

calc_prob_null(theta_cands, SSpacearr, logC, II)

Arguments

Value

A numeric vector of probabilities

Examples

sspace_3_5 <- matrix(sspace_multinom(3, 5), ncol = 3, byrow = TRUE)
theta_cands <- matrix(sample_unit_simplexn(3, 10), ncol = 3,byrow = TRUE)
calc_prob_null_fast(theta_cands, sspace_3_5,
apply(sspace_3_5, 1, log_multinom_coef, sumx = 5), II = 1:21 > 12)
# same as below but faster
calc_prob_null(theta_cands, sspace_3_5,
apply(sspace_3_5, 1, log_multinom_coef, sumx = 5), II = 1:21 > 12)

Calculate probability for given parameters

Description

Given a set of candidate parameter vectors, the enumerated sample space, and a logical vector with the same number of elements of the sample space, compute the probability for each element of the sample space and take the sum.

Usage

calc_prob_null_fast(theta_cands, SSpacearr, logC, II)

Arguments

Value

A numeric vector of probabilities

Examples

sspace_3_5 <- matrix(sspace_multinom(3, 5), ncol = 3, byrow = TRUE)
theta_cands <- matrix(sample_unit_simplexn(3, 10), ncol = 3,byrow = TRUE)
calc_prob_null_fast(theta_cands, sspace_3_5,
apply(sspace_3_5, 1, log_multinom_coef, sumx = 5), II = 1:21 > 12)
# same as below but faster
calc_prob_null(theta_cands, sspace_3_5,
apply(sspace_3_5, 1, log_multinom_coef, sumx = 5), II = 1:21 > 12)

Gradient of the multinomial likelihood sum

Description

Gradient of the multinomial likelihood sum

Usage

calc_prob_null_gradient(theta_cands, SSpacearr, II)

Arguments

Value

A matrix the same dimension as theta_cands

Examples

calc_prob_null_gradient(t(c(.28, .32, .4)),
matrix(c(2, 2, 1, 1, 2, 2, 0, 3, 2), ncol = 3),
rep(TRUE, 3))

# numerically
testenv <- new.env()
testenv$SSpacearr <- matrix(c(2, 2, 1, 1, 2, 2, 0, 3, 2), ncol = 3)
testenv$thistheta <- c(.28, .32, .4)
numericDeriv(quote(sum(exp((.colSums(t(SSpacearr) * log(thistheta), m = 3, n = 3))))),
    theta = "thistheta", rho = testenv, central = TRUE)

Arrange all combinations of rows of two matrices

Description

Given X and Y, both matrices where the rows are counts of multinomial trials, produce all combinations rowwise, concatenate the rows into a new matrix, and calculate the log multinomial coefficients for the combination.

Usage

combinate(X, Y)

Arguments

Value

A list containing Sspace, the sample space (vectors of counts), and logC, a vector of the log multinomial coefficients.

Examples

slist_2_3 <- combinate(matrix(sspace_multinom(2, 5), ncol = 2, byrow = TRUE),
   matrix(sspace_multinom(3, 6), ncol = 3, byrow = TRUE))

Like combinate but adds on to previous call

Description

Like combinate but adds on to previous call

Usage

combinate2(X, Y)

Arguments

Value

A list containing Sspace, the sample space (vectors of counts), and logC, a vector of the log multinomial coefficients.

Examples

slist_2_3 <- combinate(matrix(sspace_multinom(2, 5), ncol = 2, byrow = TRUE),
   matrix(sspace_multinom(3, 6), ncol = 3, byrow = TRUE))

sl_2_3_4 <- combinate2(slist_2_3, matrix(sspace_multinom(4, 3), ncol = 4, byrow = TRUE))

Find a univariate root of the function f

Description

This finds the value x \in [a, b] such that f(x) = 0 using the one-dimensional root finding ITP method (Interpolate Truncate Project). Also see itp.

Usage

itp_root(
  f,
  a,
  b,
  k1 = 0.1,
  k2 = 2,
  n0 = 1,
  eps = 0.005,
  maxiter = 100,
  fa = NULL,
  fb = NULL,
  verbose = FALSE,
  ...
)

Arguments

Value

A numeric vector of length 1, the root at the last iteration

References

I. F. D. Oliveira and R. H. C. Takahashi. 2020. An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality. ACM Trans. Math. Softw. 47, 1, Article 5 (March 2021), 24 pages. https://doi.org/10.1145/3423597

Examples

fpoly <- function(x) x^3 - x - 2 ## example from the ITP_method wikipedia entry
itp_root(fpoly, 1, 2, eps = .0001, verbose = TRUE)

Calculate log of multinomial coefficient

Description

Calculate log of multinomial coefficient

Usage

log_multinom_coef(x, sumx)

Arguments

Value

The vector of log multinomial coefficients

Examples

S0 <- matrix(sspace_multinom(4, 6), ncol = 4, byrow = TRUE)
logC0<- apply(S0,1,log_multinom_coef,sumx=6)

Compute a p value for the test of psi <= psi0 (lower = TRUE) or psi >= psi0 (lower = FALSE)

Description

Compute a p value for the test of psi <= psi0 (lower = TRUE) or psi >= psi0 (lower = FALSE)

Usage

pvalue_psi0(
  psi0,
  psi,
  psi_hat,
  psi_obs,
  alternative = "two.sided",
  maxit,
  chunksize,
  p_target,
  SSpacearr,
  logC,
  d_k,
  psi_is_vectorized = FALSE,
  theta_sampler = runif_dk_vects,
  ga = FALSE,
  ga_gfactor = 1,
  ga_lrate = 0.01,
  ga_restart_every = 10,
  warn = TRUE
)

Arguments

Value

A vector with two p-values, one for the lower, and one for the greater

Examples

sspace_3_5 <- matrix(sspace_multinom(3, 5), ncol = 3, byrow = TRUE)
psi <- function(theta) max(theta)
logC <- apply(sspace_3_5, 1, log_multinom_coef, sumx = 5)
psi_hat <- apply(sspace_3_5, 1, \(x) psi(x / sum(x)))
pvalue_psi0(.3, psi, psi_hat, .4, maxit = 10, chunksize = 100,
 p_target = 1, SSpacearr = sspace_3_5, logC = logC, d_k = 3, warn = FALSE)

Sample independently from Dirichlet distributions for each of d_k vectors

Description

Sample independently from Dirichlet distributions for each of d_k vectors

Usage

rdirich_dk_vects(nsamp, alpha)

Arguments

Value

A matrix with sum(d_k) columns and nsamp rows

Examples

rdirich_dk_vects(10, list(rep(1, 3), rep(1, 4), rep(1, 2)))

Sample uniformly and independently from d_k simplices

Description

Sample uniformly and independently from d_k simplices

Usage

runif_dk_vects(d_k, nsamp)

Arguments

Value

A matrix with sum(d_k) columns and nsamp rows

Examples

runif_dk_vects(c(3, 4, 2), 10)

Sample n times from the unit simplex in d dimensions

Description

Sample n times from the unit simplex in d dimensions

Usage

sample_unit_simplexn(d, n)

Arguments

Value

The grid over Theta, the parameter space. To be converted to a matrix with d columns and nsamp rows

Examples

matrix(sample_unit_simplexn(3, 10), ncol = 3, byrow = TRUE)

Enumerate the multinomial sample space

Description

Enumerate the multinomial sample space

Usage

sspace_multinom(d, n)

Arguments

Value

A vector enumerating the sample space, to be converted to a matrix with d columns and choose(n + d - 1, d - 1) rows

Examples

matrix(sspace_multinom(3, 5), ncol = 3, byrow = TRUE)

Enumerate the sample space of a multinomial

Description

We have d mutually exclusive outcomes and n independent trials. This function enumerates all possible vectors of length d of counts of each outcome for n trials, i.e., the sample space. The result is output as a matrix with d columns where each row represents a possible observation. See sspace_multinom for a faster implementation using Rust.

Usage

sspace_multinom_slow(d, n)

Arguments

Value

A matrix with d columns

Examples

d4s <- sspace_multinom_slow(4, 8)
stopifnot(abs(sum(apply(d4s, 1, dmultinom, prob = rep(.25, 4))) - 1) < 1e-12)

Improved inference for a real-valued function of multinomial parameters

Description

We consider the k sample multinomial problem where we observe k vectors (possibly of different lengths), each representing an independent sample from a multinomial. For a given function psi which takes in the concatenated vector of multinomial probabilities and outputs a real number, we are interested in computing a p-value for a test of psi >= psi0, and constructing a confidence interval for psi.

Usage

xactonomial(
  data,
  psi,
  statistic = NULL,
  psi0 = NULL,
  alternative = c("two.sided", "less", "greater"),
  psi_limits,
  theta_null_points = NULL,
  p_target = 1,
  conf_int = TRUE,
  conf_level = 0.95,
  itp_maxit = 10,
  itp_eps = 0.005,
  p_value_limits = NULL,
  maxit = 50,
  chunksize = 500,
  theta_sampler = runif_dk_vects,
  ga = TRUE,
  ga_gfactor = "adapt",
  ga_lrate = 0.01,
  ga_restart_every = 10
)

Arguments

Details

Let T_j be distributed \mbox{Multinomial}_{d_j}(\boldsymbol{\theta}_j, n_j) for j = 1, \ldots, k and denote \boldsymbol{T} = (T_1, \ldots, T_k) and \boldsymbol{\theta} = (\theta_1, \ldots, \theta_k). The subscript d_j denotes the dimension of the multinomial. Suppose one is interested in the parameter \psi = \tau(\boldsymbol{\theta}) \in \Psi \subseteq \mathbb{R}. Given a sample of size n from \boldsymbol{T}, say \boldsymbol{X} = (X_1, \ldots, X_k), which is a vector of counts obtained by concatenating the k independent count vectors, let G(\boldsymbol{X}) denote a real-valued statistic that defines the ordering of the sample space. Tne default choice of the statistic is to estimate \boldsymbol{\theta} with the sample proportions and plug them into \tau(\boldsymbol{\theta}). This function calculates a p value for a test of the null hypothesis H_0: \psi(\boldsymbol{\theta}) \neq \psi_0 for the two sided case, H_0: \psi(\boldsymbol{\theta}) \leq \psi_0 for the case alternative = "greater", and H_0: \psi(\boldsymbol{\theta}) \geq \psi_0 for the case alternative = "less". We make no assumptions and do not rely on large sample approximations. It also optionally constructs a 1 - \alpha percent confidence interval for \psi. The computation is somewhat involved so it is best for small sample sizes. The calculation is done by sampling a large number of points from the null parameter space \Theta_0, then computing multinomial probabilities under those values for the range of the sample space where the statistic is as or more extreme than the observed statistic given data. It is basically the definition of a p-value implemented with Monte Carlo methods. Some options for speeding up the calculation are available.

Value

An object of class "htest", which is a list with the following elements:

estimate

The value of the statistic at the observed data

p.value

The p value

conf.int

The upper and lower confidence limits

null.value

The null hypothesis value provided by the user

alternative

The type of test

method

A description of the method

data.name

The name of the data object provided by the user

p.sequence

A list with two elements, p.null and p.alt containing the vector of p values at each iteration for the less than null and the greater than null. Used for assessing convergence.

Specifying the function psi

The psi parameter should be a function that either: 1) takes a vector of length sum(d_j) (the total number of bins) and outputs a single number, or 2) takes a matrix with number of columns equal to sum(d_j), and arbitrary number of rows and outputs a vector with length equal to the number of rows. In other words, psi can be not vectorized or it can be vectorized by row. Writing it so that it is vectorized can speed up the calculation. See examples.

Boundary issues

It is required to provide psi_limits, a vector of length 2 giving the smallest and largest possible values that the function psi can take, e.g., c(0, 1). If the null hypothesis value psi0 is at one of the limits, it is often the case that sampling from the null parameter space is impossible because it is a set of measure 0. While it may have measure 0, it is not empty, and will contain a finite set of points. Thus you should provide the argument theta_null_points which is a matrix where the rows contain the finite set (sometimes 1) of points \theta such that \tau(\theta) = \psi_0. There is also an argument called p_value_limits that can be used to improve performance of confidence intervals around the boundary. This should be a vector of length 2 with the p-value for a test of psi_0 <= psi_limits[1] and the p-value for a test of psi_0 >= psi_limits[2]. See examples.

Optimization options

For p-value calculation, you can provide a parameter p_target, so that the sampling algorithm terminates when a p-value is found that exceeds p_target. The algorithm begins by sampling uniformly from the unit simplices defining the parameter space, but alternatives can be specified in theta_sampler. By default gradient ascent (ga = TRUE) is performed during the p-value maximization procedure, and ga_gfactor and ga_lrate control options for the gradient ascent. At each iteration, the gradient of the multinomial probability at the current maximum theta is computed, and a step is taken to theta + lrate * gradient. Then for the next iteration, a set of chunksize samples are drawn from a Dirichlet distribution with parameter ga_gfactor * (theta + ga_lrate * gradient). If ga_gfactor = "adapt" then it is set to 1 / max(theta) at each iteration. The ITP algorithm itp_root is used to find roots of the p-value function as a function of the psi0 value to get confidence intervals. The maximum number of iterations and epsilon can be controlled via itp_maxit, itp_eps.

References

Sachs, M.C., Gabriel, E.E. and Fay, M.P., 2024. Exact confidence intervals for functions of parameters in the k-sample multinomial problem. arXiv preprint arXiv:2406.19141.

Examples

psi_ba <- function(theta) {
  theta1 <- theta[1:4]
  theta2 <- theta[5:8]
  sum(sqrt(theta1 * theta2))
  }
data <- list(T1 = c(2,1,2,1), T2 = c(0,1,3,3))
xactonomial(data, psi_ba, psi_limits = c(0, 1), psi0 = .5,
  conf_int = FALSE, maxit = 15, chunksize = 200)

# vectorized by row
psi_ba_v <- function(theta) {
theta1 <- theta[,1:4, drop = FALSE]
theta2 <- theta[,5:8, drop = FALSE]
rowSums(sqrt(theta1 * theta2))
}
data <- list(T1 = c(2,1,2,1), T2 = c(0,1,3,3))
xactonomial(data, psi_ba_v, psi_limits = c(0, 1), psi0 = .5,
 conf_int = FALSE, maxit = 10, chunksize = 200)

 # example of using theta_null_points
 # psi = 1/3 occurs when all probs = 1/3
 psi_max <- function(pp) {
   max(pp)
 }

data <- list(c(13, 24, 13))

xactonomial(data, psi_max, psi_limits = c(1 / 3, 1), psi0 = 1/ 3,
  conf_int = FALSE, theta_null_points = t(c(1/3, 1/3, 1/3)))

## in this case using p_value_limits improves confidence interval performance

 xactonomial(data, psi_max, psi_limits = c(1 / 3, 1), psi0 = 1/ 3,
  conf_int = TRUE, theta_null_points = t(c(1/3, 1/3, 1/3)),
  p_value_limits = c(.1, 1e-8))