Post date: Oct 06, 2015 10:50:7 AM
<<Scott wrote to Michael>>
I have one substantive suggestion for the paper on the uniform c-boxes. It is to apply the methods to an interval data set to illustrate how a p-box shaped c-box can arise. Doing this would induce one more figure in the manuscript, and perhaps one more paragraph.
I had assumed this would be a trivial application because the c-box is a decreasing monotone function of x. It's just the envelope of the application to the left endpoints, and the application to the right endpoints. This works because the larger the values of x, the more right-ward the c-box. Right? I did some unthinking simulations in the attached R code to prove this to myself, but the simulations were not compelling and started to make me nervous.
So if I increase some x[[i]] in the data set x, am I guaranteed not to get a larger value of the CDF
nextvalue.uniform <- function(x) {
n=length(x);
w=(max(x)-min(x))/beta(n-1,2);
m=(max(x)-w/2)+(w-(max(x)-min(x)))*uniform(0,1);
return(uniform(m-w/2, m+w/2))
}
for any quantile?
I guess that length(x) and beta(n-1,2) and uniform(0,1) do not matter, so we can ignore them, and thus simplify the problem to whether
w=(max(x)-min(x))
m=(max(x)-w/2)+(w-(max(x)-min(x)))
return(uniform(m-w/2, m+w/2))
is sure to be to the left if we increase any element of x. Calculus with min and max is hard, so I tried to do this my checking with the following code
mon = function(bump = 10, x = runif(20,1, 10)) {
cat(x,'\n')
y = x
w=(max(x)-min(x))
m=(max(x)-w/2)+(w-(max(x)-min(x)))
A = m-w/2
B = m+w/2
for (i in 1:some) {
x[[i]] = x[[i]] + bump
w=(max(x)-min(x))
m=(max(x)-w/2)+(w-(max(x)-min(x)))
C = m-w/2
D = m+w/2
fail = ((C-A) < 0) | ((D-B) < 0)
if (fail) cat(C-A, D-B, A,B,C,D,i,'\n')
x = y
}
return(!fail)
}
mon()
mon(100)
mon(0.1)
This always seems to return TRUE, so I guess the function is monotone, but I still don't really 'see' it in a mathematical sense.
I am on Benedryl and cough syrup, so I am especially stupid lately.
Scott
<<attached R code>>
many = 20000 # 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)
uniform <- function(a,b) runif(many, a,b)
env <- function(x,y) c(x,y)
# x[i] ~ uniform(minimum, maximum), or
# x[i] ~ uniform(midpoint, width)
nextvalue.uniform <- function(x) {n=length(x); w=(max(x)-min(x))/beta(length(x)-1,2);
m=(max(x)-w/2)+(w-(max(x)-min(x)))*uniform(0,1); return(uniform(m-w/2, m+w/2))}
parameter.uniform.minimum <- function(x) {n=length(x); w=(max(x)-min(x))/beta(length(x)-1,2);
m=(max(x)-w/2)+(w-(max(x)-min(x)))*uniform(0,1); return(m-w/2)}
parameter.uniform.maximum <- function(x) {n=length(x); w=(max(x)-min(x))/beta(length(x)-1,2);
m=(max(x)-w/2)+(w-(max(x)-min(x)))*uniform(0,1); return(m+w/2)}
parameter.uniform.width <- function(x) return((max(x)-min(x))/beta(length(x)-1,2))
parameter.uniform.midpoint <- function(x) {w=(max(x)-min(x))/beta(length(x)-1,2);
return((max(x)-w/2)+(w-(max(x)-min(x)))*uniform(0,1))}
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)
}
m = 0
M = 1
some = 20
x = runif(some, m , M)
y = x + runif(some, 0, 1.2 - x)
b = nextvalue.uniform(x)
c = nextvalue.uniform(y)
plot(NULL, xlim=c(-1,2),ylim=c(0,1))
plotbox(b,new = FALSE, col='red')
plotbox(c,new = FALSE, col='red')
o = order(x)
xo = x[o]
yo = y[o]
so = trunc(some / 2)
z = c(xo[1:so], yo[(so+1):some])
w = c(yo[1:so], xo[(so+1):some])
zz = nextvalue.uniform(z)
ww = nextvalue.uniform(w)
plotbox(zz,new = FALSE, col='blue')
plotbox(ww,new = FALSE, col='blue')
plotbox(b,new = FALSE, col='red')
plotbox(c,new = FALSE, col='red')
###################################################
# does c-box of left bounds bound c-boxes?
plot(NULL, xlim=c(-1,2),ylim=c(0,1))
plotbox(b,new = FALSE, col='red')
o = order(x)
for (i in 1:some) {
z = x
oo = o[[i]]
z[[oo]] = y[[oo]]
zz = nextvalue.uniform(z)
plotbox(zz,new = FALSE, col='blue')
}
plotbox(b,new = FALSE, col='red')
###################################################
many = 2000 # increase for more accuracy
x = runif(some, m , M)
plot(NULL, xlim=c(-1,2),ylim=c(0,1))
x=sort(x)
a=seq(0,2,length.out=100)
c = y = NULL
for (A in a) {
x[[1]] = A
cat(x,'\n')
b = nextvalue.uniform(x)
if (!is.null(y)) plotbox(y,new = FALSE, col='gray')
if (!is.null(c)) plotbox(c,new = FALSE, col='blue')
c = b
y = x
plotbox(b,new = FALSE, col='red')
plotbox(x,new = FALSE, col='black')
}
###################################################
many = 2000 # increase for more accuracy
x = runif(some, m , M)
y = x + runif(some, 0, 1.2 - x)
plot(NULL, xlim=c(-1,2),ylim=c(0,1))
plotbox(x,new = FALSE, col='gray')
plotbox(y,new = FALSE, col='gray')
b = nextvalue.uniform(x)
c = nextvalue.uniform(y)
plotbox(b,new = FALSE, col='red')
plotbox(c,new = FALSE, col='red')
# lines(c(m,M),c(0,1),lwd=2)
# o = order(x)
# xo = x[o]
# yo = y[o]
# for (i in 1:length(xo)) lines(c(xo[[i]],yo[[i]]),rep(i/(length(xo)+1),2))
for (k in 1:1000) {
z = x + runif(some) * (y - x)
d = nextvalue.uniform(z)
plotbox(d,new = FALSE, col='blue')
}
plotbox(b,new = FALSE, col='red')
plotbox(c,new = FALSE, col='red')
###################################################