Confidence boxes - Consonant c-boxes

Confidence structures fuse the comprehensiveness and flexibility of Bayesian inference with the statistical performance and rigor of classical frequentist inference. Consonance is a consistency property that guarantees that increasing the confidence level (alpha) will yield confidence intervals that enclose those at all lower confidence levels. Drawing on ideas and methods from possibility theory, Dempster–Shafer theory, and random sets theory, consonant confidence structures emerge as a scheme for estimation appropriate for large and small data sets alike. They are especially useful in engineering applications because they are allow robust estimation with in-built statistical performance guarantees that do not require commitment to a statistical prior and are free from the false confidence problems and other paradoxes attending all additive Bayesian and frequent approaches. Additionally, they can be visualised in simple graphical depictions with obvious interpretations clear even to lay readers.

Shown below are six estimates for the probability that generated 2 successes out of 5 independent trials, i.e., k = 2 and n = 5.

The bottom, left plot shows the Bayesian estimate for this probability which is the beta distribution beta(½+k, ½+n-k) which arises from the Jeffreys prior. The bottom, right plot shows the Bayesian estimate beta(1+k, 1+n-k) which arises under the Bayes-Laplace uniform prior. These two Bayesian estimates, both shown in black, are similar, although other priors can yield still different estimates.

The middle, right estimate is Walley's imprecise beta model is a collection of beta distributions beta([k, s+k], [n-k, s+n-k]), with s=1. The middle, left estimate is the c-box env(beta(k, n-k+1), beta(k+1, n-k)), where env denotes the envelope, described in https://sites.google.com/site/confidenceboxes/isipta. These imprecise estimates, both shown in blue, are identical in form if not interepretation. The Walley estimate however will vary with different values of s.

The upper, left estimate is a consonant (but not unimodal) c-box computed using the brute-force algorithm https://sites.google.com/site/cboxbinomial/. The upper, right estimate is a consonant and unimodal confidence structure defined in Michael Balch's new paper (<<ref>>), computed with the algorithm on this page below. Note that this structure is considerably tighter than the brute-force structure to the left.

################################################################################

# see Michael's file binomialCCS.r

# Michael Scott Balch, PhD

# Alexandria Validation Consulting, LLC

# 2018-02-10

# Algorithms for binomial inference.

#

# binomial_pvalue - function for computing two-sided p-values

# binomial_ConfInt - function for computing two-sided confidence intervals

#

# These functions are based on formulas derived in Section 5 of

# Balch, M. (2019) "A Unimodal Consonant Confidence Structures [sic] for

# Binomial Inference Derived Using the P-Value Formulation."

#

# Note that these functions are designed for binomial inference

# in which the data were generated from a fixed number of trials

# experiment. The results may not be appropriate for an experiment

# with different stopping rules.

#

# FURTHER note that, while these algorithms do deliver p-values to within

# machine error and confidence intervals to within user-specified tolerance,

# it is up to the user to apply them properly.

#

# However, there is more to statistical inference than correctly executing the mathematics.

# The user takes all responsibility for the proper application of these algorithms.

print('YOU HAVE LOADED ALGORITHMS FOR BINOMIAL INFERENCE')

print('DEVELOPED AND DISTRIBUTED BY ALEXANDRIA VALIDATION CONSULTING, LLC (AVCLLC).')

print('READ DOCUMENTATION BEFORE USE.')

print('DOCUMENTATION TAKES THE FORM OF COMMENTS TO binomialCCS.R, WHICH YOU JUST LOADED.')

print('BY USING THESE ALGORITHMS, THE USER ACCEPTS ALL RESPONSIBILITY FOR THEIR APPLICATION.')

