The following is an R script that exercises and checks the Quick Bayes c-box for the normal distribution. It demonstrates that the c-boxes for the mean and standard deviation perform statistically in the sense that they produce a uniform map in a Singh plot. It also demonstrates that one can compute with confidence boxes by forming sums, differences, and―with certain limitations―products, quotients, minima, maxima, and powers. In particular, note that subtracting c-boxes for normal means solves the Behrens−Fisher problem. The results also have the confidence interpretation as shown by Singh plots. This script was created as a part of our quality assurance efforts for Quick Bayes software. You can run the script and see its output by sequentially pasting its sections (which are delimited by rows of pound signs #) into the R console.
<<THIS CODE IS CURRENTLY IN REVISION>>
############################################################################################ red
# NORMAL MEAN
# Gosset was correct...demonstrate how the scaled t-distribution is a confidence distribution,
# and thus a c-box, for the mean of a normal distribution
# c-box and confidence distribution for the normal mean
cd = function(z) {n = length(z); mean(z) + sd(z) * rt(many,n-1) / sqrt(n)}
edf = function(x) lines(sort(x), 1:length(x)/length(x),type = 's')
some = 5000
many = 5000
mu = runif(1,-100,100)
sg = runif(1,0,100)
cat('normal(',mu,',',sg,')\n',sep='')
nn = c(2,3,4,5,6,10,100,1000,10000)
par(mfrow=c(3,3))
for (n in nn) {
u = array(0, some)
for( i in 1:some ){
x = rnorm( n , mu , sg )
u[i] = sum( cd(x) < mu )/many
}
plot( c(0,1) , c(0,1) , type='l' , col=2)
edf(u)
title(paste('n',n))
}
############################################################################################ red
# GOSSET FOR GIANT VARIANCES
n = 10
mu = 2
sgsg = c(0.0001,0.01,0.5,1,2,5,10,20,100)
many = 5000
some = 5000
cd = function(z) {n = length(z); mean(z) + sd(z) * rt(many,n-1) / sqrt(n)}
edf = function(x) lines(sort(x), 1:length(x)/length(x),type = 's')
par(mfrow=c(3,3))
for (sg in sgsg) {
u = array(0, some)
for( i in 1:some ){
x = rnorm( n , mu , sg )
u[i] = sum( cd(x) < mu)/many
}
plot( c(0,1) , c(0,1) , type='l' , col=2)
edf(u)
title(paste('sigma',sg))
}
############################################################################################ red
# ADDITION OF NORMAL MEANS
f = function(a,b) a+b
cd = function(x) {n = length(x); mean(x) + sd(x) * rt(many,n-1) / sqrt(n)}
edf = function(x) lines(sort(x), 1:length(x)/length(x),type = 's')
n1 = 8
n2 = 8
sigma1 = 1
sigma2 = 1
some = 2000
many = 2000
mus = c(0.0, 0.1, 0.2, 0.3, 0.4)
m = length(mus)
par(mfrow=c(5,5))
for (j in 1:m) for (k in 1:m) {
mu1 = mus[[j]]
mu2 = mus[[k]]
u = array(0, some)
for (i in 1:some) {
x = rnorm(n1, mu1, sigma1) # sample data
y = rnorm(n2, mu2, sigma2) # sample data
d = f(cd(x), cd(y))
u[i] = sum(d < f(mu1,mu2)) / many
}
plot(c(0,1),c(0,1),type='l', col='red')
edf(u)
title(paste('means',mu1,mu2))
cat(j, k, n1,mu1,sigma1,'plus', n2,mu2,sigma2,'\n')
}
############################################################################################ red
# SUBTRACTION OF NORMAL MEANS
# this is the Behrens-Fisher problem
# The Wikipedia article on the subject says:
#
# In statistics, the Behrens–Fisher problem, named after Walter Ulrich Behrens and
# Ronald Fisher, is the problem of interval estimation and hypothesis testing concerning
# the difference between the means of two normally distributed populations [and whether
# the means can reasonably be treated as equal] when the variances of the two populations
# are not assumed to be equal, based on two independent samples.
#
# Solutions to the Behrens–Fisher problem have been presented that make use of either
# a classical or a Bayesian inference point of view and either solution would be notionally
# invalid judged from the other point of view. If consideration is restricted to classical
# statistical inference only, it is possible to seek solutions to the inference problem that
# are simple to apply in a practical sense, giving preference to this simplicity over any
# inaccuracy in the corresponding probability statements. Where exactness of the
# significance levels of statistical tests is required, there may be an additional
# requirement that the procedure should make maximum use of the statistical
# information in the dataset. It is well known that an exact test can be gained by
# randomly discarding data from the larger dataset until the sample sizes are equal,
# assembling data in pairs and taking differences, and then using an ordinary t-test
# to test for the mean-difference being zero: clearly this would not be "optimal" in any
# sense.
#
# The task of specifying interval estimates for this problem is one where a frequentist
# approach fails to provide an exact solution, although some approximations are available.
# Standard Bayesian approaches also fail to provide an answer that can be expressed as
# straightforward simple formulae, but modern computational methods of Bayesian
# analysis do allow essentially exact solutions to be found. Thus study of the problem
# can be used to elucidate the differences between the frequentist and Bayesian approaches
# to interval estimation.
f = function(a,b) a-b
cd = function(x) {n = length(x); mean(x) + sd(x) * rt(many,n-1) / sqrt(n)}
edf = function(x) lines(sort(x), 1:length(x)/length(x),type = 's')
some = 2000
many = 2000
n1 = round(runif(1,1,10))
n2 = round(runif(1,1,10))
cat('Sample sizes',n1,n2,'\n')
sigma1 = runif(1,0,10)
sigma2 = runif(1,0,20) # Behrens-Fisher does not assume equal variances
cat('Sigmas',sigma1,sigma2,'\n')
mus = c(0, 1, 5, 10, 50)
m = length(mus)
par(mfrow=c(5,5))
for (j in 1:m) for (k in 1:m) {
mu1 = mus[[j]]
mu2 = mus[[k]]
u = array(0, some)
for (i in 1:some) {
x = rnorm(n1, mu1, sigma1) # sample data
y = rnorm(n2, mu2, sigma2) # sample data
d = f(cd(x), cd(y))
u[i] = sum(d < f(mu1,mu2)) / many
}
plot(c(0,1),c(0,1),type='l', col='red')
edf(u)
title(paste('means',mu1,mu2))
cat(j, k, n1,mu1,sigma1,'plus', n2,mu2,sigma2,'\n')
}
############################################################################################ red
# Check that the IMPLEMENTATION of SUBTRACTION in R matches that in RAMAS Risk Calc
##############################
# RUN THIS PART IN R
##############################
f = function(a,b) a-b
cd = function(x) {n = length(x); mean(x) + sd(x) * rt(many,n-1) / sqrt(n)}
edf = function(x) lines(sort(x), 1:length(x)/length(x),type = 's')
some = 2000
many = 20000
n1 = 5
n2 = 10
mu1 = 4
mu2 = 9
sg1 = 3
sg2 = 12
x = rnorm(n1, mu1, sg1) # sample data
y = rnorm(n2, mu2, sg2) # sample data
x = c(3.30171812857776, 6.43717247556227, 5.16202091771867, 0.0882133027828953,
5.712631788647)
y = c(-19.7414623138356, -24.4788374324584, 8.44049162542419, -13.2825346398374,
13.1135893503267, 5.39433897900108, -2.20729228848036, 22.5911923579581,
10.2702860767471, 22.7468609896814)
d = f(cd(x), cd(y))
plot(NULL, xlim=range(d),ylim=c(0,1))
edf(d)
dput(x) # c(3.30171812857776, 6.43717247556227, 5.16202091771867, 0.0882133027828953, 5.712631788647)
dput(y) # c(-19.7414623138356, -24.4788374324584, 8.44049162542419, -13.2825346398374,13.1135893503267, 5.39433897900108, -2.20729228848036, 22.5911923579581,10.2702860767471, 22.7468609896814)
mean(x) # 4.140351
mean(y) # 2.284663
sd(x) # 2.545309
sd(y) # 16.74517
mean(d) # 1.886195
sd(d) # 6.194256
##############################
# RUN THIS PART IN RAMAS RISK CALC
##############################
#tOp = 0.99995
#bOt = 0.00005
#n1 = 5
#n2 = 10
#meanx = 4.140351
#meany = 2.284663
#sdx = 2.545309
#sdy = 16.74517
#x = meanx + sdx * t(n1-1) / sqrt(n1)
#y = meany + sdy * t(n2-1) / sqrt(n2)
#z = x |-| y
#mean(z) // 1.855688
#stddev(z) // 6.216346201
############################################################################################ red
# MULTIPLICATION OF NORMAL MEANS
f = function(a,b) a*b
cd = function(x) {n = length(x); mean(x) + sd(x) * rt(many,n-1) / sqrt(n)}
edf = function(x) lines(sort(x), 1:length(x)/length(x),type = 's')
areametric = function(u) sum(abs(sort(u)-((1:length(u))/length(u))))/length(u)
n1 = 8
n2 = 8
sigma1 = 1
sigma2 = 1
some = 2000
many = 2000
mus = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.2, 1.6, 2, 3, 4, 6)
mus = c(0.0, 0.3, 0.6, 0.9, 1.2)
m = length(mus)
A = matrix(Inf, nrow=length(mus),ncol=length(mus))
par(mfrow=c(5,5))
for (j in 1:m) for (k in 1:m) {
mu1 = mus[[j]]
mu2 = mus[[k]]
u = array(0, some)
for (i in 1:some) {
x = rnorm(n1, mu1, sigma1) # sample data
y = rnorm(n2, mu2, sigma2) # sample data
d = f(cd(x), cd(y))
u[i] = sum(d < f(mu1,mu2)) / many
}
a = areametric(u)
A[j,k] = a
plot(c(0,1),c(0,1),type='l', col='red')
edf(u)
title(paste('means',mu1,mu2))
cat(j, ' ',k, ' ',a, ' ',n1,'~N(',mu1,sigma1,') times ', n2,'~N(',mu2,sigma2,')\n',sep='')
}
#par(mfrow=c(1,1))
#persp(mus,mus,A,theta=85)
#A
############################################################################################ red
# DIVISION OF NORMAL MEANS
f = function(a,b) a/b
cd = function(x) {n = length(x); mean(x) + sd(x) * rt(many,n-1) / sqrt(n)}
edf = function(x) lines(sort(x), 1:length(x)/length(x),type = 's')
areametric = function(u) sum(abs(sort(u)-((1:length(u))/length(u))))/length(u)
some = 2000
many = 2000
m = 5
A = matrix(Inf, nrow=m, ncol=m)
par(mfrow=c(m,m))
for (j in 1:m) for (k in 1:m) {
n1 = round(runif(1,2,10))
n2 = round(runif(1,2,10))
mu1 = signif(runif(1,-15,15),3)
mu2 = signif(runif(1,-15,15),3)
sigma1 = signif(runif(1,0,10),3)
sigma2 = signif(runif(1,0,10),3)
u = array(0, some)
for (i in 1:some) {
x = rnorm(n1, mu1, sigma1) # sample data
y = rnorm(n2, mu2, sigma2) # sample data
d = f(cd(x), cd(y))
u[i] = sum(d < f(mu1,mu2)) / many
}
a = areametric(u)
A[j,k] = a
plot(c(0,1),c(0,1),type='l', col='red')
edf(u)
if (abs(mu2)<sigma2) text(0.2,0.5,'PROBLEM')
title(paste(mu1,sigma1,'/',mu2,sigma2))
cat(j, ' ',k, ' ',a, ' ',n1,'~N(',mu1,',',sigma1,') times ', n2,'~N(',mu2,',',sigma2,')\n',sep='')
}
# It seems the only restriction is not to divide by zero. That is, apparently, divisions of
# normal means are fine so long as the true mean of the divisor (denumerator) has an
# absolute value larger than its true standard deviation. The distribution of the dividend
# (numerator) can be anything; it can be positive or negative or straddle zero. The divisor
# distribution can be positive or negative, but it cannot straddle zero.
#
# Of course, analysts do not know the true parameter values, so we should devise a test
# based on data to discern whether the result of a division is likely to be reliable. Although
# the basic effect of such a test is obvious, it is not clear how it should be formally designed.
#
# When the divisor is near zero and the quotient is negative (because just one of the terms
# is negative), the quotient is likely over-estimated and the Singh plot will be above the unit
# uniform. When this bad quotient is positive (because the terms share the same sign), the
# quotient is likely under-estimated and the Singh plot will be below the unit uniform. One
# might think that these facts could be wrangled to compute crude but reliable bounds on
# the quotient. Unfortunately, the nature of the misestimation does not allow this. Whatever
# the sign of the quotient, the estimate is closer to zero than it should be. There does not
# seem to be a way to glean how much absolutely bigger the estimate should be. The
# estimates already have long, fat tails because we are essentially dividing by zero to
# estimate them.
#
# The following code snippet illustrates these conclusions. It makes three graphs which
# display (i) the c-box underestimating the true quotient m=5, (ii) the would-be imprecise
# c-box obtained by enveloping the first c-box with its negation, and (iii) the natural scale
# that R choses to display these objects.
m = 5
par(mfrow=c(3,1))
x = rnorm(8,m,1)
y = rnorm(7,1,5)
q = f(cd(x), cd(y))
p = -f(cd(x), cd(y))
plot(NULL,xlim=range(c(-10,10)),ylim=c(0,1))
edf(q)
abline(v=m)
plot(NULL,xlim=range(c(-10,10)),ylim=c(0,1))
edf(q)
edf(p)
abline(v=m)
plot(NULL,xlim=range(c(p,q)),ylim=c(0,1))
edf(q)
edf(p)
abline(v=m)