## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(xactonomial)

## -----------------------------------------------------------------------------
true_theta <- c(.45, .15, .3, .1, .05, .15, .4, .4)

sample_data <- function(n) {
  
  T1 <- rmultinom(1, n[1], prob = true_theta[1:4])
  T2 <- rmultinom(1, n[2], prob = true_theta[5:8])
  
  list(T1 = c(T1), T2 = c(T2))
  
}

psi_bc <- function(theta) {
  
  theta1 <- theta[1:4]
  theta2 <- theta[5:8]
  sum(sqrt(theta1 * theta2))
  
}

psi_bc(true_theta)

## -----------------------------------------------------------------------------

psi_bc_v <- function(theta) {
  
  theta1 <- theta[,1:4, drop = FALSE]
  theta2 <- theta[,5:8, drop = FALSE]
  rowSums(sqrt(theta1 * theta2))
  
}

## -----------------------------------------------------------------------------
data <- sample_data(n = c(10, 12))
data
results <- xactonomial(data, psi_bc_v, psi_limits = c(0,1), conf_int = TRUE, psi0 = .5,
                          maxit = 50, chunksize = 100, ga_restart_every = 10, itp_eps = .01)
results

## -----------------------------------------------------------------------------
psi_causal <- function(pp) {
  
 max(pp)
  
}

true_psi_causal <- psi_causal(c(.4, .4, .2))

sample_data2 <- function(n) {
  
  list(rmultinom(1, n, prob = c(.4, .4, .2))[, 1])
  
}

set.seed(421)
data <- sample_data2(n = c(60))
data
results <- xactonomial(data, psi_causal, psi_limits = c(1/3,1), psi0 = .55,
                          maxit = 100, chunksize = 1000, itp_eps = .01)

results

## -----------------------------------------------------------------------------
p.ll <- xactonomial(data, psi_causal, psi_limits = c(1/3,1), psi0 = 1/3, conf_int = FALSE,
   itp_eps = .01, theta_null_points = t(rep(1/3, 3)), 
   alternative = "greater")$p.value

xactonomial(data, psi_causal, psi_limits = c(1/3,1), psi0 = 1, conf_int = FALSE,
   itp_eps = .01, theta_null_points = rbind(c(1, 0, 0), c(0,1,0), c(0,0,1)), 
   alternative = "less")$p.value


## -----------------------------------------------------------------------------
xactonomial(data, psi_causal, psi_limits = c(1/3,1), p_value_limits = c(p.ll, 1e-8),
                          maxit = 100, chunksize = 1000, itp_eps = .01)