binomial_pvalue <- function( n , k , th ){

# This function returns exact two-sided p-values (i.e., plausibilities for binomial inference.

#

# INPUTS:

# n, number of trials, must be a single whole number

# k, number of "successes", must be a single non-negative integer less than or equal to n

# th, hypothetical value(s) of parameter (i.e., underlying probability)

# being assessed. Can be a singleton or array of values between zero and one

#

# OUTPUT:

# pval, pointwise plausibility at desired theta points, given k successes in fixed n trials

# Will be singleton or array with same length as input th

thMin = min(th)

thMax = max(th)

if( ( thMin < 0 ) | ( thMax > 1 ) ){

print('ERROR! ALL THETA VALUES MUST BE IN [0,1]')

}

if( (length( n ) != 1) | (length(k) != 1) ){

print('ERROR! Function will only accept a single value each, for n and k.')

}

if( (n < 1) | (n <= k) | (k < 0) | ( floor(n) != n ) | ( floor(k) != k ) ){

print('ERROR! n must be whole number. k must be integer between 0 and n (inclusive).')

}

pval = array( 0 , length(th) )

if( k > 0 ){

thLmid = critThetaWalley( n , k , k-1 )

}else{

thLmid = 0

}

if( k < n ){

thRmid = critThetaWalley( n , k , k+1 )

}else{

thRmid = 1

}

JLoc = which( ( thLmid <= th ) & ( th <= thRmid ) )

pval[JLoc] = 1

thL = thLmid

thR = thRmid

ii = k-1

while( ( ii >= 0 ) & ( thL > thMin ) ){ # working backwards from middle

thR = thL

thL = critThetaWalley( n , k , ii-1 )

JLoc = which( ( thL <= th ) & ( th < thR ) )

pval[JLoc] = pbinom( ii-1 , n , th[JLoc] ) + 1 - pbinom( k-1 , n , th[JLoc] )

ii = ii-1

}

ii=k+1

thR = thRmid

thL = thLmid

while( ( ii <= n ) & ( thR < thMax ) ){ # working forwards from middle

thL = thR

thR = critThetaWalley( n , k , ii+1 )

JLoc = which( ( thL < th ) & ( thR >= th ) )

pval[JLoc] = 1 - pbinom( ii , n , th[JLoc] ) + pbinom( k , n , th[JLoc] )

ii = ii+1

}

return( pval )

}

