Here is an approach for the computation of the intervals using a fiducial distribution. The steps are as following
Compute fiducial distributions for the two individual percentiles. This can be done by considering the probabilities that the true percentile $\theta_{90}$ is below a particular order statistic $P(\theta_{90}<x_{(i)})$. This probability can be computed with a beta distribution.
See also: How to obtain a confidence interval for a percentile?
Interpret the fiducial distributions as probabilitiy distributions (this is an approximation, the fiducial distribution does not behave exactly like a probability distribution) and compute the probabilities for the joint probabilities of the two parameters/percentiles to be inside the 2-d bins/cells created by the grid.
Compute boundaries based on the distribution and consider the bounds for the difference as boundary of regions that ensure the probability to be at least $\alpha \%$ (the computed joint probability is coarse/discrete and we can not nicely integrate it, but we can sum probabilities of all the cells through which the boundary passes)
The image below illustrates how this goes. For two samples $x$ and $y$ we draw a grid with grid lines corresponding to the values of the sample. We compute for each cell probabilities (depicted by a color in the image, purple is more probable and yellow less probable). For a line corresponding to some constant value of $x-y$ we consider all the cells that it passes through and in addition the cells above/below it (for lower/upper bound) and add up the probability to compute an estimate for the fiducial probability. The lines in the image below correspond to a 95% fiducial interval and it is computed such that the light shaded area is at least 95%. The actual coverage will be larger because the lines cover a larger area than the light shaded area.
Note 1: I created this approach quickly and do not know whether this has been described elsewhere before. What annoys me is that the actual coverage 0.9701 is different from the nominal coverage 0.95, even with large samples such as $n=500$. The method is an approximation, but I had hoped that the difference would have become smaller quickly with increasing sample size. I am not sure whether the discrepancy is because of the approximation errors (using discrete boundaries of the square grid) or because of the intrinsic error in the method which is the use of the fiducial distribution as if it is a true probability density independent from the sample.
Note 2: the problem with the discrete estimate of the density might be reduced by smoothening the distribution. I have left this out as it complicates the problem a lot and is not really necessary for a proof of principle.
Note 3: The computation with the matrices for the joint distribution takes a lot of time. As an alternative one can also use a Monte Carlo simulation after having obtained the one-dimensional fiducial distributions. Also, for repetitions with the same sample sizes and quantiles, the computations with the beta distributions can be stored and re-used.
Edit: I have implemented notes 2 and 3 and it appears that now for n=500 the actual coverage might be closer to the nominal coverage. This indicates that the discreteness of the fiducial distribution (not allowing to select exact 95% regions) plays a big role.
R-Code:
### function to compute fi% fiducial interval for difference in q-th quantiles
fi_dq = function(x, y, q = 0.9, fi = 0.95) {
nx = length(x)
ny = length(y)
x = x[order(x)]
y = y[order(y)]
kx = 1:nx
ky = 1:ny
compute one-dimensional fiducial probabilities and related probabilities
Fx = 1-pbeta(q,kx,nx+1-kx) ## P(theta < x_k)
Fy = 1-pbeta(q,ky,ny+1-ky) ## P(theta < y_k)
px = diff(c(0,Fx,1))
py = diff(c(0,Fy,1))
pxy = outer(px,py,"*") ### cell probabilities
diffxy = outer(x,y,"-") ### cell values of difference
compute fiducial lower boundaries
lowerxy = rbind(diffxy[1,],diffxy) ### add rows for edges
lowerxy = cbind(lowerxy,lowerxy[,ny])
lowerord = order(lowerxy) ### order values and add upp associated probabilities
lower_diff = lowerxy[lowerord]
lower_F = cumsum(pxy[lowerord])
lower_k = max( which(lower_F < (1-fi)/2) ) ### compute a 1-alpha/2 boundary
lower_bound = lower_diff[lower_k]
compute fiducial upper boundaries
upperxy = rbind(diffxy,diffxy[nx,]) ### add rows for edges
upperxy = cbind(upperxy[,1],upperxy)
upperord = rev(order(upperxy))
upper_diff = upperxy[upperord]
upper_F = cumsum(pxy[upperord])
upper_k = max( which(upper_F < (1-fi)/2) )
upper_bound = upper_diff[upper_k]
return(c(lower_bound, upper_bound))
}
some settings
set.seed(1)
n = 500
q = 0.9
fi = 0.95
simulations to test true coverage probability
m = 10^4
set.seed(1)
bounds = sapply(1:m, FUN = function(s) {
x = rexp(n,1/10)
y = rexp(n,1/15)
fi_dq(x, y , q, fi)
})
effective coverage
true = qexp(q,1/10) - qexp(q,1/15)
mean((bounds[1,] < true)*(true < bounds[2,]))
example from question
results in -0.8429189 2.7977309
set.seed(123)
data <- data.frame(
Group = sample(c("Before", "After"), 100, replace = TRUE),
ResponseTime = rnorm(100, mean = c(10, 15), sd = 2)
)
fi_dq(x = data$ResponseTime[data$Group == "After"],
y = data$ResponseTime[data$Group == "Before"],
q = 0.9, fi = 0.95)
illustration of joint fiducial distribution
set.seed(1)
n = 100
x = rexp(n,1/10)
y = rexp(n,1/15)
q = 0.9
fi = 0.95
nx = length(x)
ny = length(y)
x = x[order(x)]
y = y[order(y)]
kx = 1:nx
ky = 1:ny
Fx = 1-pbeta(q,kx,nx+1-kx) ## P(theta < x_k)
Fy = 1-pbeta(q,ky,ny+1-ky) ## P(theta < x_k)
fx = diff(Fx)
fy = diff(Fy)
pxy = outer(fx,fy,"*")
diffxy = outer(x,y,"-")
plot(-100,-100, xlim = range(x), ylim = range(y),
xlab = "x", ylab = "y")
bounds = fi_dq(x, y, q, fi)
color polygons based on probabilities (not densities)
p_area = 0
for (l in c(2:nx)) {
for (m in c(2:ny)) {
area1 = sum(diffxy[l+c(-1,0),m+c(-1,0)]>bounds[1])
area2 = sum(diffxy[l+c(-1,0),m+c(-1,0)]<bounds[2])
p = pxy[l-1,m-1]
rp = (p/max(pxy))
polygon(x[l+c(-1,0,0,-1)],y[m+c(-1,-1,0,0)],
col = hsv(0.15+0.55rp,
0.7-0.5rp,
1-0.3rp),
lwd = 0.1)
if (area1 < 4) {
polygon(x[l+c(-1,0,0,-1)],y[m+c(-1,-1,0,0)],
col = hsv(0,0,0,0.2),
lwd = 0.1)
}
if (area2 < 4) {
polygon(x[l+c(-1,0,0,-1)],y[m+c(-1,-1,0,0)],
col = hsv(0,0,0,0.2),
lwd = 0.1)
}
if ((area2 == 4)(area1 == 4)) {
p_area = p_area + p
}
}
}
p_area ### manual computed check p_area should be at least 'fi' or larger
lines(c(-10,70),c(-10,70)-bounds[1], lwd = 2, col = 1)
lines(c(-10,70),c(-10,70)-bounds[2], lwd = 2, col = 1)
Alternative code that combines note 2 and note 3 by computing with monte carlo simulations and making the discrete fiducial distribution a piecewise linear continuous function.
### function to compute fi% fiducial interval for difference in q-th quantiles
### using a monte carlo approach
fi_dq_mc = function(x, y, q = 0.9, fi = 0.95, m = 10^3) {
nx = length(x)
ny = length(y)
x = x[order(x)]
y = y[order(y)]
kx = 1:nx
ky = 1:ny
compute one-dimensional fiducial probabilities and related probabilities
Fx = 1-pbeta(q,kx,nx+1-kx) ## P(theta < x_k)
Fy = 1-pbeta(q,ky,ny+1-ky) ## P(theta < y_k)
px = diff(c(Fx))
py = diff(c(Fy))
draw_xy = sapply(1:m, FUN = function(k) {
draw_x = sample(1:(nx-1),1, prob = px)
draw_y = sample(1:(ny-1),1, prob = py)
samp_x = runif(1, x[draw_x],x[draw_x+1])
samp_y = runif(1, y[draw_y],y[draw_y+1])
return(samp_x-samp_y)
})
lower_bound = quantile(draw_xy,(1-fi)/2)
upper_bound = quantile(draw_xy,1-(1-fi)/2)
return(c(lower_bound, upper_bound))
}
some settings
set.seed(1)
n = 500
q = 0.9
fi = 0.95
simulations to test true coverage probability
m = 10^3
set.seed(1)
bounds = sapply(1:m, FUN = function(s) {
x = rexp(n,1/10)
y = rexp(n,1/15)
fi_dq_mc(x, y , q, fi)
})
effective coverage
result is 0.955
true = qexp(q,1/10) - qexp(q,1/15)
mean((bounds[1,] < true)*(true < bounds[2,]))
