p-box with hyriskR
pbox.RmdTo showcase the use of p-boxes with hyriskR, we use the Rosenbrock function studied by Schoebi & Sudret (2016).
library(hyriskR)Description of an example
The assessment function is coded as follows:
FUN<-function(X){
A = X[1]
B = X[2]
return(100*(B-A^2)^2+(1-A)^2)
}Step 1: uncertainty representation
The first step focuses on uncertainty representation. It aims at selecting the most appropriate mathematical tool to represent uncertainty on the considered parameter. In the case considered, the two uncertain inputs are represented by a probability distribution with imprecise parameters, i.e. a family of parametric probability distributions represented by a p-box Ferson et al. (2003).
The procedure in hyriskR first uses the create_input function to define the input variables (imprecise, random or fixed). In this example, the tow first inputs define the p-box with the faimily of normal distribution, and the subsequent inputs define the representation of the imprecise normal parameters, i.e., the mean and the standard deviation. Here they are represented with possibility distributions. Second, the create_distr function assigns the corresponding distribution (probability or possibility) to each uncertain input.
ninput<-6
input<-vector(mode="list", length=ninput)
input[[1]] = create_input(
name = "A",
type = "impr proba",
distr = "normal",
param = c(3, 5),
monoton = "dunno"
)
input[[2]] = create_input(
name = "B",
type = "impr proba",
distr = "normal",
param = c(4, 6),
monoton = "dunno"
)
input[[3]] = create_input(
name = "mu_A",
type = "possi",
distr = "interval",
param = c(-0.5, 0.5)
)
input[[4]] = create_input(
name = "mu_B",
type = "possi",
distr = "interval",
param = c(-0.5, 0.5)
)
input[[5]] = create_input(
name = "s_A",
type = "possi",
distr = "interval",
param = c(0.7, 1)
)
input[[6]] = create_input(
name = "s_B",
type = "possi",
distr = "interval",
param = c(0.7, 1)
)
input = create_distr(input)A visualisation function plot_input has also been implemented to handle the different representation forms.
#VISUALISATION of INPUT UNCERTAINTY REPRESENTATION
plot_input(input)
Step 2: uncertainty propagation
The second step aims at conducting uncertainty propagation, i.e. evaluating the impact of the uncertainty pervading the input on the outcome of the risk assessment model. To do so, the main function is propag, which implements the Monte-Carlo-based algorithm of [@Baudrit07], named IRS (Independent Random Sampling), for jointly handling possibility and probability distributions and the algorithm of [@Baudrit08] for jointly handling possibility, probability distributions and p-boxes.
choice_opt = "L-BFGS-B_MULTI" ## quasiNewton with multistart
param_opt = 25 ## multistarts
Z0_IRS = propag(N = 500, input, FUN, choice_opt, param_opt, mode = "IRS")Step 3: post-processing of the Results
The output of the propagation procedure then takes the form of N random intervals of the form \([\underline{Z}_k,\overline{Z}_k]\), with \(k=1,...,N\). This information can be summarized in the form of a pair of upper and lower cumulative probability distributions (CDFs), in the form of a p-box which is closely related to upper and lower probabilities of Dempster, and belief functions of Shafer as proposed by Baudrit et al. 2007. The following code provides the output of the propagation phase using the plot_cdf function. In some situations, the analysts may be more comfortable in deriving a unique probability distribution. In order to support decision-making with a less extreme indicator than using either probability bounds, Dubois and Guyonnet 2011. proposed to weight the bounds by an index \(w\), which reflects the attitude of the decision-maker to risk (i.e. the degree of risk aversion) so that the resulting distribution \(F\) holds as \(F^{-1}(x) = w.\overline{F}^{-1}(x)+(1-w).\underline{F}^{-1}(x)\) where \(x\) is the quantile level ranging between 0 and 1. This can be done using the summary_1cdf function as follows using two aversion weight values of respectively 30 and 50\(\%\).
#Visualisation
plot_cdf(Z0_IRS, xlab = "Z", ylab = "Cumulative Probability", main = "", lwd = 3.5)
grid(lwd = 1.5)
#Summary with one CDF
Z50 = summary_1cdf(Z0_IRS,aversion=0.5)##risk aversion of 50%
lines(ecdf(Z50), col=4, lwd = 3.5)
Z30 = summary_1cdf(Z0_IRS,aversion=0.3)##risk aversion of 30%
lines(ecdf(Z30), col = 3, lwd = 3.5)
legend("bottomright",c("upper CDF","lower CDF","Aversion weight 50%",
"Aversion weight 30%"),col=c(1, 2, 4, 3),lwd = 2.5, bty = "n")
The functionalities of hyriskR enable to summarise the p-box depending on the statistical quantity of interest supporting the decision making process.
- If the interest is a quantile at a given level (say 75\(\%\)), the function provides the corresponding interval.
#Interval of quantiles
level = 0.75##quantile level
quant = quan_interval(Z0_IRS, level)
print(paste("Quantile inf: ", round(quant$Qlow,2)))
#> [1] "Quantile inf: 3.58"
print(paste("Quantile sup: ", round(quant$Qupp,2)))
#> [1] "Quantile sup: 1201.59"- If the interest is a global measure of epistemic uncertainty affecting the whole probability distribution of \(Z\), the uncertainty function computes the area within the lower and upper probability distribution.
#Global indicator of uncertainty
unc = uncertainty(Z0_IRS)
print(paste("Epistemic uncertainty : ", round(unc,2), sep=""))
#> [1] "Epistemic uncertainty : 104728.95"- If the interest is the probability of Z being below the decision threshold at 2,000, the proba_interval function provides the corresponding upper and lower bound on this probability.