In many risk analyses, distribution tails are the focus of much attention because the extreme events represented by these tails correspond to the adverse events of most concern. The algorithms currently used in Risk Calc and our other software implementations for performing convolutions are not well suited to discerning the structure of the extreme tails of distributions and p-boxes. The algorithms, which are mostly taken from Yager (1986) and Williamson and Downs (1990), work with the quasi-inverse functions of cumulative distribution functions and use an equiprobability discretization. For instance, if there are 100 discretization steps, then we can access 101 quantile values, corresponding to the probability levels 0, 0.01, 0.02, 0.03, …, 0.98, 0.99, and 1.
In principle, we could increase the number of discretization steps to access finer details about the quantiles in the tails of distributions and p-boxes. But doubling the number of discretization steps to 200 only gives us access to quantiles at the probability levels 0, 0.005, 0.01, 0.015, …, 0.995, 1. Even if we use 2000 discretization steps (which is close to the practical limit in R for computing convolutions assuming independence), we can only get quantiles at probabilities 0, 0.0005, 0.001, 0.0015, 0.002, …, 0.9995, 1. Highly unlikely outcomes with probabilities like 10-6 or 10-12 are lost within the first discretization interval. Thus, our computational methods may be too crude for discerning tail risks when the disaggregation models or fault tree analyses are focused on very unlikely events. Mathematical rescaling or reconditionalization tricks could allow us to make the needed calculations, but they can be confusing and awkward.
We want to derive and implement tail-emphasizing algorithms to replace the current approach based on equiprobability discretization for computing mathematical operations (sums, differences, products, quotients, conjunctions, disjunctions, etc.) on uncertain numbers. For instance, it might be possible to use a mixed scheme maintaining equal-probability steps in middle ranges but order-of-magnitude steps in the tails. Such a scheme that kept track of quantile bounds for the probability levels 0, 10-20, 10-19, …, 0.0001, 0.001, 0.01, 0.02, 0.03, 0.04, …, 0.98, 0.99, 0.999, 0.9999, …, 1-10-19, 1-10-20, and 1, would be able to make quantitative assertions about risks between zero and 10-20 even though it uses only 1+18+99+18+1 = 137 discretization levels. An alternative discretization scheme might be to augment the percent levels (0, 0.01, 0.02, …, 0.98, 0.99, 1) with additional levels for the left tail at 0.005, 0.0025, 0.00125, 0.000625, …, and similar values in the right tail at 0.995, 0.9975, 0.99875, 0.999375, …. It would take 60 halvings to reach the probability level of 10-20.
This effort will likely involve three subtasks, including
designing the data structure to hold these quantiles [this may be a minor task],
rewriting the convolution algorithms to work with the new data structure, and
revising the plotting algorithm to display a p-box in the new data structure.
It may also be useful to implement a scheme to translate p-boxes in the old data structure to and from the new data structure via outward-directed rounding. The scheme might also be useful for the condensation step in convolutions assuming independence.
It might be smart, for instance, to
ii <- function() c( 0, 1:(pbox.steps-1) / pbox.steps)
jj <- function() c( 1:(pbox.steps-1) / pbox.steps, 1)
iii <- function() c( pbox.bOt, 1:(pbox.steps-1) / pbox.steps)
jjj <- function() c( 1:(pbox.steps-1) / pbox.steps, pbox.tOp)
qq = NULL; a = 0.01; for (i in 1:20) {a = a/2; qq = c(qq,a)}
q = c(0,rev(qq),1:99/100,1-qq,1)
q = c(0,rev(qq),1:99/100,1-qq,1)
y = qunif(q,0,1000)
z = qunif(q,-1000,0)
y = qnorm(q,0,1000)
z = qnorm(q,-1000,0)
qq = 10^((-10):(-3))
q = c(0,qq,1:99/100,1-rev(qq),1)
q = c( qq,1:99/100,1-rev(qq) )
x = qunif(q,0,1000)
x = qnorm(q,0,1000)
iii <- ii <- function() q[-length(q)] # c( 0, 1:(pbox.steps-1) / pbox.steps)
jjj <- jj <- function() q[-1] # c( 1:(pbox.steps-1) / pbox.steps, 1)
# iii <- function() c( pbox.bOt, 1:(pbox.steps-1) / pbox.steps)
# jjj <- function() c( 1:(pbox.steps-1) / pbox.steps, pbox.tOp)
u = qunif(ii(),0,1000)
d = qunif(jj(),0,1000)
u = qnorm(ii(),5,1)
d = qnorm(jj(),5,1)
k = is.finite(u) & is.finite(d)
steps = length(u)
old.par <- par(mfrow=c(2,2),mar=c(2,4,3,1))
# wrong plot, but it “shows” those tails
# plot(NULL,ylab='Cumulative probability', ylim=c(0,1),xlim=range(c(u,d)))
# lines(c(u,u[[steps]],d[[steps]]),c(0:(steps-1) / steps,1,1),type="S")
# lines(c(u[1],d[1],d),c(0,0,1:(steps) / steps),type="s")
plotit <- function(u,d) {
lines(c(u,u[[steps]],d[[steps]]),c(ii(),1,1),type="S")
lines(c(u[1],d[1],d),c(0,0,jj()),type="s")
}
plot(NULL,ylab='Cumulative probability', ylim=c(0,1),xlim=range(c(u,d)))
plotit(u,d)
#lines(c(0,5,5,0,0),c(0,0,0.1,0.1,0),col='red')
plot(NULL,ylab='Cumulative probability', ylim=c(0,0.1),xlim=range(c(0,5)))
plotit(u,d)
#lines(c(0,4,4,0,0),c(0,0,0.01,0.01,0),col='red')
plot(NULL,ylab='Cumulative probability', ylim=c(0,0.01),xlim=range(c(0,4)))
plotit(u,d)
#lines(c(0,3,3,0,0),c(0,0,0.001,0.001,0),col='red')
plot(NULL,ylab='Cumulative probability', ylim=c(0,0.001),xlim=range(c(0,3)))
#plotit(u,d)