binomial_ConfInt <- function( n , k , a ){

# This function returns two-sided confidence intervals

# for binomial inference.

#

# INPUTS:

# n, number of trials, must be a single whole number

# k, number of "successes", must be a single non-negative integer less than or equal to n

# a, 1-Confidence level of desired confidence intervals

# Can be a singleton or array of values between zero and one

#

# OUTPUT:

# ConfBounds, 1-a confidence interval(s)

# a 2D array with 2 columns and with each row corresponding to a value of input "a"

# ConfBounds[,1] are the lower bounds

# ConfBounds[,2] are the upper bounds

tol = a * 1e-6

JLoc = which( tol < 1e-12 )

tol[JLoc] = 1e-12

LB = array( 0 , length(a) )

UB = array( 1 , length(a) )

JLoc = which( a == 0 )

LB[JLoc] = 0

UB[JLoc] = 1

Jelse = setdiff( 1:length(a) , JLoc )

if( length(Jelse) > 0 ){

aMin = min( a[Jelse] )

}else{

aMin = 0

}

p1 = 1

if( k == 0 ){

thL = 0

}else{

thL = critThetaWalley( n , k , k-1 )

}

ii = k-1

while( ( aMin > 0 ) & ( ii >= 0 ) & ( p1 >= aMin ) ){ # Moving Left from Middle

thR = thL

thL = critThetaWalley( n , k , ii-1 )

p3 = p1

p2 = pbinom( ii-1 , n , thR ) + 1 - pbinom( k-1 , n , thR )

p1 = pbinom( ii-1 , n , thL ) + 1 - pbinom( k-1 , n , thL )

# Part One: Discrete Jump

JLoc = which( ( p2 <= a ) & ( a <= p3 ) )

LB[JLoc] = thR

# Part Two: Continuous Pointwise Plausibility Region

JLoc = which( ( p1 < a ) & ( a < p2 ) )

thLsec = array( thL , length(JLoc) )

thRsec = array( thR , length(JLoc) )

errLsec = p1 - a[JLoc]

errRsec = p2 - a[JLoc]

while( length( JLoc ) > 0 ){ # Secant method for a-values in the continuous region

thNew = ( errRsec*thLsec - errLsec*thRsec ) / ( errRsec - errLsec )

errNew = pbinom( ii-1 , n , thNew ) + 1 - pbinom( k-1 , n , thNew ) - a[JLoc]

JFin = which( abs( errNew ) <= tol[JLoc] )

LB[ JLoc[JFin] ] = thNew[JFin] # THIS IS WHERE FINISHED THETA GETS WRITTEN TO BOUND

JL = which( ( ( errNew * errLsec ) > 0 ) & ( abs(errNew) > tol[JLoc] ) )

JR = which( ( ( errNew * errRsec ) > 0 ) & ( abs(errNew) > tol[JLoc] ) )

JKeepGo = sort( c(JL,JR) )

errLsec[JL] = errNew[JL]

errRsec[JR] = errNew[JR]

thLsec[JL] = thNew[JL]

thRsec[JR] = thNew[JR]

thLsec = thLsec[JKeepGo]

thRsec = thRsec[JKeepGo]

errLsec = errLsec[JKeepGo]

errRsec = errRsec[JKeepGo]

JLoc = JLoc[JKeepGo]

}

ii = ii-1

}

p3 = 1

if( k == n ){

thR = 1

}else{

thR = critThetaWalley( n , k , k+1 )

}

ii = k+1

while( ( aMin > 0 ) & ( ii <= n ) & ( p3 >= aMin ) ){ # Moving Right from Middle

thL = thR

thR = critThetaWalley( n , k , ii+1 )

p1 = p3

p2 = 1 - pbinom( ii , n , thL ) + pbinom( k , n , thL )

p3 = 1 - pbinom( ii , n , thR ) + pbinom( k , n , thR )

# Part One: Discrete Jump

JLoc = which( ( p2 <= a ) & ( a <= p1 ) )

UB[JLoc] = thL

# Part Two: Continuous Pointwise Plausibility Region

JLoc = which( ( p3 < a ) & ( a < p2 ) )

thLsec = array( thL , length(JLoc) )

thRsec = array( thR , length(JLoc) )

errLsec = p2 - a[JLoc]

errRsec = p3 - a[JLoc]

while( length( JLoc ) > 0 ){ # Secant method for a-values in the continuous region

thNew = ( errRsec*thLsec - errLsec*thRsec ) / ( errRsec - errLsec )

errNew = 1 - pbinom( ii , n , thNew ) + pbinom( k , n , thNew ) - a[JLoc]

JFin = which( abs( errNew ) <= tol[JLoc] )

UB[ JLoc[JFin] ] = thNew[JFin] # THIS IS WHERE FINISHED THETA GETS WRITTEN TO BOUND

JL = which( ( errNew * errLsec > 0 ) & ( abs(errNew) > tol[JLoc] ) )

JR = which( ( errNew * errRsec > 0 ) & ( abs(errNew) > tol[JLoc] ) )

JKeepGo = sort( c(JL,JR) )

errLsec[JL] = errNew[JL]

errRsec[JR] = errNew[JR]

thLsec[JL] = thNew[JL]

thRsec[JR] = thNew[JR]

thLsec = thLsec[JKeepGo]

thRsec = thRsec[JKeepGo]

errLsec = errLsec[JKeepGo]

errRsec = errRsec[JKeepGo]

JLoc = JLoc[JKeepGo]

}

ii = ii+1

}

ConfBounds = array( c(LB,UB) , c( length(a) , 2 ) )

return( ConfBounds )

}

critThetaWalley <- function( n , kRaw , iRaw ){

# This function obtains the values at which the Walley IM jumps

# i.e. th_{k,i} where i < k or k < i

# designed to handle singleton kRaw

# and vector iRaw

if( length(kRaw) != 1 ){

print('ERROR! critThetaWalley ONLY DESIGNED TO HANDLE SINGLETON kRaw')

}

iVec = iRaw*0

kVec = kRaw*0

g = iRaw*0

J0 = which( iRaw < kRaw )

J1 = which( iRaw > kRaw )

J2 = which( iRaw == kRaw )

# print( J0 )

# print( J1 )

# print( J2 )

iVec[J0] = iRaw[J0]

kVec[J0] = kRaw

iVec[J1] = kRaw

kVec[J1] = iRaw[J1]

J3 = which( iVec < (kVec-1) )

J4 = which( iVec == (kVec-1) )

g=iRaw*0

k = kVec[J3]

i = iVec[J3]

g[J3] = ( lbeta( k,n-k+1 )-lbeta( i+1,n-i ) ) / ( k-i-1 )

k = kVec[J4]

i = iVec[J4]

g[J4] = digamma( k ) - digamma( n-i )

thCrit=iRaw*0

J5 = which( g>0 )

J6 = which( g<=0 )

thCrit[J5] = 1 / ( 1 + exp(-g[J5]) )

thCrit[J6] = exp( g[J6] ) / ( 1 + exp( g[J6] ) )

if( length(J2) > 0 ){

print('WARNING! i==k for at least one i-value. Returns NaN')

thCrit[J2] = 0/0

}

return( thCrit )

}

