MCDB-R: Risk Analysis in the Database Subi Arumugam1 Fei Xu3,∗ 1 University of Florida Gainesville, FL, USA {sa2,[email protected] Ravi Jampani1 Christopher Jermaine2 2 Rice University Houston, TX, USA {lp6,[email protected] ABSTRACT Enterprises often need to assess and manage the risk arising from uncertainty in their data. Such uncertainty is typically modeled as a probability distribution over the uncertain data values, specified by means of a complex (often predictive) stochastic model. The probability distribution over data values leads to a probability distribution over database query results, and risk assessment amounts to exploration of the upper or lower tail of a query-result distribution. In this paper, we extend the Monte Carlo Database System to efficiently obtain a set of samples from the tail of a query-result distribution by adapting recent “Gibbs cloning” ideas from the simulation literature to a database setting. 1. INTRODUCTION In the face of regulatory processes such as Basel II and Solvency 2, enterprises are becoming increasingly concerned with managing and assessing the credit, financial, engineering, and operational risk arising from uncertain data [16]. Examples of uncertain data include future values of financial assets, customer order quantities under hypothetical price changes, and transportation times for future shipments under alternative shipping schemes. Data uncertainty is usually modeled as a probability distribution over possible data values, and such probabilities are often specified via complex stochastic models. E.g., we might specify the foregoing uncertain financial, retail, and logistics data sets using Euler approximations to stochastic differential equations, Bayesian demand models, and stochastic network models, respectively. This specification in turn leads to the representation of data uncertainty as a probability distribution over possible database instances (“possible DBs,” also called “possible worlds”). Running an aggregation query such as “select sum of sales” over uncertain data therefore does not yield a single, deterministic result for total sales; rather, there is a probability distribution over all possible query results. In this setting, risk assessment typically corresponds to computing interesting properties of the upper or lower tails of the queryresult distribution; for example, computing the probability of a ∗Fei Xu performed this work while at U. Florida. Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Articles from this volume were presented at The 36th International Conference on Very Large Data Bases, September 13-17, 2010, Singapore. Proceedings of the VLDB Endowment, Vol. 3, No. 1 Copyright 2010 VLDB Endowment 2150-8097/10/09... $ 10.00. 3 Microsoft Corporation Redmond, WA, USA [email protected] Luis L. Perez2 Peter J. Haas4 4 IBM Research - Almaden San Jose, CA, USA [email protected] large investment loss. The problem is made more complex by the fact that the tails of the distribution are hard to specify a priori, and so risk analysis frequently focuses on the inverse problem of determining an extreme quantile, e.g., determining the value γ such that there is a 0.1% probability of seeing a loss of γ or more. Such a “value at risk” can be viewed as a probabilistic worst-case scenario, and might be used to specify precisely where the “upper tail” of the loss distribution starts. Proceeding further, one might be interested in computing a more sophisticated “coherent” risk measure such as “expected shortfall” [16], which in this example is defined as the expected total loss, given that this loss exceeds γ. More generally, the entire conditional distribution of the loss—given that the loss exceeds γ—might be of interest. MCDB In prior work, the Monte Carlo Database System (MCDB) of Jampani et al. [13] was designed for flexible exploration of queryresult distributions under arbitrary SQL queries and a broad range of complex, user-defined stochastic models. MCDB uses (possibly user-defined) “variable generation” (VG) functions to pseudorandomly generate instances of each uncertain data value in a database, yielding a sample from the possible-DB distribution. Repeating this process multiple times (i.e., executing multiple “Monte Carlo repetitions”) generates a set of independent and identically distributed (i.i.d.) samples from this distribution. Given an SQL aggregation query of interest, MCDB executes the query on each sampled DB instance, thereby generating i.i.d. samples from the query-result distribution. MCDB uses Monte Carlo techniques to estimate interesting features of the query-result distribution—the expected value, variance, and quantiles of the query answer—along with probabilistic error bounds on the estimates. Importantly, a VG function takes as input one or more parameter tables (ordinary relations) that control the function’s behavior, and produces as output a table containing one or more correlated data values. For n Monte Carlo repetitions, MCDB does not actually materialize an uncertain database n times; the costs for such a naive approach would be exorbitant. Instead, MCDB executes a query plan only once over a set of tuple bundles rather than over ordinary tuples—no matter how many Monte Carlo repetitions are required—often leading to significant time savings. A tuple bundle encapsulates the instantiations of a tuple over a set of generated DB instances, and carries along the pseudorandom number seeds used by the VG functions to instantiate the uncertain data values. Limitations of MCDB for Risk Analysis Unfortunately, naive Monte Carlo, as implemented in the original MCDB prototype, is not the best tool for exploring the tails of a query-result distribution. Consider the query SELECT SUM(loss) AS totalLoss FROM t, where t.loss is an uncertain attribute, perhaps representing a future financial loss. Suppose that the query-result distribution is normally distributed, with a mean of $10 million and a standard deviation of $1 million. Suppose that interest centers on the upper tail of this distribution that corresponds to totalLoss values of $15 million or more. On average, roughly 3.5 million Monte Carlo repetitions are required before such an extremely high loss is observed even once. If our goal is simply to estimate the area of the tail (i.e., the probability of seeing a totalLoss value of $15 million or more), then 130 billion repetitions are required to estimate the desired probability to within ±1% with a confidence of 95%. If the user instead wants to define the tail indirectly starting at the value γ, where γ is the 0.999 quantile of the totalLoss distribution, then standard quantile-estimation techniques [19, Sec. 2.6] require roughly ten million Monte Carlo repetitions to estimate γ to within ±1% with a confidence of 95%, Our Contributions In this paper, we describe an enhanced version of MCDB for risk analysis called MCDB-R. MCDB-R inherits the key strengths of MCDB. Chief among these is generality. MCDB can handle almost any stochastic model, as well as database operations such as aggregation, grouping, and so forth. However, unlike MCDB, MCDB-R supports efficient risk analysis. This presents several technical challenges: efficiently locating a tail of interest and efficiently sampling from the tail, all in the presence of black-box VG functions and complex SQL queries. Our contributions are as follows. 1. We adapt “cloning” and Gibbs-sampling ideas from the rareevent simulation literature to develop a statistical method for both estimating a user-specified quantile on a query-result distribution and then generating a set of samples from the tail. 2. We show how the tail-sampling method can be integrated into MCDB’s current tuple-bundle processing infrastructure. 3. In the Appendix, we provide guidance on both setting the sampling parameters for MCDB-R and identifying situations in which MCDB-R will be effective. Related Work The original MCDB system was partially inspired by several efforts to push analytic capabilities into the database or “close to the data,” such as MauveDB [10], which was oriented toward sensor data and integrated smoothing and interpolation methods. These ideas have been further developed in subsequent systems such as FunctionDB [21] and the SciDB project [20]. A recent discussion of analytics in the database can also be found in [6], and there have been a number of efforts to push statistical analyses into Map-Reduce data processing environments; see, e.g., [5, 12]. MCDB is also related to “probabilistic databases” (PrDBs) [7, 8]. Uncertainty in a PrDB is typically represented not by general stochastic models as in MCDB, but rather by explicit probabilities on alternative tuple values. MCDB can capture the functionality of a PrDB, but will not be as efficient or precise in cases where a PrDB can compute answers exactly and efficiently. Overall, MCDB is quite different from PrDBs in terms of motivation, design, and application space. The recent PIP system [14] combines PrDB and Monte Carlo techniques, yielding superior performance for certain MCDB-style queries with simple VG-function parameterizations. The current paper is also related to the general problem of “conditioning” a probabilistic database, as discussed in Koch and Olteanu [15]. The work in [15] focuses on exact probability calculations in the setting of PrDBs, whereas the current paper emphasizes Monte-Carlo-based approximations. 2. SPECIFYING QUERIES IN MCDB-R In this section, we indicate what MCDB-R looks like to a user, via a simple example. Suppose that we want to assess potential total financial loss over a specified subset of our clients. The steps in the analysis are to (1) define a stochastic loss model, (2) execute a SUM query over the uncertain loss values specified by the model, and then (3) report pertinent features of the total-loss distribution. For step (1), the user defines an uncertain loss table, just as in MCDB. Suppose for simplicity that we “model” the loss for a customer as being distributed according to a normal distribution with variance 1 and a mean that is customer specific. For this simple model, we use the built-in Normal VG function that generates normally distributed random numbers. The data needed to parameterize this function is stored on disk in a “parameter table” (an ordinary SQL table) called means(CID,m), which stores the mean loss for each customer. The uncertain table Losses(CID,val) that we wish to query is not stored on disk; only the schema is stored. The schema is specified by the following SQL-like statement: CREATE TABLE Losses (CID,val) AS FOR EACH CID IN means WITH myVal AS Normal(VALUES(m, 1.0)) SELECT CID, myVal.* FROM myVal This statement is identical to a standard SQL CREATE TABLE statement, except for the FOR EACH clause. The statement gives a recipe for constructing a sample instance of Losses. Scan through the CIDs in the means table. For each CID, create a one-row, onecolumn table myVal by invoking the Normal VG function, parameterized with the mean loss for the CID and the default variance value of 1.0. Form a row of the DB instance by joining myVal with the CID from the outer FOR EACH loop, as specified in the final SELECT clause. See [13] for more details on schema specification. For steps (2) and (3), we write a query as follows: SELECT SUM(val) as totalLoss FROM Losses WHERE CID < 10010 WITH RESULTDISTRIBUTION MONTECARLO(100) DOMAIN totalLoss >= QUANTILE(0.99) FREQUENCYTABLE totalLoss The first three lines correspond to a standard SQL (aggregation) query. Since this query is being executed over uncertain data, there is a probability distribution over the query result. The DOMAIN clause indicates that we want to condition the query-result distribution: we discard all possible query results not in the specified domain and then we renormalize the query-result distribution so that the total probability over the domain is 1. In our example, the domain coincides with the upper tail of the loss distribution corresponding to the highest 1% of losses. The MONTECARLO keyword indicates that we will use an approximate (conditioned) query-result distribution estimated from 100 Monte Carlo samples. The remaining lines of the WITH clause specify the statistical features of the query-result distribution that are to be computed. The FREQUENCYTABLE keyword causes creation of an ordinary table FTABLE(totalLoss,FRAC) in which the first column contains the distinct values of totalLoss observed over the 100 tail samples, and FRAC is the fraction of Monte Carlo samples in which each value was observed. By running the query SELECT MIN(totalLoss) FROM FTABLE we obtain an estimate of the lower tail-boundary, i.e., the extreme quantile of interest, which is accurate when the number l of Monte Carlo samples is large.1 Similarly, by executing a query such as 1 As discussed in the next section, our tail-sampling algorithm also produces an estimate of the extreme quantile, but the two estimates will be almost identical for large l. Algorithm 1 Systematic Gibbs sampler 1: Inputs: 2: X (0) : initial random element of X r 3: k: number of Gibbs updating steps 4: Output: 5: X (k) : updated value of X (0) 6: 7: G IBBS(X (0) , k): 8: x ← X (0) // Initialize 9: for j ← 1 to k do 10: // Perform one systematic updating step: X (j−1) → X (j) 11: for i ← 1 to r do 12: xi ← G EN C OND(x−i , i) // Generate from h∗i ( · |x−i ) 13: end for 14: end for 15: Return x // = X (k) SELECT SUM(totalLoss * FRAC) FROM FTABLE we can estimate the expected shortfall, i.e., the expected loss, given that the loss was among the highest 1% of losses. 3. MONTE CARLO BACKGROUND We next review results from the Monte Carlo literature and show how they are relevant to the problem of estimating an extreme quantile and efficiently sampling from the tail defined by the quantile. The key methods that we adapt are Gibbs sampling and cloningbased methods for rare-event simulation. 3.1 Gibbs Sampling Gibbs sampling [11] is a technique that is designed for generating samples from a high-dimensional probability distribution function that is known only up to a normalizing constant. For simplicity of notation, we describe the method in the setting of simple discrete random variables; the technique extends straightforwardly to continuous and mixed random variables. Let X = (X1 , X2 , . . . , Xr ) be an r-dimensional random vector taking values in a discrete set X r and distributed according to h, so that h(x1 , . . . , xr ) = P (X1 = x1 , . . . , Xr = xr ) for x = (x1 , . . . , xr ) ∈ X r . A key requirement of Gibbs sampling is that we can efficiently generate samples from the conditional distributions h∗i (u|v) = P (Xi = u | Xj = vj for j 6= i), where 1 ≤ i ≤ r. According to this definition, h∗i is the marginal distribution of Xi , given the values of Xj for j 6= i. For a vector x ∈ X r , we set x−i = (x1 , x2 , . . . , xi−1 , xi+1 , . . . , xr ) and then we can define h and h∗i more concisely as h(x) = P (X = x) and h∗i (u|v) = P (Xi = u | X−i = v) for u ∈ X and v ∈ X r−1 . (0) (0) Given an initial value X (0) = (X1 , . . . , Xr ), the “systematic” Gibbs sampler in Algorithm 1 generates a sequence X (0) , . . . , X (k) of random vectors. The function G EN C OND(x−i , i) invoked in Step 12 generates a sample from h∗i ( · |x−i ). Since each new sample is generated recursively from the previous sample, the sequence {X (j) } forms a Markov chain. The Gibbs sampler is a special case of a Markov Chain Monte Carlo (MCMC) sampling technique. In our database setting, the Gibbs sampler is more appealing that other MCMC methods because it is extremely generic in nature, requiring no special knowledge about the distribution h. Crucially, if the initial sample X (0) is generated from h, then the chain will be stationary in that every subsequent sample will be distributed according to h [1, Th. XIII.5.1]. Although the samples are not statistically independent, under mild regularity conditions Algorithm 2 Rejection algorithm for example Gibbs sampler 1: G EN C OND(v, i): 2: repeat 3: Generate u according to hi 4: until Q(u ⊕i v) ≥ c 5: Return u the random vectors X (0) and X (k) become increasingly independent as k increases.2 This “convergence to independence” is usually exponentially fast, so that k need not be very large. (In our experiments taking k = 1 sufficed.) Thus, if two chains {X (j) } and {Y (j) } start from the same state X (0) = Y (0) but the respective Gibbs updates are performed independently for the chains, then X (k) and Y (k) will become approximately independent as k increases, again typically at an exponentially fast rate. We illustrate the method further with an example that is pertinent to our tail-sampling problem. Suppose that r is large and that the components of X are statistically independent, so that h(x) = Q r i=1 hi (xi ), where hi (u) = P (Xi = u) is the marginal distribution of Xi for 1 ≤ i ≤ r. We assume that it is straightforward to generate samples from each hi . Let Q be a real-valued function defined on X r , and suppose that for some value c we want to generate samples from the conditional distribution h(x; c) = P (X = x | Q(X) ≥ c) = h(x)I Q(x) ≥ c /pc . Here I(A) = 1 if event APoccurs and I(A) = 0 otherwise, and pc = P (Q(X) ≥ c) = x∈X r :Q(x)≥c h(x). If the cardinality of X is very large, then pc is extremely hard to compute, so that h(x; c) is known only up to the constant of proportionality pc , and direct generation of samples from h(x; c) is nontrivial. We can, however, apply the Gibbs sampler. Note that, by the independence of the components of X, we have h∗i (u|v) = P Xi = u | Q(Xi ⊕i v) ≥ c . Here u ⊕i v denotes the vector (v1 , . . . , vi−1 , u, vi , . . . , vk−1 ). Thus, samples from h∗i in Step 12 can be generated by a rejection algorithm (Algorithm 2). To illustrate one especially simple case, Psuppose that Q(x) = x1 + · · · + xr . Then Q(u ⊕i v) = u + j vj , so that h∗i (u|v) = P (Xi = u | P Xi ≥ c − j vj ). An efficient implementation of Algorithm 2 would operate on an existing Gibbs iterateP x ∈ X r by subtracting − xi from Q(x), thereby computing q = j6=i xj , and then generating samples u from hi until u + q − ≥ c. In the following, we write G IBBS(X (0) , k, c) and G EN C OND(v, i, c) to indicate the dependence of these two algorithms on the parameter c. To see the relevance of the above example to MCDB-R, consider a database consisting of a single table R having a single, random attribute A and exactly r tuples in all possible worlds (e.g., like the Losses table in Section 2 with the CID attribute dropped). Suppose that we interpret Xi above as the random variable corresponding to the value of the ith tuple (1 ≤ i ≤ r) and Q as an aggregation query over the values in R. Then the foregoing discussion shows that if we have somehow obtained a DB instance D(0) sampled from h(x; c)—i.e., an instance that is “large” in the sense that Q(D(0) ) ≥ c—we can run the Gibbs sampler with the rejection algorithm for k steps to obtain another random instance D(k) that is (approximately) independent of D(0) and also ‘large”—i.e., also distributed according to h(x; c). That is, once we have obtained a DB instance corresponding to an upper tail of the query-result 2 More precisely, the Markov chain associated with the Gibbs sampler is usually “irreducible” and “aperiodic” in an appropriate sense, and the chain “mixes” at an exponential rate; see [3]. distribution, we can generate additional, independent instances that also yield query results lying in the tail. In the general MCDB-R setting, we can still interpret X as a random vector in which each component corresponds to an uncertain value in the database, but the components of X decompose into mutually independent blocks, where the variables within a block are dependent and are all generated via a call to a specified VG function. Moreover, the number of variables in a block—i.e., the number of values returned by a VG function—may vary over DB instances. Although formal discussion of Gibbs sampling at this level of generality would lead to very cumbersome notation, the basic idea is the same as above. 3.2 Rare-Event Simulation The foregoing discussion shows that Gibbs sampling can be used to generate essentially independent samples from a tail of the queryanswer distribution, starting from a “large” DB instance. However, such a starter instance can be hard to obtain, even when the quantile that defines the tail is known a priori. To address this challenging problem, we adapt ideas from the literature on “rare-event simulation”—see, e.g., [17] for an overview. The primary goal of rareevent simulation is to estimate the probabilities of very infrequent events; quantile estimation has received very little attention. The two main techniques for rare-event simulation are importance sampling (IS) and cloning, also called splitting. The idea in IS is to sample from a modified possible-DB distribution (by suitably modifying the VG functions) such that extreme DB instances are generated with high frequency. Due to numerical instability, however, IS is known to behave badly for high dimensional problems [2]. In our setting, the effective dimension typically is roughly equal to the number of tuples in a relation, rendering IS unusable. We therefore look to cloning techniques. The cloning approach has been of interest to the simulation community for a long time, but only recent work [2, 4, 18] has considered the type of “static” Monte Carlo simulations that are of interest in our setting. These algorithms start with a set of n randomly generated “particles.” Each particle gives rise to a Markov chain that is obtained by repeatedly applying a random MCMC update (“perturbation”) to the particle. The set of particles is “enriched” by deleting “non-extreme” particles and cloning “extreme” particles. The “Gibbs cloning” framework of Rubinstein [18] provides a general formulation of these ideas that we adapt for our purposes. 3.3 Adapting the Gibbs Cloner Adapting and specializing the Gibbs cloning framework to our setting—the application to quantiles is novel—yields the basic procedure of Algorithm 3. In the algorithm, we represent a random database D as a vector of random variables, as per our prior discussion. (We treat each deterministic data value c as a random variable that is equal to c with probability 1.) In the algorithm, we maintain a set S of DB instances over a sequence of bootstrapping steps, so called because we are “bootstrapping” our way out to the tail. At the start of the ith step, S contains ni instances, and the algorithm proceeds by (1) purging all but the top 100pi % “elite” (most extreme) instances, (2) cloning the elite instances to increase the size of S up to ni+1 , and then (3) perturbing each instance in S using the Gibbs sampler, in order to re-establish (approximate) independence while ensuring that all instances yield query results that lie in the current tail. The function C LONE(S, n) operates by simply duplicating each DB instance in S approximately n/|S| times. The Appendix describes how to choose the algorithm parameters to efficiently control the relative difference between the desired upper-tail probability p and the actual probability of the upper tail Algorithm 3 Basic tail-sampling algorithm 1: Inputs: 2: p: target upper-tail probability 3: l: desired number of tail samples 4: Outputs: 5: γ̂: estimate of (1 − p)-quantile 6: S: set of l samples from h( · ; γ̂) 7: Parameters 8: k: number of Gibbs updating steps 9: m: number of bootstrapping steps 10: n1 , n2 , . . . , nm : intermediate sample sizes 11: p1 , p2 , . . . , pm : intermediate tail probabilities 12: 13: // Initialize 14: Generate databases D(1) , . . . , D(n1 ) i.i.d. according to h( · ) 15: S ← { D(1) , D(2) , . . . , D(n1 ) } 16: nm+1 ← l 17: // Execute m bootstrapping steps 18: for i ← 1 to m do 19: γ̂i ← the (pi |S|)-largest element of { Q(D) : D ∈ S } 20: Discard all elements D ∈ S with Q(D) < γ̂i 21: S ← C LONE(S, ni+1 ) 22: for D ∈ S do // Gibbs-update step 23: D ← G IBBS(D, k, γ̂i ) // G IBBS defined as in Sec. 3.1 24: end for 25: end for 26: return γ̂ = γ̂m and S defined by the quantile estimator γ̂. In particular, we show that this relative error is minimized by setting ni = n and pi = p1/m for 1 ≤ i ≤ m. Here n = N/m, where N is the total number of samples over all bootstrapping steps; increasing N increases cost but improves accuracy. We assume these values throughout. After the ith purge, the smallest query-result value for the retained elite samples is an estimate of γi , defined to be the (1 − pi/m )-quantile of the query-result distribution; thus the γi ’s are increasing and γm corresponds to the extreme quantile of interest. At each step i, we estimate a 1−p1/m quantile of Fi−1 , the conditional distribution of the query result Q(D), given that Q(D) ≥ γi−1 . (Here we take γ0 = −∞.) For typical values of, say, p = 0.001 and m = 4, we see that even though we ultimately need to estimate an extreme 0.999-quantile, at each step we merely need to estimate a 0.82-quantile, a much easier problem. 4. IMPLEMENTATION OVERVIEW In the following sections, we describe how we implemented tail sampling (or TS for short) as in Algorithm 3. The goal is to preserve, to the extent possible, the computational efficiency of the MCDB tuple-bundle processing mechanism as discussed previously. 4.1 The Gibbs Looper and Data Streams In the MCDB-R setting, we call the iterative process in Algorithm 3 the Gibbs Looper. There are two key variables that track the progress of the this algorithm: cutoff and curQuantile. In terms of our previous notation, cutoff is the γi value that defines the current tail and, after i iterations, curQuantile = pi/m is the probability that a query-result instance falls in the current tail. At each iteration of the Gibbs Looper, each of the n versions (i.e., realized instances) of the database is randomly perturbed by the Gibbs sampler to form n new DB versions. Initially, a cloned database has the same contents as its parent, but as each DB ver- sion is perturbed by its Gibbs sampler, parent and child will drift apart. To “fuel” the perturbation, the Gibbs sampler needs access to a stream of random data. There is a data stream associated with every uncertain data value (or correlated set of uncertain data values) in the database. E.g., consider the random Losses table from Section 2. The tuple bundle corresponding to a given customer is initially enhanced with a seed for use by the pseudorandom number generator (PRNG). Repeated execution of the Normal VG function, parameterized with the customer’s mean loss value m, produces a stream of realized loss values for the customer. In MCDB, the first n elements of the stream simply correspond to the customer’s losses in n Monte Carlo repetitions, i.e., in n generated DB instances. In MCDB-R, because of the rejection algorithm (Algorithm 2) used in the Gibbs sampler, there is a more complicated mapping from DB instance to position in the stream, but the Gibbs sampler nonetheless goes to the stream whenever it needs a loss value for the customer. We will often refer to a stream of data via the PRNG seed that was used to produce it. In general, a VG function may produce a table of correlated data values when called, so an “element” of a stream may actually comprise a set of values. Also, a given PRNG seed may occur in multiple tuples—due to a 1-to-m join of the random table to another table—or multiple times in a tuple bundle due to a self-join. 4.2 An Example We illustrate the Gibbs Looper in more detail using the total-loss query of Section 2 as an example; for brevity, we drop the CID column from the Losses and means tables, as well as the selection predicate on CID from the query itself. Suppose that there are three customers and that the corresponding three rows in the means table have values 3.0, 4.0, and 5.0, respectively. Also suppose that p = 1/32 and n = 4, so that we want to generate a sample of four total losses that are each in the top 3.125% of the total-loss (i.e., query-result) distribution. We will use m = 5 iterations of the Gibbs Looper to generate the samples. First, the underlying VG function is used to generate three streams of random data, as depicted in Fig. 1. The first four values produced by each VG function are assigned, in sequence, to each of the n = 4 DB versions, as in Fig. 1(a). After this initialization, iteration i = 1 of the Gibbs Looper begins. First, the two DB versions with final SUM values in the lowest 100(1 − p1/m )% percentile (50% in our example) are discarded, and the remaining two elite DB versions are cloned. “Cloning” is accomplished by duplicating the assignment of random values from each stream of random data to the various DB versions, as depicted in Fig. 1(b). The current value for cutoff at this point is set to be 12.07, since this is the aggregate value associated with the least-extreme remaining DB version. cutoff also serves as an estimate for the point at which there is only a 100p1/m % (= 50%) chance of seeing a larger value for the answer to the SUM query. Since two of the DB versions have been duplicated, it is necessary to perturb all of the versions so that they are sufficiently differentiated. In our example, the Gibbs sampler begins the perturbation by attempting to replace the value 3.26 that is associated with the first stream of loss values (the stream associated with mean value 3.0). This is done by considering the next unassigned random value associated with mean 3.0—indicated in bold in Fig. 1(b). Unfortunately, this value (which is 3.24) cannot be used to replace 3.26, because it would result in an overall aggregate value of 12.05, which does not meet or exceed the current value of cutoff = 12.07. Thus, the 3.24 is rejected. Next, the 5.13 is considered. Since the 5.13 would result in an overall aggregate value of 13.94 (which exceeds cutoff), the 5.13 is accepted. The procedure is repeated using the other two streams. The 4.59 associated with mean.val = 4.0 is replaced with 4.55 (acceptable since this decreases the aggregate value to 13.90, which exceeds the cutoff of 12.07), and the 4.22 associated with mean.val = 5.0 is replaced with 4.29, resulting in the state of Fig. 1(c). DB version two is updated in the same way. The value 3.26 is replaced with 4.07 (acceptable, since the new aggregate value is 12.88). An attempt is made to replace the 4.59 associated with mean.val = 4.0 using 3.68. However, this would reduce the overall aggregate for DB version two to 11.97, which does not meet or exceed the current value for cutoff. So the 3.68 is rejected, and 5.16 (the next unassigned value associated with mean.val = 4.0) is chosen. Finally, the 4.22 is replaced with 5.77, resulting in Fig. 1(d). After DB versions three and four are updated, the result is Fig. 1(e). Then iteration i = 2 begins, and cloning the top two DB versions gives us Fig. 1(f). This whole process is repeated until the end of iteration i = m. 4.3 Implementing A Gibbs Looper In practice, there are many possible ways to implement a Gibbs Looper in the presence of complicated database operations such as joins. The simplest would maintain n different random versions of the database. Then, in each iteration of the Gibbs looper (where an “iteration” is the process of perturbing all n DB versions) and for each DB version, we would undertake the following steps. For each stream of random data, we would (1) replace the current stream element that is being used in the DB version with a new stream element, (2) run the query over the modified DB version, (3) accept the new stream element if the new query result exceeds cutoff and otherwise try steps (1) and (2) again. The drawback of such an implementation is that it requires running the underlying query many, many times over many, many random databases. For example, imagine that we will maintain n = 100 DB versions, each requiring a modest one million streams of random data. Imagine that the Gibbs Looper needs ten iterations to complete and that we need to reject (on average) ten random values before one is accepted. In this case, we will need to run an entire query plan 100 × 106 × 10 × 10 = 1010 different times to complete the TS procedure, which is obviously unacceptable. We therefore use MCDB’s tuple-bundle trick of running the underlying query only once, no matter how many DB versions are needed. The extensions to MCDB’s tuple-bundle mechanism are many and intricate; we now sketch the key ideas. 5. GIBBS TUPLE BUNDLES As with MCDB, query execution in MCDB-R is organized around the idea of bundling together the different versions of a tuple over the generated DB instances. To execute a query, MCDB-R runs a plan that pushes a stream of instantiated tuple bundles into a special GibbsLooper operation. These bundles must have all of the information necessary to implement the process illustrated in Fig. 1, i.e., given a DB version and a stream of random data, the Gibbs Looper must be able to determine which random value from the stream is assigned to the DB version, and then efficiently determine the effect of an update to the value on the final query result. This requires that MCDB tuple bundles must be augmented with a notion of “lineage” that links each random value in the database to the stream from which it came. Gibbs tuples therefore differ from MCDB tuple bundles. Instead of containing a PRNG seed, a Gibbs tuple contains a pointer to a “tail-sampling seed”, or TS-seed for short. This object, described in Section 6 below, augments the basic PRNG seed with “bookkeeping” information needed by the Gibbs Looper. In an MCDB bun- SUM = 14.17 (e) 3.08 3.91 7.64 SUM = 14.63 DB Version 1 DB Version 2 means.val = 5.0 assigned to? assigned to? means.val = 4.0 DB Version 3 DB Version 1 SUM = 15.00 4.07 5.16 5.77 DB Version 2 assigned to? means.val = 5.0 assigned to? means.val = 4.0 assigned to? means.val = 3.0 DB Version 1 DB Version 2 means.val = 5.0 assigned to? SUM = 13.11 4.07 5.16 5.77 2.59 3.18 4.87 3.26 4.59 4.22 2.23 3.10 6.09 SUM = 15.00 4.56 3.96 5.65 3.24 4.55 4.29 5.13 3.68 5.77 V12 3.08 4.07 V12 5.16 V12 4.04 3.91 2.24 5.56 7.64 V34 0.13 3.91 V34 4.54 7.64 3.51 4.01 3.69 2.44 4.02 5.85 SUM = 14.63 3.08 V34 3.84 6.58 3.41 4.13 5.36 3.08 4.25 4.55 5.53 1.22 5.16 3.77 3.91 2.82 3.17 6.18 7.64 DB Version 3 3.51 5.56 4.04 SUM = 14.17 (c) DB Version 4 V4 V2 V3 V4 SUM = 15.00 V1 V2 V3 V4 SUM = 13.97 DB Version 4 DB Version 4 V3 V1 4.87 4.22 6.09 5.65 4.29 5.77 4.04 7.64 4.54 3.69 5.85 6.58 5.36 5.53 3.77 6.18 assigned to? means.val = 3.0 DB Version 1 DB Version 2 assigned to? means.val = 5.0 means.val = 4.0 assigned to? assigned to? V1 V2 3.18 4.59 3.10 3.96 4.55 3.68 5.16 5.56 3.91 4.01 4.02 3.84 4.13 4.55 5.16 3.17 4.07 5.16 5.77 5.13 4.55 4.29 3.26 4.59 2.59 3.18 4.87 4.22 3.26 V2 4.59 V2 4.22 V2 2.23 3.10 6.09 4.56 V34 3.96 V34 5.65 V34 SUM = 12.07 3.24 4.55 V1 4.29 V1 5.13 V1 3.68 5.77 4.56 4.07 5.16 4.04 3.96 2.24 5.56 7.64 0.13 3.91 4.54 5.65 3.51 4.01 3.69 2.44 4.02 5.85 SUM = 14.17 3.08 3.84 6.58 3.41 4.13 5.36 4.56 4.25 4.55 5.53 1.22 5.16 3.77 3.96 2.82 3.17 6.18 5.65 SUM = 13.97 DB Version 3 2.59 3.26 2.23 4.56 3.24 5.13 4.07 2.24 0.13 3.51 2.44 3.08 3.41 4.25 1.22 2.82 5.13 4.55 4.29 DB Version 4 DB Version 4 DB Version 3 2.59 3.18 4.87 3.26 4.59 4.22 2.23 3.10 6.09 4.56 V34 3.96 V34 5.65 V34 SUM = 15.00 3.24 4.55 V1 4.29 V1 5.13 V1 3.68 5.77 V2 4.56 4.07 V2 5.16 V2 4.04 3.96 2.24 5.56 7.64 0.13 3.91 4.54 5.65 3.51 4.01 3.69 2.44 4.02 5.85 SUM = 14.17 3.08 3.84 6.58 3.41 4.13 5.36 4.56 4.25 4.55 5.53 1.22 5.16 3.77 3.96 2.82 3.17 6.18 5.65 assigned to? 4.07 5.16 5.77 means.val = 4.0 SUM = 13.97 DB Version 2 means.val = 5.0 assigned to? means.val = 4.0 assigned to? assigned to? means.val = 3.0 5.13 4.55 4.29 3.26 4.59 4.22 SUM = 14.17 (b) DB Version 1 SUM = 14.17 means.val = 3.0 4.56 3.96 5.65 SUM = 12.07 DB Version 3 SUM = 11.42 (a) (d) means.val = 3.0 DB Version 1 SUM = 12.07 2.23 3.10 6.09 3.26 4.59 4.22 2.59 3.18 4.87 3.26 V12 4.59 V12 4.22 V12 2.23 3.10 6.09 4.56 V34 3.96 V34 5.65 V34 SUM = 12.07 3.24 4.55 4.29 5.13 3.68 5.77 4.56 4.07 5.16 4.04 3.96 2.24 5.56 7.64 0.13 3.91 4.54 5.65 3.51 4.01 3.69 2.44 4.02 5.85 SUM = 14.17 3.08 3.84 6.58 3.41 4.13 5.36 4.56 4.25 4.55 5.53 1.22 5.16 3.77 3.96 2.82 3.17 6.18 5.65 assigned to? V1 V2 V3 V4 DB Version 2 means.val = 5.0 4.87 4.22 6.09 5.65 4.29 5.77 4.04 7.64 4.54 3.69 5.85 6.58 5.36 5.53 3.77 6.18 3.26 4.59 4.22 DB Version 3 V1 V2 V3 V4 assigned to? means.val = 4.0 3.18 4.59 3.10 3.96 4.55 3.68 5.16 5.56 3.91 4.01 4.02 3.84 4.13 4.55 5.16 3.17 SUM = 10.64 DB Version 4 V1 V2 V3 V4 assigned to? assigned to? means.val = 3.0 2.59 3.26 2.23 4.56 3.24 5.13 4.07 2.24 0.13 3.51 2.44 3.08 3.41 4.25 1.22 2.82 2.59 3.18 4.87 (f) SUM = 14.63 Figure 1: Running the tail-sampling process to obtain four DB instances dle, PRNG seeds may be discarded if they are no longer needed. In MCDB-R, a TS-seed handle can never be discarded, because each random value must always point to its associated PRNG (i.e., associated stream). Another important difference is that, whereas an instantiated MCDB tuple bundle contains exactly n stream elements per random value—one for each Monte Carlo repetition—an instantiated Gibbs bundle contains many more elements, because the Gibbs Looper can use up many stream elements during the Gibbs perturbation process.3 Indeed, the Gibbs Looper may run out of data during processing and need to re-execute part of the query plan (see Section 9). The number of stream elements to instantiate in a Gibbs tuple is chosen to trade off the cost of carrying lots of data throughout the query plan with the cost of re-executing part of the query if the data runs out. For an example of how Gibbs tuples are pushed through a query 3 Because instantiated values are often needed to prepare tuple bundles for the Gibbs Looper—see the example below—they cannot be simply generated on the fly during Gibbs-Looper operation. plan, consider the following query, which determines a company’s total salary “inversion,” resulting from certain employees who earn more than their managers: SELECT SUM(emp2.sal - emp1.sal) FROM emp AS emp1, emp AS emp2, sup WHERE sup.boss = emp1.eid AND emp1.sal < 90K AND sup.peon = emp2.eid AND emp2.sal > 25K AND emp2.sal > emp1.sal Imagine that the attribute emp.sal is random, and we want three samples from the tail of the query-result distribution. The resulting query plan and execution of this query is illustrated in Fig. 2. This plan operates over Gibbs tuples rather than ordinary tuples. An instantiated Gibbs tuple can have one or more arrays of random values; every array of random values is associated with exactly one PRNG seed. A Gibbs tuple can also have one or more isPres attributes. An array of isPres values is created when a selection predicate is applied to a random attribute, and indicates for each DB instance whether or not the predicate is satisfied. If the predicate emp(sal,eid) ??? Joe ??? Sue ??? Ann ??? Jim ??? Sid temp(sal,eid) (Seed1) (Seed2) (Seed3) (Seed4) (Seed5) ??? Joe ??? Sue ??? Ann ??? Jim ??? Sid Seed sup(boss,peon) Sue Joe Jim Sue Jim Ann Sid Jim eid=boss temp(sal,eid1,eid2) (Seed2) (Seed4) ??? Jim Sue ??? Sue Joe (Seed4) (Seed5) ??? Jim Ann ??? Sid Jim temp(sal,eid1,sal2,eid2) (Seed2) (Seed1) (Seed4) (Seed2) (Seed4) (Seed3) (Seed5) (Seed4) ??? Sue ??? Joe ??? Jim ??? Sue ??? Jim ??? Ann ??? Sid ??? Jim eid=peon Instantiatesa1 temp(sal,eid1,sal2,eid2) (Seed2) (Seed1) (Seed4) (Seed2) 24K Sue ??? Joe 77K Jim ??? Sue 76K 25K 72K 26K 79K 22K 70K 23K (Seed4) (Seed3) (Seed5) (Seed4) 77K Jim ??? Ann 92K Sid ??? Jim 76K 94K 72K 99K 79K 97K 70K 91K temp(sal,eid1,isPres1,sal2,eid2) (Seed2) (Seed2) (Seed4) (Seed1) 24K Sue Yes ??? Joe 77K Jim Yes ??? Sue 25K Yes 76K Yes 26K Yes 72K Yes 22K Yes 79K Yes 23K Yes 70K Yes σ sal1<90K σ sal2>25K temp(sal,eid1,isPres1,sal2,eid2) (Seed2) (Seed4) (Seed1) 24K Sue Yes 28K Joe Yes 77K Jim Yes 25K 76K Yes Yes 26K Yes 26K 72K Yes Yes 25K No 22K 79K Yes Yes 21K No 23K 70K Yes Yes 27K Yes (Seed2) 24K Sue 25K 26K 22K 23K Yes Yes Yes Yes Yes (Seed3) ??? Ann Instantiatesal2 temp(sal,eid1,isPres1,sal2,eid2) (Seed2) (Seed2) (Seed4) (Seed1) 24K Sue Yes 28K Joe 77K Jim Yes 24K Sue 25K Yes 26K 76K Yes 25K 26K Yes 25K 72K Yes 26K 22K Yes 21K 79K Yes 22K 23K Yes 27K 70K Yes 23K No No Yes No No (Seed4) 77K Jim 76K 72K 79K 70K (Seed4) 77K Jim 76K 72K 79K 70K Yes Yes Yes Yes Yes (Seed3) 45K Ann Yes Yes 39K Yes 43K Yes 47K Yes 44K (Seed4) 77K Jim 76K 72K 79K 70K Yes Yes Yes Yes Yes GibbsLooper( 0.001, sal2>sal1, sum(sal2-sal1)) (Seed3) 45K Ann 39K 43K 47K 44K res 14K 11K 10K Figure 2: Gibbs-tuple processing over a multi-table query plan is not satisfied in any DB instance, then the entire Gibbs tuple is dropped. (Unlike MCDB, we cannot simply record isPres values for the entire tuple, since attribute values change individually during Gibbs sampling.) Finally, note that the query plan makes use of the MCDB-R Seed and Instantiate operations. The former operation attaches the handle for a TS-seed to each Gibbs tuple, and the latter operation uses a PRNG seed to attach a subsequence from a stream of random values to the Gibbs tuple. 6. TAIL-SAMPLING SEEDS A TS-seed contains both the PRNG seed needed to instantiate a stream of random values and the information needed to map the current value of an attribute in a DB version to a position in the stream. This mapping is needed to compute a query-result value for the current DB version. E.g., in Fig. 2, the query plan produces a set of three Gibbs tuples that are piped into the GibbsLooper operation. Suppose that we have the assignment (seed1, 1), (seed2, 1), (seed3, 2), (seed4, 2) for a DB version. That is, Joe’s salary is $28K (the first value in his stream), Ann’s salary is $39K (the second value in her stream), and so forth. Then the query result for this DB version is ($28K - $24K) = $4K, since Joe makes $4K more than Sue and no other employee makes more than their boss. In MCDB-R, a TS-seed contains (1) a TS-seed identifier, (2) the actual PRNG seed used to produce a stream of random data, (3) the range of stream values currently materialized and present within the Gibbs tuples, (4) the last random value in that range that has previously been assigned to any DB version for this TS-seed, and (5) the random value currently assigned to each DB version for this TS-seed. Item (5) is the mapping described above, and item (4) is needed during Gibbs rejection sampling to obtain a new random value. Item (3) is needed to quickly determine whether the Gibbs Looper has run out of data and, if so, the position in the stream from which to start replenishing. A call to the Seed operator both adds a seed handle to a Gibbs tuple and creates the actual TS-seed data structure for the tuple, which is written to disk. 7. THE GIBBS LOOPER The Appendix describes the GibbsLooper operator in detail. A key feature is that it switches the inner and outer for loops of Algorithm 3. Rather than perturb DB versions one at a time, looping through the random data values (i.e., the TS-seeds), the algorithm perturbs data values one at a time, looping through the DB versions, thereby amortizing expensive data scans. Specifically, GibbsLooper inserts the Gibbs tuples into a diskbased priority queue; the initial sort key for a Gibbs tuple is the smallest TS-seed handle in the tuple. GibbsLooper then iterates through the overall set of TS-seed handles in increasing order. For each TS-seed handle, it loads the TS-seed—and all “associated” Gibbs tuples that contain at least one handle for the TS-seed—into memory. Then, for the first DB version, it consults the TS-seed to determine the first unused stream value to try in the rejection algorithm. It uses the associated Gibbs tuples to compute the changed value of the query result for the DB version due to the new stream value, and then checks if the new query-result value still exceeds cutoff. If so, then GibbsLooper accepts the new stream value and moves on to update the next DB version. If not, it tries again until an appropriate stream value is found. Once all DB versions have been updated, it reinserts all of the Gibbs tuples into the priority queue—first updating the sort key for each reinserted tuple to be the next largest TS-seed handle in the tuple—and moves on to the next overall TS-seed handle. (Thus a Gibbs tuple will be processed multiple times if it contains multiple TS-seed handles.) GibbsLooper is efficient because it essentially merges Gibbs tuples in the disk-based priority queue with a sorted file containing all of the TS-seeds. When a TS-seed is loaded into memory, the priority queue is used to quickly obtain all of the Gibbs tuples that could possibly be affected by an update of the TS-seed. 8. JOINS ON RANDOM ATTRIBUTES One issue that demands more explanation is how GibbsLooper handles join predicates on random attributes. E.g., suppose R1 joins with R2 on a random attribute, and R2 joins with R3 on a deterministic attribute. Some t2 in R2 is modified when one of its random attributes is updated, so it joins with both a t1 from R1 and a t3 from R3 —the result is that a new tuple “pops” into existence. Fortunately, this case is easily handled by MCDB-R. When the underlying query plan is run (before GibbsLooper is ever invoked), original MCDB’s Split operation is run on any random attribute or attributes that are to be joined. This operation changes the attributes so that they are no longer random and transfers any non-determinism to the isPres attribute. E.g., consider a Gibbs tuple g with schema (fname, age) where attribute fname is constant, and attribute age is non-constant. Specifically, age = 20 for the first and third values in the stream of random data stored in the Gibbs tuple, and age = 21 for the second and fourth values: g = (Jane,(20,21,20,21),(Y,Y,Y,Y)), where the last nested vector contains the isPres values, and indicates that Jane appeared in all four Monte Carlo repetitions (though with varying ages). An application of the Split operation to t with Atts = {age} yields two Gibbs tuples: g1 = (Jane,20,(Y,N,Y,N)) and g2 = (Jane,21,(N,Y,N,Y)). Then, any join on the age attribute is a join on a deterministic attribute, and we will not have the problem of new Gibbs tuples “popping” into existence, since two sets of Gibbs tuples (each associated with one of Jane’s possible ages) would both be created and sent to GibbsLooper. Furthermore, only one of them can be valid at a time: if, e.g., Jane’s age is updated to be 21, then the mapping information within the TS-seed associated with Jane’s age is updated as well, causing any Gibbs tuple associated with an age of 20 to become invalid. 9. RE-RUNNING THE QUERY PLAN Since each Gibbs tuple contains only a finite amount of data, the Gibbs Looper may run out of data while perturbing a DB version. If so, MCDB-R discards all existing Gibbs tuples, empties the GibbsLooper priority queue, and then runs the entire query plan again to produce the required data. The newly-created Gibbs tuples that contain the data are piped into the GibbsLooper operation. During such a “replenishing” run, query-plan execution differs from the first time around in a couple of ways. First, the result of each deterministic part of the query plan is materialized and saved during the first run, avoiding the need to repeat these computations during replenishment. Second, the Instantiate operation never adds stream values to a Gibbs tuple that have already been processed; it only adds new or currently assigned values. 10. CONCLUSION MCDB-R extends the MCDB system to permit risk evaluation in the database. Space constraints have forced us to exclude some important material from the paper. Some very brief benchmarking results, given in the Appendix, provide evidence for the accuracy of the method and show a reduction in processing time from roughly 18 hours in MCDB to 11 minutes in MCDB-R. Acknowledgements This work was supported in part by the US DoE under grant no. DE-SC0001779 and the US NSF under grant no. 0915315, and by a Sloan Foundation Fellowship. 11. REFERENCES [1] S. Asmussen and P. W. Glynn. Stochastic Simulation: Algorithms and Analysis. Springer, 2007. [2] Z. I. Botev and D. P. Kroese. An efficient algorithm for rare-event probability estimation, combinatorial optimization, and counting. Methodol. Comput. Appl. Prob., 10:471–505, 2008. [3] R. C. Bradley. Basic properties of strong mixing conditions: A survey and some open questions. Probab. Surveys, 2:107–144, 2005. [4] F. Cérou, P. D. Moral, T. Furon, and A. Guyader. Rare event simulation for a static distribution. INRIA Research Report 6792, Rennes, France, 2009. [5] C.-T. Chu, S. K. Kim, Y.-A. Lin, Y. Yu, G. R. Bradski, A. Y. Ng, and K. Olukotun. Map-reduce for machine learning on multicore. In NIPS, pages 281–288, 2006. [6] J. Cohen, B. Dolan, M. Dunlap, J. M. Hellerstein, and C. Welton. MAD skills: New analysis practices for big data. PVLDB, 2(2):1481–1492, 2009. [7] N. N. Dalvi, C. Ré, and D. Suciu. Probabilistic databases: diamonds in the dirt. Commun. ACM, 52(7):86–94, 2009. [8] A. Das Sarma, O. Benjelloun, A. Y. Halevy, S. U. Nabar, and J. Widom. Representing uncertain data: models, properties, and algorithms. VLDB J., 18(5):989–1019, 2009. [9] H. A. David and H. N. Nagaraja. Order Statistics. Wiley, third edition, 2003. [10] A. Deshpande and S. Madden. MauveDB: supporting model-based user views in database systems. In SIGMOD, pages 73–84, 2006. [11] S. Geman and D. Geman. Stochastic relaxation, Gibbs distribution and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intelligence, 6(6):721–741, 1984. [12] S. Guha. RHIPE - R and Hadoop Integrated Processing Environment. http://ml.stat.purdue.edu/rhipe/. [13] R. Jampani, F. Xu, M. Wu, L. L. Perez, C. M. Jermaine, and P. J. Haas. MCDB: a Monte Carlo approach to managing uncertain data. In ACM SIGMOD, pages 687–700, 2008. [14] O. Kennedy and C. Koch. PIP: A database system for great and small expectations. In ICDE, pages 157–168, 2010. [15] C. Koch and D. Olteanu. Conditioning probabilistic databases. PVLDB, 1(1):313–325, 2008. [16] A. J. McNeil, R. Frey, and P. Embrechts. Quantitative Risk Management: Concepts, Techniques, and Tools. Princeton University Press, 2005. [17] G. Rubino and B. Tuffin, editors. Rare Event Simulation Using Monte Carlo. Wiley, 2009. [18] R. Rubinstein. The Gibbs cloner for combinatorial optimization, counting, and sampling. Methodol. Comput. Appl. Prob., 11(4):491–549, 2009. [19] R. J. Serfling. Approximation Theorems of Mathematical Statistics. Wiley, 1980. [20] M. Stonebraker, J. Becla, D. J. DeWitt, K.-T. Lim, D. Maier, O. Ratzesberger, and S. B. Zdonik. Requirements for science data bases and SciDB. In CIDR, page 26, 2009. [21] A. Thiagarajan and S. Madden. Querying continuous functions in a database system. In SIGMOD, pages 791–804, 2008. APPENDIX A. THE GIBBS LOOPER IN DETAIL A.1 Initialization Aside from n, m, and p (see Section 3.3), GibbsLooper takes as input five principal arguments: 1. The pipeline producing the input stream of Gibbs tuples. 2. The aggregate operation from whose tail distribution GibbsLooper will sample, e.g., sum(sal2-sal1) in Fig. 2. 3. The final selection predicate (e.g., sal2 > sal1 in Fig. 2) that must be applied to each tuple prior to inclusion of the tuple’s value in the aggregate. In particular, any selection predicate that references random attributes generated by more than one PRNG seed must be pulled up into GibbsLooper. Such a predicate cannot be applied lower in the query plan, because the Gibbs sampling will arbitrarily combine random values from the streams associated with different PRNG seeds, and it is not possible to know the outcome of the predicate before GibbsLooper runs its Gibbs sampler. 4. The grouping attributes (if any) for the query.4 5. A file containing all of the TS-seed objects, sorted on each object’s handle. Given these inputs, GibbsLooper initially processes the stream of incoming Gibbs tuples by computing the initial aggregate value for each of the n DB versions and inserting the tuples into a diskbased priority queue, sorted by the smallest TS-seed handle in each Gibbs tuple (as discussed in Section 7). The initial aggregate values are computed for the TS-seed mapping in which the ith value in each stream is mapped to the ith DB version for 1 ≤ i ≤ n. Next, GibbsLooper discards the “non-elite” DB versions whose aggregate totals are in the bottom 100(1 − p1/m )% by overwriting them with clones of the elite versions. The cutoff value is set equal to the smallest aggregate value. The overwriting of a nonelite DB version by a clone of an elite version is accomplished during a single read/write pass through the TS-seed file (see below). A.2 Performing the Perturbation At this point, the GibbsLooper is ready to begin perturbing the various DB versions using a Gibbs sampler. The manner in which the perturbation is accomplished is best illustrated using an example, given as Fig. 3. Fig. 3(a) begins with the three Gibbs tuples that were output from the query plan of Fig. 2, inserted into the GibbsLooper’s priority queue. (For simplicity, we drop the isPres attributes and look at only two DB instances.) The four TS-seeds are also depicted, as are the two DB versions that correspond to the initial TS-seed mapping (described above). The goal at this point to perturb these versions via Gibbs sampling, subject to the constraint that the total aggregate value must meet or exceed cutoff = $1K. Note that while the evolving DB versions are depicted in the figure, they are never actually materialized; they are completely determined, however, by the current state of the Gibbs tuples and the TS-seeds. To obtain two new DB versions with SUM values that meet or exceed cutoff, the first TS-seed handle, (seed1), is considered, and the lone associated Gibbs tuple is removed from the priority queue, as in Fig. 3(b). DB version one is updated by assigning the first unassigned random value associated with seed1 to the DB version, transforming (24K, Sue, 28K, Joe) into (24K, Sue, 25K, 4 Joe). This decreases the final aggregate result from $4 to $1K, since Joe’s salary exceeds Sue’s by $1K. Because the total salary inversion still meets the cutoff of $1K, the update is accepted. Next, we attempt to update DB version two. The next unassigned random value in the stream associated with seed1 is the fourth one, which transforms (25K, Sue, 26K, Joe) into (25K, Sue, 21K, Joe). This changes the final aggregate value to $0, and so this proposed update is rejected. Next we try the fifth random value in the stream. This gives us (25K, Sue, 27K, Joe), which increases the final aggregate value to $2K, and so it is accepted. We are now done with seed1, and so we reinsert the one Gibbs tuple associated with this seed back into the priority queue, after first updating its sort key to its next largest seed handle, (seed2); see Fig. 3(c). As shown in Fig. 3(d), the two Gibbs tuples associated with (seed2) are then removed from the priority queue. The assignments for seed2’s stream are updated in much the same way as for seed1; note, however, that seed2 influences two tuples in each DB version. Once the update of seed2 is complete and both of the associated Gibbs tuples have been reinserted back into the priority queue, we are left with Fig. 3(e). Note that the keys for these reinserted Gibbs tuples are now (seed4) and infinity. The value infinity is used because the random values associated with seed1 and seed2 have been each been updated—there is no reason to consider this Gibbs tuple again during this particular perturbation, so it is forced to the tail of the queue. As illustrated in Fig. 3(f), we now move to (seed3). The two DB instances are updated with respect to seed3; because seed3 control’s Ann’s salary and Ann has no impact on the query result, both updates are trivially accepted. The process continues with (seed4), which again has no effect on the final query result, at which point both DB versions are fully perturbed (Fig. 4(a)). Since all of the TS-seed handles have been processed, it is now time to clone the elite DB versions and overwrite the non-elite versions. This process is illustrated by the transition from Fig. 4(a) to Fig. 4(b). In a single read/write pass, the column in each TSseed that records the assignment for DB version two is simply copied to the column for version one. Now the first iteration of GibbsLooper has completed, and iteration two can begin. Grouping is handled by, in effect, treating a GROUP BY query over g groups as g separate, simultaneous queries, each with a selection predicate that limits the query to a specific group. B. WHEN WILL MCDB-R WORK BEST? Algorithm 1—and hence Algorithm 3—will work best when, for each most-likely extreme database D, the query result Q(D) is relatively insensitive to changes in the value of a single data value x. This scenario fails, e.g., for a SUM query in which data values X1 , . . . , Xr are i.i.d. according to a “subexponential” distribution such as lognormal or Pareto. Under such a heavy-tailed distribution, P (X1 + · · · + Xk > c) ≈ kP (X1 > c) for large c; i.e., with high probability, the query result for an extreme database D— viewed as a vector y of values—will lie in the upper tail because some individual yi is extremely large. After generating a new candidate value u for yi in Algorithm 2, the value of Q(u ⊕i y) will likely drop sharply, and many candidates will be required prior to acceptance in Line 4. Our approach is thus unlikely to be universally applicable. That said, it is still useful for a wide class of riskestimation problems, as evidenced by the extensive applied literature that pertains to extreme-value theory for light-tailed distributions and by the prevalence of aggregates, such as SUM and AVG, that have the insensitivity property under a light-tailed regime. C. SETTING PARAMETERS In this section, we address the key question of how to choose the parameters in Algorithm 3. Because the cost of a Gibbs-updating 1 1 1 1 5 5 5 5 2 2 2 2 1 1 1 1 DB Version 1 24K Sue 28K Joe 77K Jim 45K Ann 77K Jim 24K Sue SUM = 4K 2 2 2 2 seed1 62 seed2 27 seed3 91 seed4 14 DB Version 2 25K Sue 26K Joe 77K Jim 27K Ann 76K Jim 25K Sue (a) PRNG seed low iter. num. high iter. num. max. used iter. V1 iter. num. V2 iter. num. Head (seed3) 45K Ann 39K 43K 47K 44K Priority Queue (seed2) (seed4) (seed2) 77K Jim 24K Sue 76K 25K 72K 26K 79K 22K 70K 23K seed identifier PRNG seed low iter. num. high iter. num. max. used iter. V1 iter. num. V2 iter. num. Head seed identifier seed1 62 seed2 27 seed3 91 seed4 14 (seed3) (seed4) 77K Jim 76K 72K 79K 70K 1 1 1 1 5 5 5 5 5 2 2 2 3 1 1 1 5 2 2 2 (seed1) (seed2) (seed1) 24K Sue 28K Joe 25K 26K 26K 25K 22K 21K 27K 23K SUM = 2K 3 1 1 1 seed1 62 seed2 27 seed3 91 seed4 14 DB Version 2 25K Sue 27K Joe 77K Jim 39K Ann 76K Jim 25K Sue seed identifier PRNG seed low iter. num. high iter. num. max. used iter. V1 iter. num. V2 iter. num. 5 5 5 5 5 5 2 2 3 4 1 1 5 5 2 2 DB Version 1 SUM = 3K 23K Sue 27K Joe 77K Jim 39K Ann 76K Jim 23K Sue 5 5 2 2 3 4 1 1 5 5 2 2 (seed2) (seed4) 77K Jim 76K 72K 79K 70K seed1 62 seed2 27 seed3 91 seed4 14 1 1 1 1 5 5 5 5 5 5 4 2 3 4 3 1 5 5 4 2 (seed3) (seed4) 77K Jim 76K 72K 79K 70K DB Version 1 22K Sue 25K Joe 77K Jim 45K Ann 77K Jim 22K Sue SUM = 3K DB Version 2 (seed2) 24K Sue 25K 26K 22K 23K Priority Queue (seed4) (seed4) (seed2) 77K Jim 24K Sue 76K 25K 72K 26K 79K 22K 70K 23K 22K Sue 25K Joe 77K Jim 45K Ann 77K Jim 22K Sue DB Version 2 5 5 5 5 (seed2) (seed2) (seed1) 24K Sue 28K Joe 25K 26K 26K 25K 21K 22K 27K 23K (seed3) 45K Ann 39K 43K 47K 44K 23K Sue 27K Joe 77K Jim 39K Ann 76K Jim 23K Sue SUM = 4K (infinity) (seed2) (seed1) 24K Sue 28K Joe 25K 26K 26K 25K 22K 21K 23K 27K Tail (e) 1 1 1 1 (infinity) (seed2) (seed1) 24K Sue 28K Joe 25K 26K 26K 25K 22K 21K 23K 27K 1 1 1 1 (d) Tail seed1 62 seed2 27 seed3 91 seed4 14 Priority Queue (seed4) (seed4) (seed2) 77K Jim 24K Sue 76K 25K 72K 26K 79K 22K 70K 23K PRNG seed low iter. num. high iter. num. max. used iter. V1 iter. num. V2 iter. num. seed identifier SUM = 1K 5 2 2 2 (seed3) 45K Ann 39K 43K 47K 44K Head DB Version 1 SUM = 2K (seed3) (seed4) 77K Jim 76K 72K 79K 70K (seed3) (seed4) 77K Jim 76K 72K 79K 70K 24K Sue 25K Joe 77K Jim 45K Ann 77K Jim 24K Sue (c) Head (seed3) 45K Ann 39K 43K 47K 44K PRNG seed low iter. num. high iter. num. max. used iter. V1 iter. num. V2 iter. num. 5 2 2 2 (seed3) (seed4) 77K Jim 76K 72K 79K 70K seed identifier 5 5 5 5 Priority Queue Head Head seed identifier PRNG seed low iter. num. high iter. num. max. used iter. V1 iter. num. V2 iter. num. 1 1 1 1 SUM = 1K Tail seed1 62 seed2 27 seed3 91 seed4 14 DB Version 1 24K Sue 25K Joe 77K Jim 45K Ann 77K Jim 27K Sue DB Version 2 (b) Tail (seed2) (seed2) (seed1) 24K Sue 28K Joe 25K 26K 26K 25K 22K 21K 23K 27K (seed3) 45K Ann 39K 43K 47K 44K 25K Sue 27K Joe 77K Jim 27K Ann 76K Jim 25K Sue SUM = 1K Priority Queue (seed2) (seed4) (seed2) 77K Jim 24K Sue 76K 25K 72K 26K 79K 22K 70K 23K (seed3) (seed4) 77K Jim 76K 72K 79K 70K Tail Priority Queue (seed2) (seed4) (seed2) 77K Jim 24K Sue 76K 25K 72K 26K 79K 22K 70K 23K Tail (seed1) (seed2) (seed1) 24K Sue 28K Joe 25K 26K 26K 25K 22K 21K 23K 27K DB Version 1 (seed3) 45K Ann 39K 43K 47K 44K (f) 22K Sue 25K Joe 77K Jim 43K Ann 77K Jim 22K Sue SUM = 3K DB Version 2 23K Sue 27K Joe 76K Jim 47K Ann 76K Jim 23K Sue SUM = 4K SUM = 3K Figure 3: The perturbation process for the “salary inversion” query iteration is proportional to the data size and hence extremely expensive, we choose the number of iterations k to be as small as possible, subject to the requirement that, after the Gibbs-update step in line 23, the ni+1 databases in S are approximately statistically in- dependent. The MCMC literature provides a variety of procedures for determining k. Recall, however, that a sequence of Gibbs samples typically “forgets” its initial state exponentially fast, so that very small values of k suffice in practice; as mentioned previously, 1 1 1 1 5 5 5 5 5 5 4 4 3 4 3 3 5 5 4 4 (seed4) (seed4) 77K Jim 76K 72K 79K 70K (a) DB Version 1 (seed3) 45K Ann 39K 43K 47K 44K 22K Sue 25K Joe 72K Jim 43K Ann 72K Jim 22K Sue SUM = 3K DB Version 2 (seed2) 24K Sue 25K 26K 22K 23K 23K Sue 27K Joe 79K Jim 47K Ann 79K Jim 23K Sue (infinity) (seed4) (seed2) 77K Jim 24K Sue 76K 25K 72K 26K 79K 22K 70K 23K PRNG seed low iter. num. high iter. num. max. used iter. V1 iter. num. V2 iter. num. seed identifier PRNG seed low iter. num. high iter. num. max. used iter. V1 iter. num. V2 iter. num. seed identifier seed1 62 seed2 27 seed3 91 seed4 14 (infinity) (seed4) (seed3) 77K Jim 45K Ann 76K 39K 72K 43K 79K 47K 70K 44K seed1 62 seed2 27 seed3 91 seed4 14 1 1 1 1 5 5 5 5 5 5 4 4 5 5 4 4 5 5 4 4 (b) SUM = 4K (infinity) (seed2) (seed1) 24K Sue 28K Joe 25K 26K 26K 25K 22K 21K 23K 27K Tail (seed4) (seed4) 77K Jim 76K 72K 79K 70K Tail (infinity) (seed2) (seed1) 24K Sue 28K Joe 25K 26K 26K 25K 22K 21K 23K 27K DB Version 1 23K Sue 27K Joe 79K Jim 47K Ann 79K Jim 23K Sue SUM = 4K DB Version 2 23K Sue 27K Joe 79K Jim 47K Ann 79K Jim 23K Sue SUM = 4K Figure 4: Perturbation example, continued taking k = 1 worked well in experiments. Now suppose that k has been fixed and is sufficiently large so that Gibbs updating yields samples that can be viewed as independent. Also suppose that the total number N of samples over all bootstrapping steps has been fixed. We consider the problem of selecting m, n1 , . . . , nm and p1 , . . . , pm to minimize statistical error. One approach is to minimize the relative difference between the desired upper-tail probability p and the actual probability of the upper tail defined by the quantile estimator γ̂m . I.e., we solve the following optimization problem: minimize m,n1 ,...,nm ,p1 ,...,pm such that m X ni = N i=1 E and h F̄ (γ̂ ) − p 2 i 0 m p m Y pi = p (1) (2) i=1 pi ∈ [0, 1] and ni ∈ { 0, 1, . . . , N } , 1≤i≤m (3) Here F0 (x) = P Q(D) < x and F̄0 = 1 − F0 . Given a solution to the minimization problem (1)–(3), the overall parameterselection task will reduce to selection of N . We refer to the minimand in (1) as the mean-squared relative error (MSRE). Note that an alternative goal might be to control the relative error of γ̂m directly. For light-tailed query-result distributions, which are the class of primary interest in this paper, minimizing MSRE rather than minimizing the relative error of γ̂m typically leads to more stringent error control. E.g., suppose that, unknown to us, F0 is the standard normal distribution, and that our desired tail probability is p = 0.001, so that γ ≈ 3.090. If we allow the actual tail probability to deviate from p by no more than ±1%, then the corresponding estimate of γ will lie approximately in the range (3.087, 3.093), which represents a maximum error of about ±0.1%. In the following, we fix k and N as described above, and show how to select values of m, n1 , . . . , nm , and p1 , . . . , pm that approximately solve the minimization problem (1)–(3). We assume that F0 is continuous and strictly increasing. We first abstract the process by which Algorithm 3 computes γ̂m . For real numbers c1 , c2 , . . . , cm , set F1 (x; c1 ) = F0 (x) − F0 (c1 ) = P Q(D) < x | Q(D) ≥ c1 , 1 − F0 (c1 ) and for 2 ≤ i ≤ m recursively define Fi (x; c1 , . . . , ci ) = Fi−1 (x; ; c1 , . . . , ci−1 ) − Fi−1 (ci ; c1 , . . . , ci−1 ) 1 − Fi−1 (ci ; c1 , . . . , ci−1 ) for any sequence c1 , c2 , . . . , cm . An inductive argument establishes the identity Fi (x; c1 , . . . , ci ) = F1 (x; ci ) for 2 ≤ i ≤ m. For a set X1 , X2 , . . . , Xn of i.i.d. samples from a strictly increasing continuous distribution function F , the order statistics are denoted as X(1) < X(2) < · · · < X(n) . Setting ri = ni (1 − pi ) (rounded to the nearest integer), and assuming that k is sufficiently large as discussed previously, we can view Algorithm 3 as follows. The algorithm proceeds by generating n1 i.i.d. samples X0,1 , X0,2 , . . . , X0,n1 from F0 . Then, for i = 1, 2, . . . , m, the algorithm (1) sets γ̂i = Xi−1,(ri ) and (2) generates ni+1 samples Xi,1 , . . . , Xi,ni+1 from Fi ( · ; γ̂1 , . . . , γ̂i ) = F1 ( · ; γ̂i ). To analyze the behavior of F̄0 (γ̂m ), and hence the behavior of the MSRE, we employ the common tactic of reducing the analysis to the setting of uniform random variables. Specifically, let { Ui,j : 0 ≤ i ≤ m, 1 ≤ j ≤ ni+1 } be a collection of i.i.d. random variables uniformly distributed on [0, 1]. Consider a procedure that sets V0,j = U0,j for 1 ≤ j ≤ n1 , and then, for i = 1, 2, . . . , m, sets α̂i = Vi−1,(ri ) and Vi,j = (1 − α̂i )Ui,j + α̂i for 1 ≤ j ≤ ni+1 . Using a standard result [1, p. 37] on uniform random variables together with an inductive argument, it can be shown that, for 1 ≤ i ≤ m, (1) F0−1 (Vi−1,j ) has the same distribution as Xi−1,j for 1 ≤ j ≤ ni and (2) F0−1 (α̂i ) has the same distribution as γ̂i . In particular, F̄0 (γ̂m ) is distributed as 1 − α̂m . An easy argument shows that 1 − α̂i = Zi (1Q− α̂i−1 ), where Zi = 1−Ui−1,(ri ) , which implies that 1− α̂m = m i=1 Zi . By construction, the random variables {Ui−1,(ri ) }, and hence the random variables {Zi }, are mutually independent. Moreover, each Ui,(ri ) has a Beta(ri , ni − ri + 1) distribution [9]. Fixing m ∈ { 1, 2, . . . , N }, η = (n1 , . . . , nm ), and π = (p1 , . . . , pm ) satisfying (2)–(3), a calculation using both the independence of the Zi ’s and standard properties of the Beta distribution shows that MSRE = u(η, π, m), where u(η, π, m) = h1 (η, π, m) h2 (η, π, m)p−2 − 2p−1 + 1. Q Here hc (η, π, m) = m i=1 (ni pi + c)/(ni + c) for c = 1, 2. We now consider how to choose η, π, and m so as to approximately minimize u, subject to (2)–(3). It follows from the results given below that p ≤ hc (η, π, m) ≤ 1 for c = 1, 2 and feasible choices of the parameters. Moreover, h1 (η, π, m) ≈ h2 (η, π, m) for m ∈ [1, N ]. We can establish the following result. T HEOREM 1. For c ∈ { 1, 2 }, hc (η, π, m) is minimized by setting m∗c = min { m ≥ 1 : gm (N, p, c) < gm+1 (N, p, c) }, ηc∗ = ∗ ∗ (N/m∗c , . . . , N/m∗c ), and πc∗ = (p1/mc , . . . , p1/mc ). To summarize, we can approximately minimize the MSRE by (1) computing m∗1 and m∗2 as indicated in the theorem, (2) choosing θ = arg minc∈{1,2} u(ηc∗ , πc∗ , m∗c ) and setting m∗ = m∗θ , and (3) ∗ setting ni = N/m∗ and pi = p1/m for 1 ≤ i ≤ m∗ . In practice, ∗ ∗ it is often the case that m1 = m2 . Thus we need only choose N . For an MSRE target value , we can select N by numerically finding min { N : w(N ) ≤ }, where w(N ) = gm∗ (N, p, 1) gm∗ (N, p, 2)p−2 − 2p−1 + 1. Note that limN →∞ w(N ) = 0, which implies that our quantile estimator converges in mean square to the true value as the total number of DB versions increases. D. BENCHMARKING We have implemented a C++ prototype of MCDB-R with the MCDB code base as a starting point, including all of the features described in the paper. As with our MCDB implementation, MCDBR does not yet have an optimizer or SQL compiler; instead, we use an MCDB-specific language to specify a query plan directly. We now describe a small benchmark of our prototype. Consider the following random version of the TPC-H orders table: CREATE TABLE random ord (o orderkey, o yr, o tot) AS FOR EACH o IN (SELECT * FROM orders) WITH VAL AS Normal (VALUES(o mean, o var)) SELECT o.o orderkey, year(o.o orderdate), v.value FROM VAL v This version of orders attaches a normally-distributed random “loss” value to each tuple (we use a mean and variance of one). Now, consider the following query: SELECT SUM(val) as totalLoss FROM random ord, lineitem WHERE o orderkey = l orderkey AND (o yr = ’1994’ OR o yr = ’1995’) Because the sum of a set of normal variables is also normal, we know that the result of this query must itself have a normal distribution with mean and variance computed as follows: SELECT SUM(grpsize * o mean) AS mean, SUM(grpsize * grpsize * o var) AS var FROM (SELECT o mean, o var, COUNT(*) AS grpsize FROM orders, lineitem WHERE year(o orderdate) in (’1994’, ’1995’) AND o orderkey = l orderkey GROUP BY o orderkey, o mean, o var) Thus, if we run the former query using adding the following: WITH RESULTDISTRIBUTION MONTECARLO(100) DOMAIN totalLoss >= QUANTILE(0.999) then we can use the result of the latter query to test the accuracy of the resulting tail sample. Using the TPC-H database (scalefactor = 10) on an 8-core server, we ran the former query using 1 0.9 Cumulative Probability for the range of p and N values encountered in practice. Thus, to a good approximation, the problem of minimizing ulooks like the problem of minimizing v(x) = x (x/p2 ) − (2/p) such that x ∈ [p, 1]. It is easy to see that v is strictly increasing on [p, 1], and so we want to make x as small as possible. That is, we want to make h1 and h2 as small as possible. To this end, set (N/m)p1/m + c m gm (N, p, c) = (N/m) + c 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 4.975E+05 5.025E+05 5.075E+05 Query Result Figure 5: Observed empirical CDFs versus analytic CDF GibbsLooper parameters m = 5, p1/m = 0.25, N = 500, and l = 100. This should (if everything is correct) compute a set of 100 samples from the 1.0 − (0.25)5 = 0.99902 ≈ 0.999 quantile of the query-result distribution. Each TS-seed is used to produce 1000 random values initially. The first three GibbsLooper iterations require 156, 124, and 134 seconds to complete. At that point, the first TS-seed runs out of random data, and the query plan is run again to fuel the last two iterations. These iterations take 122 seconds and 115 seconds to run, for a total time of around eleven minutes to complete the tail sampling process and obtain 100 database instances in the upper 0.999-quantile. By comparison, MCDB would require about 18 hours for this task. Next we tested the accuracy of MCDB-R. We re-ran the above query, but changed how the underlying database is generated to make the example more interesting (and difficult). We modified random ord so that it is generated by selecting 100,000 tuples from the orders parameter table. The mean (resp., variance) of each normal used to generate random ord was itself generated by sampling from an inverse gamma-distribution with shape 3 and scale 1 (resp., shape 3 and scale 0.5). The lineitem table was generated so that one million of the tuples join with some tuple from random ord, and the rest find no mate. The probability that one of those tuples from lineitem will mate with the first of the 100,000 tuples in random ord is 2 × 10−5 − 10−10 . The probability that the tuple will mate with the ith tuple in random ord is equal to the probability that it will mate with the (i − 1)th tuple, minus 2 × (10−5 − 10−10 )/(105 − 1). We ran MCDB-R 20 times on this query using GibbsLooper parameters m = 5, p1/m = 0.25, N = 1000, and l = 100. Figure 5 shows the 20 empirical tail CDFs that we observed (each based on 100 samples). The true CDF for the query-result distribution (calculated analytically as described above) appears as a thick black line on the left side of the plot. The true tail CDF corresponding to the 0.99902 quantile of the query-result distribution appears as a thick black line on the right. It is clear from the plot that the 20 empirical tail CDFs cluster closely around the true tail CDF. We also obtained 20 estimates of the 0.99902 quantile by recording the value of the minimum tail sample in each of the 20 runs. The mean value of these 20 estimates was 5.0728e5, whereas the true quantile value was 5.0738e5. The empirical standard error of the estimates was 265. To put this in perspective, the middle 99% of the query-result distribution has a width of approximately 2503; i.e., by generating only 1000 intermediate database instances (efficiently, using Gibbs tuples), MCDB-R is able to “walk” out to the upper 0.99902 quartile and incur a standard error of only 10%.

- Similar pages