Embedded Files
THIS PAGE IS OBSOLETE AND SCHEDULED FOR REMOVAL. THE NEW SITE IS HERE.
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)
Page updated
Google Sites
Report abuse