b = binomial_pvalue(n,k,(1:10)/10)

b # [1] 0.08146 0.59040 0.63985 1.00000 1.00000 0.65440 0.33115 0.05792 0.00856 0.00000

b = binomial_pvalue(n,k,(1:100)/100)

plot((1:100)/100, b, type='l')

bp = binomial_ConfInt( n , k , 1-p )

lines(bp,1-rep(p,2))

################################################################################

# Balch's original consonant c-box, see https://sites.google.com/site/cboxbinomial/

#####################################################

# brute-force c-box for the binomial rate problem

#####################################################

# Euclidean metric in mean,sigma-space

d <- function(m,s, mm,ss) return(sqrt((m-mm)^2 + (s-ss)^2))

cbox.binomialp <- function(

x, # data in the form of counts

n = max(x), # (fixed) number of trials that resulted in the counts given in x

many = 400, # number of bootstrap replications

pp = 1 / (1 + exp(rev((-50:50)/8))) # possible p-values to explore

) {

s = length(x) # number of samples

sm = mean(x) # sample mean

ss = sd(x) # sample standard deviation

ii = 1:length(pp)

pls = rep(0,length(pp))

for (i in ii) {

p = pp[i]

# theoretical moments

tm = n * p

ts = sqrt(n * p * (1-p))

# difference between theoretical moments and observed moments

dobserved = d(tm,ts, sm,ss)

dsampling = NULL

for (k in 1:many) {

y = rbinom(s,n,p)

dsampling = c(dsampling, d(tm,ts, mean(y),sd(y)))

}

pls[[i]] = length(dsampling[dobserved <= dsampling])

}

# normalize the plausibilities

pls = pls / many

if (max(pls)<1) cat('Plausibilities do not achieve unity\n')

list(pp=pp, pls=pls / max(pls))

}

#a = cbox.binomialp(k,n,100) # doesn't seem to work

#plot(a$pp,a$pls,ylim=c(0,1),xlab='Probability', ylab='Plausibility',type='b')

showcbox = function(a,new=FALSE,col=1,lwd=1) {

if (new) plot(NULL,xlim=c(0,1),ylim=c(0,1),xlab='', ylab='Plausibility',col=col,lwd=lwd)

lines(a$pp,a$pls,type='l',col=col,lwd=lwd)

}

a1 = cbox.binomialp(c(rep(1,k), rep(0,n-k)),1,1000); showcbox(a1, new=TRUE)

a2 = cbox.binomialp(c(rep(1,k), rep(0,n-k)),1,1000); showcbox(a2)

a3 = cbox.binomialp(c(rep(1,k), rep(0,n-k)),1,1000); showcbox(a3)

a4 = cbox.binomialp(c(rep(1,k), rep(0,n-k)),1,1000); showcbox(a4)

a5 = cbox.binomialp(c(rep(1,k), rep(0,n-k)),1,1000); showcbox(a5)

q = cbox.binomialp(c(rep(1,k), rep(0,n-k)),1,10000);

# ap should be the p% confidence interval

################################################################################

# c-boxes, see https://sites.google.com/site/confidenceboxes/software

# x[i] ~ Bernoulli(parameter), x[i] is either 0 or 1

nextvalue.bernoulli <- function(x) {n <- length(x); k <- sum(x); return(env(bernoulli(k/(n+1)), bernoulli((k+1)/(n+1))))}

