Estimating both binomial parameters

The function below, designed for use with R, analyzes random data assumed to have been sampled from a binomial(n, p) probability distribution to create a joint confidence structure for the unknown parameters n and p, where n is a positive integer, and p is a probability between zero and one.  See also the simpler problem of estimating just the p.

The graphs created depict non-consonant confidence structures, thus expressing them as p-boxes will lose some information. The regions in n,p-space with plausibilities greater than α represent a 100(1−α)% confidence region for the binomial parameters n and p. 

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

# c-boxes for the binomial(n,p) problem 

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

# Euclidean metric in mean,sigma-space 

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

cbox.binomialnp <- function( 

     x,                      # data in the form of counts 

     maxn = 20*max(x),       # largest n we're considering 

     many = 4000             # number of bootstrap replications 

     ) { 

s = length(x)                # number of samples 

sm = mean(x)                 # sample mean

ss = sd(x)                   # sample standard deviation

nn = max(x):maxn 

some = -30:30

ii = 1:length(some) 

pp = 1 / (1 + exp(rev(some/8))) 

pls = matrix(0,nrow=max(nn),ncol=length(some)) 

for (n in nn) 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[n,i] = length(dsampling[dobserved <= dsampling]) 

   } 

# normalize the plausibilities 

pls = pls / many 

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

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

} 

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

# depict and rotate the two-dimensional distribution 

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

show <- function(phi,theta, binomialnp.answer) { 

   maxn = binomialnp.answer$maxn 

   pp = binomialnp.answer$pp 

   pls = binomialnp.answer$pls 

   persp(1:maxn,pp,pls,ticktype='detailed', phi=phi, theta=theta, ylim=c(0,1),

      xlab='Number of trials', ylab='Probability', zlab='Plausibility',

      col='lightblue', shade = 0.85, ltheta=20) -> res 

   }

# example

x = c(6, 6, 4, 8, 5, 4, 6, 7)                                # generated as rbinom(8, 10, 0.6213)

answer = cbox.binomialnp(x,40)                 # will often take several minutes to compute

for (s in -120:120) show(s,-20, answer)      # gratuitous animation

for (s in -120:120) show(20,s, answer)       # gratuitous animation

show(20,155,answer) 

# another example

x = c(3, 3, 4, 3, 4, 4, 3, 3, 11, 4, 

         5, 2, 3, 2, 3, 6, 4, 4, 6, 4, 1, 

         3, 4, 3, 7, 5, 3, 5, 2, 6, 3, 5, 

         1, 5, 3, 1, 4, 4, 5, 4, 5, 3, 4, 

         1, 5, 6, 4, 5, 3, 5, 1, 5)                   # generated as x = rbinom(52, 20, 0.2) 

a = cbox.binomialnp(x,40)                  # will often take several minutes to compute

for (s in -120:120) show(s,-20, a)       # gratuitous animation

for (s in -120:120) show(20,s, a)        # gratuitous animation

show(20,-20,a)