Using Supervised Learning to Improve Monte Carlo Integral Estimation
Monte Carlo (MC) techniques are often used to estimate integrals of a multivariate function using randomly generated samples of the function. In light of the increasing interest in uncertainty quantification and robust design applications in aerospace engineering, the calculation of expected values of such functions (e.g. performance measures) becomes important. However, MC techniques often suffer from high variance and slow convergence as the number of samples increases. In this paper we present Stacked Monte Carlo (StackMC), a new method for post-processing an existing set of MC samples to improve the associated integral estimate. StackMC is based on the supervised learning techniques of fitting functions and cross validation. It should reduce the variance of any type of Monte Carlo integral estimate (simple sampling, importance sampling, quasi-Monte Carlo, MCMC, etc.) without adding bias. We report on an extensive set of experiments confirming that the StackMC estimate of an integral is more accurate than both the associated unprocessed Monte Carlo estimate and an estimate based on a functional fit to the MC samples. These experiments run over a wide variety of integration spaces, numbers of sample points, dimensions, and fitting functions. In particular, we apply StackMC in estimating the expected value of the fuel burn metric of future commercial aircraft and in estimating sonic boom loudness measures. We compare the efficiency of StackMC with that of more standard methods and show that for negligible additional computational cost significant increases in accuracy are gained.
|Bias of M|
|Set of samples of|
|Subset of Samples in the Testing Set|
|Subset of Samples in the Training Set|
|True Expected Value of|
|Estimate of from M|
|Function Value at the sample|
|Prediction of Fit at|
|Expected Value of|
|Number of Folds|
|Likelihood of the Expected Value of the Fold|
|An Estimator of|
|Number of Samples in the Testing Set|
|Number of Samples|
|Number of Samples in the Training Set|
|Probability Distribution of|
|Alternate Sample Distribution|
|Fitting Algorithm for|
|Variance of M|
|Free Parameter of StackMC|
|Free Parameter in the Fitting Algorithm|
A system optimized only for best performance at nominal operating conditions can see severe performance degradation at even slightly off-design conditions. Aerospace systems rarely operate at precisely the design condition and a robust design approach dictates trading some performance at the nominal condition for improved performance over a wide range of input conditions. However, certain input conditions will occur rarely or never, and adding robustness for these conditions will degrade average system performance. Consequently, we are often interested in optimizing the expected performance of a system rather than the performance at any specific operating point.
Formally, we are interested in an integral of the form
where is the (potentially multivariate) objective function to be optimized, and is the probability density function (or probability distribution, in which case (1) is actually a summation) from which is generated.
Evaluating (1) is non-trivial, especially when the objective function is high dimensional and/or expensive to evaluate. When standard quadrature techniques (such as Simpson’s rule) become intractable, Monte Carlo (MC) methods are powerful techniques which have become the standard for integral estimation. However, MC techniques also have their drawbacks; the required computational expense may be prohibitive in most situations. The question is if a technique exists to significantly lower the cost of estimating (1). This paper describes such a technique which combines MC with machine learning and statistical techniques and can lead to significant computational savings over traditional methods.
We begin this paper with a brief review of integral estimation techniques including Monte Carlo, fitting algorithms, and stacking. We then introduce a new technique we call Stacked Monte Carlo (StackMC), which uses stacking to reduce the estimation error of any Monte Carlo method. Finally, we apply StackMC to a number of analytic example problems and to two problems from the aerospace literature. We show that StackMC can significantly reduce estimation error and never has higher estimation error than the best of Monte Carlo and the chosen fitting algorithm.
2 Integral Estimation Techniques
2.1 Estimation Error
When using any method returning as an estimate of (1), there are two sources of estimation error: bias () and variance (). Bias is the expected difference between the method and the truth: , where the expectation is over all data sets. A method whose average over many runs gives the correct answer has zero bias, whereas a method that estimates high or low on average has non-zero bias. As an example, the Euler equations have significant bias in their estimate of viscous drag. Variance is a measure of how much varies between different runs of : . If is deterministic, it has zero variance because every run returns the same answer, whereas if is stochastic multiple runs have different outputs leading to a non-zero variance. The total expected squared error is the combination of these two factors
Any method seeking to reduce estimation error must keep both of sources of error low.
2.2 Monte Carlo Techniques
In “Simple Sampling” Monte Carlo (MC), a set of samples, , is generated randomly according to . The estimate of the expected value of the function based on is the average of the function values of the samples
where refers to the generated sample. Simple Monte Carlo has two extremely useful properties. First, MC is guaranteed to be unbiased and thus, on average, will return the correct answer. Second, MC is guaranteed to converge to the correct answer at a rate of independent of the number of dimensions. That is, for any problem, increasing the number of samples by a factor of one hundred decreases the expected error by a factor of ten. As Simple Monte Carlo has zero bias, all of the error comes from the variance of Monte Carlo runs. This variance is due a different kind of variance concerning itself; fluctuations in the Monte Carlo estimate are directly related to fluctuations in . As a result, the smaller the variance of , the smaller the expected error of Monte Carlo.
Many different sample generation techniques exist to help reduce the variance of Monte Carlo. Importance sampling  generates samples from an alternate distribution (instead of the true distribution ) and estimates integral (1) as a weighted average of sample values
Importance sampling is often used when the tails of a distribution have a measurable effect on , but occur very infrequently. Using Simple MC, several million samples are needed to accurately capture a once-in-a-million event, but by using importance sampling unlikely outcomes can be made to occur more frequently reducing the total number of samples needed.
Quasi-Monte Carlo (QMC) techniques reduce variance by choosing sample locations more carefully. Sample points generated from simple Monte Carlo will inevitably cluster in some locations and leave other locations void of samples. QMC methods are usually deterministic and reduce variance by spreading points evenly throughout the sample space. Two such methods are the scrambled Halton sequence  and the Sobol sequence . While often effective, due to deterministic sample generation they are not guaranteed to be unbiased. It is also difficult to generate points from an arbitrary ; most QMC algorithms generate sample points from a uniform distribution over the unit hypercube.
2.3 Supervised Learning
A second class of techniques for estimating integrals from data seeks to use the data samples more efficiently through the use of a fitting algorithm and supervised learning techniques. The fitting algorithm uses data samples to make a “fit” to the data, and the integral estimate is the integral of the fit. Examples of fitting algorithms include splines, high-order polynomials and Fourier expansions; even Simpson’s rule computes the quadrature by fitting a piecewise quadratic polynomial to the data samples. Fits incorporate the spatial distribution of sample points and thus often give more accurate estimates of than MC. However, using a fitting algorithm can induce bias and may or may not exhibit convergence to the correct answer as the number of points increases. Additionally, when the number of data points is “too small” many fitting algorithms exhibit higher variance than MC, and, as a result, using a fitting algorithm can be worse than not using one, especially with low numbers of sample points and in high-dimensional spaces. It is usually impossible to know a priori how many points is “too small” making it difficult to know when to use a fitting algorithm.
Additionally, it is difficult to know if a fit to the data samples is an accurate representation of at values of not in the data set. Many fitting algorithms exhibit a property known as “overfitting” where a fit to the data is very accurate at locations that are among the data samples, but is very inaccurate at locations not among the data samples. As a result, one cannot use the data samples themselves to both make a fit and to evaluate its accuracy. The standard technique for addressing this issue is known as cross-validation, where a certain percentage of the data is used to make the fit and the rest of the data is used to evaluate its performance. Usually this process is repeated several times, and the best performing fit is used as the output of the fitting algorithm (a.k.a. winner-takes-all).
Stacking is a technique first introduced by Wolpert  which improves upon cross-validation by combining the results of different fits. Wolpert describes stacking in his paper as “ a strategy … which is more sophisticated than winner-takes-all. Loosely speaking, this strategy is to combine the [fits] rather than choose one amongst them. … [Stacking] can be viewed as a more flexible version of nonparametric statistics techniques like cross-validation, … therefore it can be argued that for almost any generalization or classification problem … one should use [stacking] rather than any single generalizer by itself.” Interested readers can see stacking applied to improving regression , probability density estimation , confidence interval estimation , and optimization . Stacking recently gained fame as being a major part of the first and second place algorithms in the Netflix competition .
In addition to combining the predictions of multiple fitting algorithms, stacking can also be used to improve the prediction of a single fitting algorithm. In this variant of stacking, one again partitions the full set of data into a subset used to make fits and the remainder of the data, and a subset . However now one does not use to determine how to combine the predictions of multiple fitting algorithms that were all run on . Instead one uses to correct the prediction of a single algorithm that was run on . This variant of stacking is closest to the algorithm we use below. For a comparison of stacking to other generalizers, see Clarke.
3 Stacked Monte Carlo
The main idea of Stacked Monte Carlo (StackMC) is to incorporate the variance reduction potential of a fitting algorithm while avoiding introducing bias and additional variance from poor fits. It makes no assumptions about the function nor the sample generation method, and therefore StackMC can be used to augment nearly any Monte Carlo application.
Let us assume that we have a function that is a reasonable (though not necessarily perfect) fit to . Equation (1) can be re-written as
Since is assumed to be a reasonable fit, for a properly chosen , has lower variance than alone (see Fig. 1), and so a MC estimate of (3) has lower expected error than an estimate of (1). In fact, it can be shown  that the optimal value of to minimize expected error is
where and are the standard deviations of and respectively, and is the correlation between and . Intuitively, if is a good fit to , (and correspondingly ) will be high and will be trusted as a good estimate for . If is a poor fit to , and will both be low and will be downplayed. Looking at it from a different perspective, (3) takes the expected value estimated from , and corrects it with a term estimating the bias of . Either way, by using (3) the error should be lower than using either Monte Carlo or the fitting function alone. Since the expected value of the fit is both added and subtracted, (4) remains an unbiased estimate of (1). Thus, (4) incorporates a fitter while remaining unbiased, and allows us to emphasize good fits while de-emphasizing poor ones.
The obvious question is how to obtain and find from data samples. The first step is to pick a functional form for a fitting algorithm , such as a polynomial with unknown coefficients. can be nearly anything, the only restrictions are that it can make predictions at new values. For now we also require that can be analytically integrated over the volume, i.e. we can calculate
analytically (see further discussion later in the paper). By comparing the output of a fit to the true values, an estimate for the “goodness” of the fit (and by extention ) could be obtained, and (4) could be applied. However, doing this directly would cause over-fitting and would lead to an inaccurate estimate of .
Over-fitting can be mitigated by using a technique known as -fold cross-validation . The data samples are randomly partitioned into testing sets, , which are mutually exclusive and contain the same number of samples (differing slightly if is not an integer). Training sets, , contain all of the data samples not in . We call the number of samples in and the number of samples in (such that ). One fit per fold, , is created using only the samples in (for a total of fitters). The samples in a training set are called “held in” points because they are used to generate the fit, whereas points in the corresponding testing set are “held out” points because the fit is generated without using these samples. Each data sample is in a held out data set once and in the held in data set times. Using a prediction of the function values for the points in is made ( ). Standard cross-validation evaluates the accuracy of each of the fits, and chooses the best fit to use as a single .
Instead, we adopt the approach of stacking and use (3) to get an estimate of from each of the fitters, and we can use the held out data points as an estimate of the integral term.
The mean of these estimates is taken as the final prediction
We can also use the predictions at the held out points to estimate . is estimated from the variance of the sample function values themselves, from the variance of the predictions, and from the correlation between the predictions and the truth.
Finally, we plug into (9)
One final correction is needed to complete StackMC. It is possible that some or all of the differ greatly from but is calculated high because of good predictions at held out data points. If left alone, this would cause StackMC to return a low-quality prediction on some data sets.
However, the error in the mean (EIM) of the Monte Carlo samples can be used as a second metric to evaluate the goodness of the fit. EIM is defined as
and represents the uncertainty in the mean of a set of Monte Carlo samples. Specifically, it gives us a likelihood bound on based on and . We can find a “likelihood” for each fold by calculating
The higher , the less likely it is that , and the more likely it is that the fitter is bad and should be ignored. A heuristic is set that if , , i.e. all of the fits are ignored entirely and the Monte Carlo average is used. If is set too low, then too many reasonable fits are ignored, and if is set too high, then too many unreasonable fits are kept. A value of was chosen based on experimentation; it was clear that setting as low as or as high as were inferior.
Finally, a note about the bias of StackMC. Because the expected value of the fit is being both added and subtracted, it is true that (3) (and by extension, (7) and (9)) are also unbiased for a fixed . However, using the held-out data to both set and estimate the integral can introduce bias. In practice, we have found that this is only a problem for very small numbers of sample points where the output fit of the fitting algorithm changes significantly depending on the held-in samples.
3.1 Generalized Stacked Monte Carlo
The discussion above was for the case where the samples are generated according to a known , this is only a special case of a class of estimation scenarios. While the overall methodology remains approximately the same for these other scenarios, some specifics (such as the calculation of ) change with the choice for and the sample generation method.
3.1.1 Importance sampling
If importance sampling methods are used, samples are generated from instead of . As a result, (1) is expanded as
and so should be a fitting algorithm for instead of just . As a result a few modifications are needed including the fact that,
and that in the calculation of and (10), is used instead of just .
3.1.2 Unable to integrate over
If the samples are generated from , but the choice for cannot be integrated over , there are two options. The first option is to use another function as a fitting algorithm for over which is integrable. We expand (1) as
Similarly to the above case, when calculating and (10), is replaced with .
Alternately, when is significantly less computationally expensive than , itself can be estimated by Monte Carlo techniques. When estimating from additional samples, it should be the case that so that errors in the estimate of do not cause significant errors in the estimate of .
3.2 A Simple Example
We attempt to find the expected value of , where is generated according to a uniform distribution between and . Thus, the integral of concern is
The actual value of can be calculated to be 0.7069. Twenty samples are generated and divided into five folds, with each fold’s testing set containing four points and each corresponding training set containing the other sixteen points. While the function of interest is a sixth order polynomial, was chosen to be a third-order polynomial: . is found by choosing the ’s which minimize the squared error of the data in (using the standard pseudo-inverse). - are set similarly. Next, is evaluated for each test point in .
The following parameters are calculated to find :
Finally, equation (10) is applied to calculate .
Using the Monte Carlo estimate alone gives , and using a fit to all of the data points alone gives . By using StackMC, we find , much closer to the true value than either of the individual estimates. Details of the calculations can be seen in Tables 1 and 2, and the fit from the first fold can be seen in Fig. 2.
We test the performance of StackMC on a number of different problems. In the first set of example cases the problems have known analytic results and an exact analysis of the performance of StackMC is examined. In the second set, StackMC is applied to two aerospace problems in the literature. While an exact solution to the aerospace problems is not available, an approximate answer is obtained from a Monte Carlo estimate using 100,000 samples. For each example problem, a range of number of samples was tested using two thousand simulations at each value of . In each case, 10-fold cross-validation was used. Unless otherwise noted, the plots in this section have three lines which show the mean squared error versus the number of sample points. The lines represent the average squared error from using just the prediction of Monte Carlo (green), the fit to all the samples (red), and StackMC (blue).
4.1 Analytic Test Cases
4.1.1 Ten-dimensional Rosenbrock function under a uniform , polynomial fitter.
The D-dimensional Rosenbrock function is given by
and is drawn from a uniform distribution over the [-3,3] hypercube.
The fitting algorithm is chosen as a third-order polynomial in each dimension whose form is
where the are free parameters that are set from the data samples.
A comparison of the error of StackMC can be seen in Fig. 3 for the ten-dimensional version of the Rosenbrock function. At low numbers of sample points, Monte Carlo is more accurate, on average, than the polynomial fit to all the data points, but at higher numbers of points the polynomial outperforms Monte Carlo. Throughout the entire range of number of samples, StackMC performs at least as well as the best of the two. Additionally, as can be seen in Fig. 4, the polynomial fitter is actually a biased fitter for this example problem. StackMC is able to use cross-validation to remove the bias of the fitter while keeping the accuracy of its estimate.
4.1.2 Ten dimensional Rosenbrock function under a Gaussian .
This is the same set-up as above, except that is generated according to a multivariate Gaussian distribution, each dimension independent with and .
Much like the case above, in Fig. 5 it can be seen that at low numbers of sample points Monte Carlo outperforms the polynomial fit, whereas at high sample points the polynomial does better than Monte Carlo alone. However, for all numbers of sample points, StackMC outperforms the other two algorithms.
4.1.3 20 Dimensional BTButterfly, uniform , Fourier fitter
In the previous examples, the Rosenbrock function is a 4 order polynomial, and the fitting algorithm is a 3 order polynomial, so StackMC had reasonable fits to use to improve upon Monte Carlo. In order to challenge StackMC, a function we call the BTButterfly was created so that it would be very difficult to fit accurately.
Like the Rosenbrock function, its general form is given by
with the contour plot of shown in Fig. 6. is generated uniformly over the [-3,3] hypercube.
A Fourier expansion for was chosen whose form is given by
Results for this function are shown in Fig.7, and exhibit the same trends seen previously; StackMC has as low of an error as the lowest of the two methods individually.
4.2 Aerospace Applications
4.2.1 Future Aircraft Uncertainty Quantification (UQ)
In Ref., the authors use the Program for Aircraft Synthesis Studies  (PASS), a conceptual aircraft design tool, to predict the fuel burn of future aircraft given certain assumptions about technology advancement in the 2020 and 2030 time frames. In their predictions for single-aisle aircraft in 2020, the authors model eight probabilistic variables representing different effects of improvements in aircraft technology (propulsion, structures, aerodynamics). Each of the variables is represented by a unique beta distribution. The authors generated 15,000 samples (each representing one optimized aircraft) from which they measured the expected fuel-burn metric and the standard deviation of the expected fuel-burn metric.
To apply StackMC, we choose a third-order polynomial fitter. A polynomial fitter is convenient for this case because the over a beta distribution is easily found analytically as follows:
where and are the two parameters of the beta distribution, (not to be confused with the StackMC parameter ).
A set of 100,000 samples were generated and the mean and standard deviation of their function values were taken as the “true” expected value and standard deviation. The standard deviation is defined as , and therefore, to find the standard deviation we can run StackMC twice, once to find and a second time to find . The results of applying StackMC to find the expected value and standard deviation are shown in Fig. 8 and Fig. 9, respectively.
The exact same trend is seen here as in all of the analytic problems. At low numbers of samples, where the fit to all the data points performs badly, StackMC does as well as Monte Carlo. At high numbers of samples, where the fit greatly outperforms Monte Carlo, StackMC does as well as the fitter. In between these two extremes, StackMC outperforms both algorithms.
At 1500 samples, StackMC reduces the error by roughly an order of magnitude compared with the Monte Carlo result used by the authors. Each sample takes about 6 seconds to generate, whereas StackMC takes less than a half a second to form the 10 fits, calculate , and calculate (10) . For a computational cost of less than one additional sample, StackMC has reduced the number of samples needed by 90% for the same expected error.
4.2.2 Sonic Boom Uncertainty Quantification
The uncertainty quantification of sonic boom noise is used as the second application case. Unlike the aircraft design test case, the response surface for the sonic boom noise signature is not smooth; the output can vary significantly with slight adjustments to the input parameters. Colonno and Alonso recently created a new sonic boom propagation tool, SUBoom, and used it to analyze the robustness of several aircraft pressure signatures optimized for minimal boom noise. Specifically, in their high-fidelity near-field case, they have 62 uncertain variables: 4 representing aircraft parameters such as cruise Mach number and roll angle, and 58 representing uncertainties in the near-field pressure signal. Like in the Aircraft UQ case, 100,000 samples were generated and used as the true mean. A third-order polynomial was used as . The results can be seen in Fig. 10.
Even with the discontinuities in the function space, the same performance for StackMC is seen. StackMC does as well as the best of Monte Carlo and the fitting algorithm. The reduction in error is not as great here as in the Aircraft UQ case because a polynomial is not a good fit to , though using domain knowledge to improve the accuracy of the fitting algorithm would further increase the performance of StackMC. Despite the relatively poor performance of the fitting algorithm, it is still better to use StackMC than to not regardless of the number of sample points used. In some senses, this test case demonstrates one of the strengths of StackMC: at virtually no additional computation cost the user is guaranteed the highest accuracy without having to worry about the appropriateness of the methodology for a given number of samples.
In this paper, we have introduced a new technique, Stacked Monte Carlo, which reduces error of Monte Carlo sampling by blending several different fitters of . StackMC is unbiased, (thus retaining a major advantage of Monte Carlo), and is able to use cross validation to effectively incorporate the use of a fitting function.
StackMC is an extremely generic post-processing technique. It is applied after the samples are generated and makes no assumptions about the generation method, and therefore StackMC can be used not only with simple MC, but also with other sample generation techniques such as Importance Sampling, Quasi-Monte Carlo, and Markov-Chain Monte Carlo. Furthermore, StackMC makes no assumptions about ; it not only applies to smooth functions, but also to discontinuous functions, and even functions with discrete variables. In a future paper,  we will present results of the application of StackMC to different input regimes and different sample generation methods. We will examine higher dimensional spaces, explore the application of StackMC to multi-fidelity methods, and extended StackMC to incorporate multiple fitting functions.
The computation time of StackMC is dominated by forming the fits . As a result, StackMC is only affected by the curse of dimensionality insofar as the fitting algorithm is. In both of the aerospace applications, the additional cost of the entire StackMC algorithm was virtually non-existant; there are significant improvements in accuracy for less computation time than the cost of one additional function evaluation. The only assumptions made about are that it be able to predict the value at new sample locations and we are able to evaluate analytically (and even this second assumption is relaxed). As shown in the BTButterfly and sonic boom example cases, the fit does not need to be particularly good to see improvement, and so the choice for can be modified as computational effort and domain knowledge allows. StackMC does not eliminate the need for finding better fitters; a better fitter will always lead to an improved result. With the exception of the error in the mean test (whose value we did not change for any of the example cases), StackMC requires no user-set parameters that need to be heuristically tuned to see good results.
The version of StackMC implemented in this paper is relatively naive. Instead of taking the mean of the results from the different fits, taking a weighted combination of the fits could improve the estimate. StackMC currently only partitions the data into folds once, but we could imagine re-partitioning the same samples several times to have more fitters. We use Monte Carlo to estimate the integral term in (7), but using a fitter, or even using StackMC recursively, could improve the estimates of . Furthermore, was set as a constant, but in general could vary between the folds or could even be a function of so the confidence in the fit varies spatially.
Despite this simplistic implementation, StackMC performs at least as well as the better of MC and the fitting function, and for a range of sample points outperforms both. Normally, there is a transition number of samples at which using a fit to all the data samples has a smaller average error than MC, and it is hard to know a priori if it is better to use the fit. StackMC eliminates the need to decide whether or not to use a fit; it will never be harmful to do so. While it is not true that for any set of data samples is closer to than , on average StackMC reduces the expected error. StackMC is a very generic method for the post-processing of data samples; it can be used by anyone trying to estimate an integral or expected value of a function where is known.
-  Haugh, M., “Variance Reduction Methods II,” Monte Carlo Simulation: IEOR E4703, 2004.
-  Kocis, L. and Whiten, W. J., “Computational Investigations of Low-Discrepancy Sequences,” ACM Transactions on Mathematical Software, Vol. 23, No. 2, 1997, pp. 266–294.
-  Sobol, I., “Uniformly distributed sequences with additional uniformity properties,” USSR Comput. Math. Math. Phys, Vol. 16, No. 5, 1976, pp. 236–242.
-  Wolpert, D., “Stacked Generalizations,” Neural Networks, Vol. 5, No. 2, 1992, pp. 241–260.
-  Breiman, L., “Stacked regressions,” Machine Learning, Vol. 24, pp. 49–64.
-  Smyth, P. and Wolpert, D., “Stacked density estimation,” Proceedings of the 1997 conference on Advances in neural information processing systems 10, NIPS ’97, MIT Press, Cambridge, MA, USA, 1998, pp. 668–674.
-  Kim, K. and Bartlett, E. B., “Error Estimation by Series Association for Neural Network Systems,” Neural Computation, Vol. 7, No. 4, 1995, pp. 799–808.
-  Rajnarayan, D. and Wolpert, D., “Exploiting Parametric Learning to Improve Black-Box Optimization,” .
-  Sill, J., Takács, G., Mackey, L., and Lin, D., “Feature-Weighted Linear Stacking,” CoRR, Vol. abs/0911.0460, 2009.
-  Clarke, B., “Comparing Bayes Modeling Averaging and Stacking when Model Aprroximation Error Cannot Be Ignored,” Journal of Machine Learning Research, Vol. 4:683-712, 2003.
-  Haugh, M., “Variance Reduction Methods I,” Monte Carlo Simulation: IEOR E4703, 2004.
-  Duda, R. O., Hart, P. E., and Stork, D. G., Pattern Classification, Wiley, New York, 2nd ed., 2001.
-  Economon, T., Copeland, S., Alonso, J., Zeinali, M., and Rutherford, D., “Design and Optimization of Future Aircraft for Assessing the Fuel Burn Trends of Commercial Aviation,” 49th AIAA Aerospace Sciences Meeting, Orlando, Fl., January 2011, AIAA Paper 2011-267.
-  Kroo, I., “An Interactive System for Aircraft Design and Optimization,” Aerospace Design Conference, Irvine, Ca., 1992, AIAA 1992-1190.
-  Colonno, M. and Alonso, J., “Sonic Boom Minimization Revisited: The Robustness of Optimal Low-Boom Designs,” 13th AIAA/ISSMO Multidisciplinary Analysis Optimization Conference, Fort Worth, Texas, September 2010, AIAA-2010-9364.
-  Tracey, B., Wolpert, D., and Alonso, J., “Variance Reduction of Monte Carlo Methods Using Stacking,” In Submission.