3.8 Computational Implementation

  

Figure 3.2: Schematic depiction of the steps involved in the application of conventional Monte Carlo method

  
3.8.1 Implementation of the Conventional Uncertainty Propagation Methods

Modules for the application of conventional Monte Carlo and the Latin Hypercube Sampling methods were developed as part of this work. Though these modules are mainly used in this work for the evaluation of alternative methods, they can be used in isolation.

As shown in Figure 3.2, the conventional Monte Carlo method involves the following steps: (a) obtaining random samples from the probability distributions of the inputs, (b) performing model simulations for the combination of the sampled inputs, and (c) statistically analyzing the model outputs.

The random numbers for sampling the input distributions were generated depending on the case study considered. For all the cases, the normal random numbers were generated through the function GASDEV provided with the mathematical routines of the Numerical Recipes [163]. Gamma random variables were generated by first generating a normal random number, and applying the transformation shown in Table 3.1. In the case of Dirichlet distributions, the independent random numbers were generated first, and the corresponding ratios were calculated using Equation 3.2.4. In one specific case, in the study involving a physiologically based pharmacokinetic model, described in Section 5.1, the model was implemented on the SimuSolv modeling platform [58], and the Monte Carlo method involved sampling from random numbers generated by the pseudo-gaussian random number generator function in SimuSolv, R-GAUSS.

In all the case studies described in this work, the original numerical models were modified accordingly to accept the uncertain model inputs and parameters as inputs from a file, and model simulations were performed at the sampled values of the inputs. The outputs of the simulations were then analyzed by ``binning'' them into a set of intervals; a bin number associated with an interval contains the number of outputs realizations that belong to that interval, from all the model simulations performed. The ratio of the bin number to the total number of samples, gives the probability that the output belongs to that interval. The ratio of this probability to the interval size gives the probability density at the mid-point of the interval. The probability density is then plotted against the range of values that the output can belong to, as presented in the case studies. Further, FORTRAN routines to calculate the mean, variance, and higher moments were also implemented as part of the SRSM; these were based on equations presented in Section 3.5.

The Latin Hypercube Sampling (LHS) method was developed and implemented from basic principles, mainly because of the lack of availability of a standardized routine readily usable on Unix operating system. The EPA approved routine for the LHS, developed by Iman et al. [108], is developed for the VAX/VMS operating system, and other commercial software such as Crystal ball [46] and @RISK [218], for the PCs have LHS routines embedded in the modeling platform (typically spreadsheet based), and thus cannot be used in conjunction with stand-alone numerical codes. The main aspects of the LHS implementation developed here are presented in the following.

In order to generate samples from input distributions, the range of probable values for each uncertain input parameter is divided into $M$ segments of equal probability. Thus, the whole parameter space, consisting of $N$parameters, is partitioned into $M^N$ cells, with each cell corresponding to equal probability regions from the input sampling space. For example, for the case of 3 input parameters and 4 segments, the parameter space is divided into $4\times4\times4$ cells. The next step involves the generation of $M$ samples from $M^N$ cells. This is accomplished as follows: first, a random sample is generated, and its cell number is calculated. The cell number indicates the segment number the sample belongs to, with respect to each of the parameters. For example, a cell number (2,1,2) indicates that the sample lies in the segment 2 with respect to first parameter, segment 1 with respect to second parameter, and segment 2 with respect to third parameter. At each successive step, a random sample is generated, and is accepted only if it does not agree with any previous sample on any of the segment numbers. Thus, no two samples have any input corresponding to the same segment. The advantage of this approach is that the random samples are generated from all the ranges of possible values, thus providing samples from the tails of the probability distributions (i.e., regions that have very low probabilities of occurrence, but correspond to extreme values). Further, for every $M$ samples, each parameter is sampled from each of its $M$ subranges.

The LHS FORTRAN subroutine developed here requires the following inputs: (a) the total number of samples to be generated, $N_{\mbox{tot}}$, (b) the number of random inputs, $N$, and (c) the input probability distributions, consisting of the distribution type (a number), and the distribution parameters (dependent on the distribution type). Further, the number of ranges the parameter range of each input should be divided ($M$) is also required as an input. The routine generates a set of $M$ sample points in each iteration, till $N_{\mbox{tot}}$ samples are generated. A typical value of $M$ used in the case studies is about 50. A much higher value results in too small a range for each sampled input, requiring a large number of trial samples before an acceptable sample is generated, whereas much smaller values do not adequately represent all ranges of possible values, especially the tails of the distributions. The routine, along with all other routines is included in the CDROM accompanying this document.

  

Figure 3.3: Schematic depiction of the steps involved in the application of the Stochastic Response Surface Method

3.8.2 Implementation of the Stochastic Response Surface Method

The steps involved in the application of the SRSM are presented in Figure 3.3. The implementation of the SRSM is based on the following routines:

• A symbolic calculation routine in Maple [94], to generate polynomial chaos expansions of arbitrary order, consisting of an arbitrary number of parameters.
• A FORTRAN subroutine to generate polynomial chaos expansions of up to fifth order for an arbitrary number of parameters. The extension to higher orders is straightforward, but not pursued because in all the case studies considered, the third order expansion was the highest required. This routine is developed using the results from the Maple routine.
• A FORTRAN subroutine to generate the sample points based on the input distributions and the order of the polynomial chaos expansion selected.
• A set of FORTRAN utilities to sample distributions in the form of the srvs, to calculate the corresponding transformed random variables, and to calculate the numerical values of the polynomial chaos expansions.
• A FORTRAN subroutine to perform collocation or regression (based on the scenario) on a set of sample points and the corresponding model outputs. These are used to estimate the coefficients of a polynomial chaos expansion.
• A FORTRAN subroutine to generate the probability density functions from the coefficients of a polynomial chaos expansion in srvs, by obtaining a large number of samples of srvs, and calculating the statistics of the corresponding outputs.
• FORTRAN subroutines to calculate the moments, percentiles, and other statistical properties of both the inputs and outputs.