parameter.bernoulli <- function(x) {n <- length(x); k <- sum(x); return(env(beta(k, n-k+1), beta(k+1, n-k)))}

# these functions may be replaced with better representations of distributions that make subsequent calculations easier

many = 2000 # increase for more accuracy

beta <- function(v,w) if (v==0) rep(0,many) else if (w==0) rep(1,many) else rbeta(many, v, w)

bernoulli <- function(p) runif(many) < p

betabinomial <- function(size,v,w) rbinom(many, size, beta(v,w))

env <- function(x,y) c(x,y)

plotbox <- function(b,new=TRUE,col='blue',lwd=2,xlim=range(b[is.finite(b)]),ylim=c(0,1),xlab='',ylab='Prob',...) {

edf <- function (x, col, lwd, ...) {

n <- length(x)

s <- sort(x)

lines(c(s[[1]],s[[1]]),c(0,1/n),lwd=lwd,col=col,...)

for (i in 2:n) lines(c(s[[i-1]],s[[i]],s[[i]]),c(i-1,i-1,i)/n,col=col,lwd=lwd,...)

}

b = ifelse(b==-Inf, xlim[1] - 10, b)

b = ifelse(b==Inf, xlim[2] + 10, b)

if (new) plot(NULL, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

if (length(b) < many) edf(b,col,lwd) else

edf(c(min(b),max(b),b[1:min(length(b),many)]),col,lwd)

if (many < length(b)) edf(c(min(b),max(b),b[(many+1):length(b)]),col,lwd)

}

ci <- function(b, c=0.95, alpha=(1-c)/2, beta=1-(1-c)/2) {

left = sort(b[1:many])[round(alpha*many)]

if (many < length(b)) right = sort(b[(many+1):length(b)])[round(beta*many)] else

right = sort(b[1:many])[round(beta*many)]

c(left,right)

}

# p-boxy c-box from

c = parameter.bernoulli(c(rep(1,k), rep(0,n-k)))

plotbox(c)

cp = ci(c,c=p)

lines(cp,rep(p,2))

# c = c(c, c(c,1))

################################################################################

#Walley's IBM

#s = 1

#t = 0;

#w1 = beta(s*t+k, s*(1-t)+n-k)

#t = 1;

#w2 = beta(s*t+k, s*(1-t)+n-k)

#w = env(w1,w2)

#plotbox(w)

s = 1

W = env(beta(k, s+n-k), beta(s+k, n-k))

plotbox(W)

plotbox(c, col='lightgray', new= FALSE)

################################################################################

par(mfrow=c(3,2))

xx=(1:100)/100

# A brute-force estimate

showcbox(a1, new=TRUE,col='gray90')

showcbox(a2,col='gray95')

showcbox(a3,col='gray95')

showcbox(a4,col='gray95')

showcbox(a5,col='gray95')

lines(xx, b, type='l',xlab='',ylab='Plausibility',col='gray90')

showcbox(q,lwd=4,col='red')

# Balch

showcbox(q, new=TRUE,col='gray90')

lines(xx, b, type='l',xlab='',ylab='Plausibility',lwd=4,col='red')

# C-box

plotbox(c,ylab='Probability',xlim=c(0,1))

# Walley

plotbox(W,ylab='Probability')

bj = pbeta(xx,0.5+k, 0.5+n-k)

bu = pbeta(xx,1+k, 1+n-k)

# Jeffreys prior beta

plot(xx, bu, type='l',xlab='',ylab='Plausibility',col='gray80')

#plotbox(c,ylab='Probability',xlim=c(0,1),new=FALSE,col='gray80',lwd=1)

#lines(xx, bj, type='l',xlab='',ylab='Plausibility',col='gray80')

lines(xx, bj, type='l',xlab='',ylab='Plausibility',lwd=4)

# Uniform prior beta

plot(xx, bj, type='l',xlab='',ylab='Plausibility',col='gray80')

#plotbox(c,ylab='Probability',xlim=c(0,1),new=FALSE,col='gray80',lwd=1)

#lines(xx, bu, type='l',xlab='',ylab='Plausibility',col='gray80')

lines(xx, bu, type='l',xlab='',ylab='Plausibility',lwd=4)