18 Uncertainty Propagation
As discussed in the previous chapter, probabilistic approaches to optimization under uncertainty model some of the inputs to the objective function as a probability
distribution. This chapter discusses how to propagate known input distributions
to estimate quantities associated with the output distribution, such as the mean
and variance of the objective function. There are a variety of approaches to uncertainty propagation, some based on mathematical concepts such as Monte Carlo, theTaylor series approximation, orthogonal polynomials, and Gaussian processes.
These approaches differ in the assumptions they make and the quality of theirestimates.

18.1 Sampling Methods
The mean and variance of the objective function at a particular design point canbe approximated using Monte Carlo integration,1 which approximates the integral
using m samples z(1) , . . . , z(m) , from the distribution p over Z . These estimates
are also called the sample mean and sample variance:Ez∼ p [ f (z)] ≈ µ̂ =Varz∼ p [ f (z)] ≈ ν̂ =1 m
f (z( i ) )
m i∑
=1
1 m
f (z( i ) )2
m i∑
=1

(18.1)!

− µ̂2(18.2)In the equation above, and for the rest of this chapter, we drop x from f (x, z) for
notational convenience, but the dependency on x still exists. For each new design
point x in our optimization process, we recompute the mean and variance.

Alternatively, quasi Monte Carlointegration can be used to produce
estimates with faster convergence
as discussed in chapter 13.

1

322 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

A desirable property of this sampling-based approach is that p does not need
to be known exactly. We can obtain samples directly from simulation or real-world
experiments. A potential limitation of this approach is that many samples may
be required before there is convergence to a suitable estimate. The variance ofthe sample mean for a normally distributed f is Var[µ̂] = ν/m, where ν is thetrue variance of f . Thus, doubling the number of samples m tends to halve thevariance of the sample mean.

18.2 Taylor Approximation
Another way to estimate µ̂ and ν̂ is to use the Taylor series approximation for f at
a fixed design point x.2 For the moment, we will assume that the n componentsof z are uncorrelated and have finite variance. We will denote the mean of the
distribution over z as µ and the variances of the individual components of z as
ν.3 The following is the second-order Taylor series approximation of f (z) at the
point z = µ:n
1 n n ∂2 f
∂f
( zi − µi ) + ∑ ∑(z − µi )(z j − µ j )fˆ(z) = f (µ) + ∑∂zi
2 i=1 j=1 ∂zi ∂z j ii =1(18.3)For a derivation of the mean andvariance of a general function of n
random variables, see H. Benaroya
and S. M. Han, Probability Models
in Engineering and Science. Taylor &Francis, 2005.
3
If the components of z are uncorrelated, then the covariance matrix
is diagonal and ν would be the vector composed of the diagonal elements.
2

From this approximation, we can analytically compute estimates of the mean
and variance of f :1 n ∂2 f
ν
2 i
2 i∑
z=µ
=1 ∂zi
!2
2
n n
n 
1
∂2 f
∂f
νi + ∑ ∑νi νjν̂ = ∑∂zi
2 i=1 j=1 ∂zi ∂z j
z=µ
i =1

µ̂ = f (µ) +(18.4)(18.5)The higher-order terms can be neglected to obtain a first-order approximation:µ̂ = f (µ)

n ν̂ = ∑

i =1

∂f
∂zi

2

(18.6)

νiz=µWe can relax the assumption that the components of z are uncorrelated, but it
makes the mathematics more complex. In practice, it can be easier to transformthe random variables so that they are uncorrelated. We can transform a vector ofn correlated random variables c with covariance matrix C into m uncorrelated© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

18 . 3. p ol y n omi a l c ha os 323

random variables z by multiplying by an orthogonal m × n matrix T containing
the eigenvectors corresponding to the m largest eigenvalues of C. We have z = Tc.4
The Taylor approximation method is implemented in algorithm 18.1. First- andsecond-order approximations are compared in example 18.1.
using ForwardDiff
function taylor_approx(f, μ, ν, secondorder=false)
μhat = f(μ)∇ = (z -> ForwardDiff.gradient(f, z))(μ)νhat = ∇.^2⋅ν
if secondorderH = (z -> ForwardDiff.hessian(f, z))(μ)μhat += (diag(H)⋅ν)/2
νhat += ν⋅(H.^2*ν)/2
end
return (μhat, νhat)
end

18.3

It is also common to scale the outputs such that the covariance matrix becomes the identity matrix.
This process is known as whitening.
J. H. Friedman, “Exploratory Projection Pursuit,” Journal of the American Statistical Association, vol. 82,
no. 397, pp. 249–266, 1987.
4

Algorithm 18.1. A method forautomatically computing the Taylor approximation of the mean
and variance of objective function f at design point x withnoise mean vector μ and variancevector ν. The Boolean parametersecondorder controls whether the
first- or second-order approximation is computed.

Polynomial Chaos

Polynomial chaos is a method for fitting a polynomial to evaluations of f (z) andusing the resulting surrogate model to estimate the mean and variance. We will
begin this section by discussing how polynomial chaos is used in the univariatecase. We will then generalize the concept to multivariate functions and show
how to obtain estimates of the mean and variance by integrating the function
represented by our surrogate model.

18.3.1 Univariate
In one dimension, we approximate f (z) with a surrogate model consisting of kpolynomial basis functions, b1 , . . . , bk :kf (z) ≈ fˆ(z) = ∑ θi bi (z)(18.7)

i =1© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

324 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

Consider the objective function f ( x, z) = sin( x + z1 ) cos( x + z2 ), where z1
and z2 are zero-mean Gaussian noise with variances 0.1 and 0.2, respectively.
The first and second partial derivatives of f with respect to the zs are
∂f
=
∂z1

∂2 f
= − sin( x + z1 ) cos( x + z2 )∂z22cos( x + z1 ) cos( x + z2 )∂f
= − sin( x + z1 ) sin( x + z2 )∂z2∂2 f
= − sin( x + z1 ) cos( x + z2 )∂z21

Example 18.1. The Taylor approximation applied to a univariate design problem with twodimensional Gaussian noise.

∂2 f= − cos( x + z1 ) sin( x + z2 )∂z1 ∂z2

which allow us to construct the Taylor approximation:µ̂( x ) = 0.85 sin( x ) cos( x )ν̂( x ) = 0.1 cos4 ( x ) + 0.2 sin4 ( x ) + 0.045 sin2 ( x ) cos2 ( x )We can use taylor_approx for a given x using:taylor_approx(z->sin(x+z[1])*cos(x+z[2]), [0,0], [0.1,0.2])

mean

0.5
0

variance

−0.5
0.2
0.15
0.1
0.05
0

−3

true

−2

−1

0
x

first-order Taylor approx.

1

2

3

second-order Taylor approx.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

18 . 3. p ol y n omi a l c ha os 325

In contrast with the Monte Carlo methods discussed in section 18.1, our samplesof z do not have to be randomly drawn from p. In fact, it may be desirable toobtain samples using one of the sampling plans discussed in chapter 13. We will
discuss how to obtain the basis coefficients in section 18.3.2.
We can use the surrogate model fˆ to estimate the mean:h i(18.8)µ̂ = E fˆ

=

=

Z

Z

(18.9)

fˆ(z) p(z) dz

Zk(18.10)∑ θi bi (z) p(z) dz
Z
i =1

k

Z

= ∑ θi
i =1

= θ1

Z

(18.11)bi (z) p(z) dz

Z

b1 (z) p(z) dz + . . . + θkZ

We can also estimate the variance:h i2 ν̂ = E fˆ − E fˆZZ

bk (z) p(z) dz(18.13)h i2 
h i= E fˆ2 − 2 fˆ E fˆ + E fˆ

(18.14)h i
h i2= E fˆ2 − E fˆ

=

Z

=

Z

=

Z

Z

(18.15)(18.16)fˆ(z)2 p(z) dz − µ2kk(18.17)∑ ∑ θi θ j bi (z)bj (z) p(z) dz − µ2

Z i =1 j =1
k

∑

Z

i =1

k

Z

= ∑ θi2i =1

θi2 bi (z)2 + 2

k i −1!

(18.18)bi (z)b j (z) p(z) dz − µ2(18.19)∑ ∑ θi θ j bi (z)bj (z) p(z) dz − µ2

i =2 j =1

k i −1Z(18.12)bi (z)2 p(z) dz + 2 ∑ ∑ θi θ j
i =2 j =1

Z

Z

The mean and variance can be efficiently computed if the basis functions arechosen to be orthogonal under p. Two basis functions bi and b j are orthogonal withrespect to a probability density p(z) ifZZbi (z)b j (z) p(z) dz = 0

if i 6= j(18.20)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

326 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

If the chosen basis functions are all orthogonal to one another and the firstbasis function is b1 (z) = 1, the mean is:µ̂ = θ1

= θ1
= θ1
= θ1

Z

ZZ
Z

Z
Z

b1 (z) p(z) dz + θ2Z

b1 (z)2 p(z) dz + θ2ZZ

b2 (z) p(z) dz + · · · + θk

Z

Z

Z

bk (z) p(z) dz

b1 (z)b2 (z) p(z) dz + · · · + θkZZ

b1 (z)bk (z) p(z) dz(18.21)(18.22)(18.23)p(z) dz + 0 + · · · + 0(18.24)

Similarly, the variance is:kν̂ = ∑ θi2
i =1k= ∑ θi2
i =1

= θ12

Z

k

Z

Z

Z

bi (z)2 p(z) dz + 2 ∑ ∑ θi θ j
i =2 j =1

Z

b1 (z)2 p(z) dz + ∑ θi2
i =2

Z

Z

Z

Z

bi (z)b j (z) p(z) dz − µ2bi (z)2 p(z) dz(18.25)(18.26)bi (z)2 p(z) dz − µ2k= ∑ θi2i =2

k i −1ZZ

Z

bi (z)2 p(z) dz − θ12(18.27)(18.28)The mean thus falls immediately from fitting a surrogate model to the observed data, and the variance can be very efficiently computed given the valuesR2
5
Z bi ( z ) p ( z ) dz for a choice of basis functions and probability distribution. Example 18.2 uses these procedures to estimate the mean and variance with different
sample sizes.
Polynomial chaos approximates the function using kth degree orthogonal
polynomial basis functions with i ∈ {1, . . . , k + 1} and b1 = 1. All orthogonal
polynomials satisfy the recurrence relation:(z − α )b (z)for i = 1i i(18.29)bi + 1 ( z ) =(z − αi )bi (z) − β i bi−1 (z) for i > 1

5Integrals of this form can be efficiently computed using Gaussian quadrature, covered in appendix C.8.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

18 . 3. p ol y n omi a l c ha os 327

Example 18.2. Estimating the expected value of an unknown objective function using polynomialchaos.

Consider optimizing the (unknown) objective function2f ( x, z) = 1 − e−( x+z−1) − 2e−( x+z−3)

2

E[ f | x]

E[ f | x]

E[ f | x]with z known to be drawn from a zero-mean unit-Gaussian distribution.
The objective function, its true expected value, and estimated expectedvalues with different sample counts are plotted below. The estimated expectedvalue is computed using third-order Hermite polynomials.
noise-free
exact
95% interval
mean

1
0

−1

10 samples

1
0

−1

30 samples

1
0

−1

−2

50 samples

−1

0

1

2
x

3

4

5

6

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

328 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

with b1 (z) = 1 and weightsRz bi (z)2 p(z) dz
αi = RZbi (z)2 p(z) dz
RZbi (z)2 p(z) dz
βi = R Z
2
Z bi −1 ( z ) p ( z ) dz

(18.30)

The recurrence relation can be used to generate the basis functions. Each basis
function bi is a polynomial of degree i − 1. The basis functions for several common probability distributions are given in table 18.1, can be generated using themethods in algorithm 18.2, and are plotted in figure 18.1. Example 18.3 illustrates
the effect the polynomial order has on the estimates of the mean and variance.
Distribution

DomainDensity

Name

Recursive Formhk

Uniform[−1, 1]

1
2

LegendredLek ( x ) = 2k1k! dx
kExponential[0, ∞)e− x

Laguerreddx Lak ( x ) =Unit Gaussian(−∞, ∞)

2√1 e− x /2
2π

Hermitex2 − 1ddx − 1



 i
k

Lak−1

d
Hk −1
Hk ( x ) = xHk−1 − dx

Closed Formi −1 
1−x j−i − 2i−1bi ( x ) = ∑
j
j
2
j =0
i −1 

bi ( x ) = ∑j =0bi ( x ) =
i − 1 (−1) j jxj!
j⌊(i −1)/2⌋

Basis functions for arbitrary probability density functions and domains canbe constructed both analytically and numerically.6 The Stieltjes algorithm7 (algorithm 18.3) generates orthogonal polynomials using the recurrence relationin equation (18.29). Example 18.4 shows how the polynomial order affects theestimates of the mean and variance.

18.3.2 Coefficients
The coefficients θ1 , . . . , θk in equation (18.7) can be inferred in two different ways.
The first way is to fit the values of the samples from Z using the linear regression
method discussed in section 14.3. The second way is to exploit the orthogonalityof the basis functions, producing an integration term amenable to Gaussian
quadrature.

∑

j =0

( i − 1) !

i −1(−1) 2 − j(2x )2j(2j)!( i−2 1 − j)!

Table 18.1. Orthogonal polynomial
basis functions for several common
probability distributions.
The polynomials can be scaled bya nonzero factor. It is convention to
set b1 ( x ) = 1.

6

T. J. Stieltjes, “Quelques Recherches sur la Théorie des Quadratures Dites Mécaniques,” Annales Scientifiques de l’École Normale
Supérieure, vol. 1, pp. 409–426, 1884.
in French. An overview in English
is provided by W. Gautschi, Orthogonal Polynomials: Computation andApproximation. Oxford University
Press, 2004.
7© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

18 . 3. p ol y n omi a l c ha os 329

Algorithm 18.2. Methods for constructing polynomial orthogonalbasis functions, where i indicatesthe construction of bi .

using Polynomials
function legendre(i)n = i-1
p = Polynomial([-1,0,1])^n
for i in 1 : n
p = derivative(p)endreturn p / (2^n * factorial(n))
endfunction laguerre(i)p = Polynomial([1])
for j in 2 : i
p = integrate(derivative(p) - p) + 1
end
return pendfunction hermite(i)p = Polynomial([1])
x = Polynomial([0,1])
for j in 2 : ip = x*p - derivative(p)endreturn p
end

Laguerre

Legendre

Hermite

20

b1
b2
b3
b4
b5
b6

1
10
10
0

0

0

−10
−1
−1

−0.5

0
z

0.5

1

−20

−10
0

5

10
z

15

−4

−2

0
z

2

4

Figure 18.1. Orthogonal basis functions for uniform, exponential, andunit Gaussian distributions.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

330 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

Consider the function f (z) = sin(πz) with input z drawn from a uniform
distribution over the domain [−1, 1]. The true mean and variance can becomputed analytically:
Z b

Z 1

1
sin(πz) dz = 0
2
a
−1
Z b
Z 1
1
1
ν=
f (z)2 p(z) dz − µ2 =sin2 (πz) dz − 0 =
2
2
a
−1

µ=

f (z) p(z) dz =Example 18.3. Legendre polynomials used to estimate the mean andvariance of a function with a uniformly distributed input.

(18.31)(18.32)

Suppose we have five samples of f at z = {−1, −0.2, 0.3, 0.7, 0.9}. We
can fit a Legendre polynomial to the data to obtain our surrogate model fˆ.
Polynomials of different degrees yield:
true
i=1
i=2
i=3
i=4
i=5

1

0

µ = +0.268µ = +0.207µ = +0.170µ = −0.056µ = −0.136

ν = 0.000
ν = 0.063
ν = 0.075
ν = 0.573
ν = 0.757

−1

−1

−0.5

0
z

0.5

1

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

18 . 3. p ol y n omi a l c ha os 331

Algorithm 18.3. The Stieltjes algorithm for constructing the next
polynomial basis function bi+1 according to the orthogonal recurrence relation, where bs contains{b1 , . . . , bi }, p is the probability distribution, and dom is a tuple containing a lower and upper boundfor z. The optional parameter ϵ controls the absolute tolerance of thenumerical integration. We makeuse of the Polynomials.jl package.

using Polynomials
function orthogonal_recurrence(bs, p, dom, ϵ=1e-6)
i = length(bs)c1 = quadgk(z->z*bs[i](z)^2*p(z), dom..., atol=ϵ)[1]c2 = quadgk(z-> bs[i](z)^2*p(z), dom..., atol=ϵ)[1]α = c1 / c2
if i > 1
c3 = quadgk(z->bs[i-1](z)^2*p(z), dom..., atol=ϵ)[1]β = c2 / c3
return Polynomial([-α, 1])*bs[i] - β*bs[i-1]else
return Polynomial([-α, 1])*bs[i]endend

We multiply each side of equation (18.7) by the jth basis and our probability
density function and integrate:k(18.33)f ( z ) ≈ ∑ θ i bi ( z )
i =1

Z

Z

f (z)b j (z) p(z) dz ≈

Z

k

i =1

k

Z

= ∑ θii =1

Z(18.34)bi (z)b j (z) p(z) dz(18.35)∑ θi bi (z) bj (z) p(z) dz

Z

= θj

!

Z

Z

b j (z)2 p(z) dz(18.36)

where we made use of the orthogonality property from equation (18.20).
It follows that the jth coefficient is:R
f (z)b j (z) p(z) dz
θ j = ZR
(18.37)2
Z b j ( z ) p ( z ) dz
The denominator of equation (18.37) typically has a known analytic solution
or can be inexpensively precomputed. Calculating the coefficient thus primarily
requires solving the integral in the numerator, which can be done numericallyusing Gaussian quadrature.8

Gaussian quadrature is implemented in QuadGK.jl via thequadgk function, and is covered inappendix C.8. Quadrature rules
can also be obtained using the
eigenvalues and eigenvectors of atri-diagonal matrix formed using
the coefficients αi and β i fromequation (18.30). G. H. Golub andJ. H. Welsch, “Calculation of Gauss
Quadrature Rules,” Mathematics
of Computation, vol. 23, no. 106,
pp. 221–230, 1969.
8© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

332 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

Consider the function f (z) = sin(πz) with input z drawn from a truncatedGaussian distribution with mean 3 and variance 1 over the domain [2, 5]. The
true mean and variance are:
µ=
ν=

Z b
a

Z b
a

f (z) p(z) dz =Z 5
2

Example 18.4. Legendre polynomials constructed using the Stieltjes
method to estimate the mean andvariance of a function with a random variable input.

sin(πz) p(z) dz ≈ 0.104

f (z)2 p(z) dz − µ2 =Z 52

sin2 (πz) p(z) dz − 0.1042 ≈ 0.495

where the probability density of the truncated Gaussian is: R 5 N (z|3,1)
if z ∈ [2, 5]
2 N ( τ |3,1) dτ
p(z) =
0otherwise

Suppose we have five samples of f at z = {2.1, 2.5, 3.3, 3.9, 4.7}. We canfit orthogonal polynomials to the data to obtain our surrogate model fˆ.
Polynomials of different degrees yield:
true
i=1
i=2
i=3
i=4
i=5

2

1

µ = +0.200µ = +0.204µ = −0.036µ = +0.037µ = +0.055

ν = 0.000
ν = 0.001
ν = 0.159
ν = 0.273
ν = 0.478

0

−1
2

3

4

5

z

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

18 . 3. p ol y n omi a l c ha os 333

18.3.3 Multivariate
Polynomial chaos can be applied to functions with multiple random inputs. Multivariate basis functions over m variables are constructed as a product over univariate orthogonal polynomials:mbi ( z ) = ∏ b a j ( z j )(18.38)j =1

where a is an assignment vector that assigns the a j th basis function to the jthrandom component. This basis function construction is demonstrated in example 18.5.
Consider a three-dimensional polynomial chaos model of which one of themultidimensional basis functions is b(z) = b3 (z1 )b1 (z2 )b3 (z3 ). The corresponding assignment vector is a = [3, 1, 3].
A common method for constructing multivariate basis functions is to generateunivariate orthogonal polynomials for each random variable and then to constructa multivariate basis function for every possible combination.9 This procedure
is implemented in algorithm 18.4. Constructing basis functions in this mannerassumes that the variables are independent. Interdependence can be resolved
using the same transformation discussed in section 18.2.

Example 18.5. Constructing a multivariate polynomial chaos basis
function using equation (18.38).

Here the number of multivariate
exponential basis functions grows
exponentially in the number ofvariables.
9

Algorithm 18.4. A method for constructing multivariate basis functions where bases1d contains listsof univariate orthogonal basis functions for each random variable.

function polynomial_chaos_bases(bases1d)bases = []for a in Iterators.product(bases1d...)push!(bases,z -> prod(b(z[i]) for (i,b) in enumerate(a)))
endreturn bases
end

A multivariate polynomial chaos approximation with k basis functions is stilla linear combination of termskf (z) ≈ fˆ(z) = ∑ θi bi (z)(18.39)i =1© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

334 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

where the mean and variance can be computed using the equations in section 18.3.1, provided that b1 (z) = 1.

18.4

Bayesian Monte Carlo

Gaussian processes, covered in chapter 16, are probability distributions overfunctions. They can be used as surrogates for stochastic objective functions. We canincorporate prior information, such as the expected smoothness of the objectivefunction, in a process known as Bayesian Monte Carlo or Bayes-Hermite Quadrature.
Consider a Gaussian process fit to several points with the same value forthe design point x but different values for the uncertain point z. The Gaussian
process obtained is a distribution over functions based on the observed data.
When obtaining the expected value through integration, we must consider the
expected value of the functions in the probability distribution represented by theGaussian process p( fˆ):hh ii
Ez∼ p [ f ] ≈ E fˆ∼ p( fˆ) Ez∼ p fˆ(18.40)Z Z=
fˆ(z) p(z) dz p( fˆ) d fˆ(18.41)
F̂
Z


Z
Z
=
fˆ(z) p( fˆ) d fˆ p(z) dz(18.42)

=

Z

Z

Z

F̂

µ̂(z) p(z) dz(18.43)where µ̂(z) is the predicted mean under the Gaussian process and F̂ is the spaceof functions. The variance of the estimate ishh ii
Varz∼ p [ f ] ≈ Var fˆ∼ p( fˆ) Ez∼ p fˆ
2
Z
Z Z
fˆ(z) p(z) dz −E fˆ∼ p( fˆ) [ fˆ(z′ )] p(z′ ) dz′ p( fˆ) d fˆ=ZZF̂Z Z Z hiihfˆ(z) − E fˆ∼ p( fˆ) [ fˆ(z)] fˆ(z′ ) − E fˆ∼ p( fˆ) [ fˆ(z′ )] p( fˆ) d fˆ p(z) p(z′ ) dz dz′==

ZZ ZZ F̂Z

Z

Cov( fˆ(z), fˆ(z′ )) p(z) p(z′ ) dz dz′

where Cov is the posterior covariance under the Gaussian process:Cov( fˆ(z), fˆ(z′ )) = k(z, z′ ) − k (z, Z )K( Z, Z )−1 k ( Z, z′ )(18.48)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com(18.44)(18.45)(18.46)(18.47)

18 . 4 . ba y e si a n mon te c a rl o 335

where Z contains the observed inputs.
Analytic expressions exist for the mean and variance for the special case wherez is Gaussian.10 Under a Gaussian kernel,
2 !
1 n xi − xi′′(18.49)k(x, x ) = exp − ∑2 i =1wi2
the mean for Gaussian uncertainty z ∼ N (µz , Σz ) ishh ii
Ez∼ p [ f ] ≈ E fˆ∼ p( fˆ) Ez∼ p fˆ = q⊤ K−1 ynk(x, x′ ) = ∏ r ( xi − xi′ )i =1(18.50)with

qi = |W−1

Σz + I|

−1/2


⊤

1
−1
(i )
(i )
µz − z(18.51)
exp − µz − z( Σz + W)2

where W = diag[w12 , . . . , w2n ], and we have constructed our Gaussian processusing samples (z(i) , yi ) for i ∈ {1, . . . , m}.11
The variance isVarz∼ p [ f ] = |2W−1 Σz + I|−1/2 − q⊤ K−1 q(18.52)

Even when the analytic expressions are not available, there are many problems
for which numerically evaluating the expectation is sufficiently inexpensive that
the Gaussian process approach is better than a Monte Carlo estimation.
Bayesian Monte Carlo is implemented in algorithm 18.5 and is worked out inexample 18.6.
function bayesian_monte_carlo(GP, w, μz, Σz)
W = Matrix(Diagonal(w.^2))
invK = inv(K(GP.X, GP.X, GP.k))q = [exp(-((z-μz)⋅(inv(W+Σz)*(z-μz)))/2) for z in GP.X]
q .*= (det(W\Σz + I))^(-0.5)
μ = q'*invK*GP.yν = (det(2W\Σz + I))^(-0.5) - (q'*invK*q)[1]
return (μ, ν)
endIt is also required that the covariance function obey the product correlation rule, that is, it can be written as the product of a univariate
positive-definite function r:

10

Analytic results exist for polynomial kernels and mixtures ofGaussians. C. E. Rasmussen and
Z. Ghahramani, “Bayesian Monte
Carlo,” in Advances in Neural Information Processing Systems (NIPS),
2003.

See A. Girard, C. E. Rasmussen, J. Q. Candela, and R.
Murray-Smith, “Gaussian Process
Priors with Uncertain Inputs—
Application to Multiple-Step
Ahead Time Series Forecasting,”
in Advances in Neural Information
Processing Systems (NIPS), 2003.
11

Algorithm 18.5. A method for obtaining the Bayesian Monte Carloestimate for the expected value ofa function under a Gaussian process GP with a Gaussian kernel withweights w, where the variables are
drawn from a normal distribution
with mean μz and covariance Σz.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

336 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

Consider again estimating the expected value and variance of f ( x, z) =sin( x + z1 ) cos( x + z2 ), where z1 and z2 are zero-mean Gaussian noise withvariances 1 and 1/2, respectively: µz = [0, 0] and Σz = diag([1, 1/2]).
We use Bayesian Monte Carlo with a Gaussian kernel with unit weights
w = [1, 1] for x = 0 with samples Z = {[0, 0], [1, 0], [−1, 0], [0, 1], [0, −1]}.
We compute:
#
"
1 0
W=
0 1


1
0.607 0.607 0.607 0.607
0.60710.135 0.368 0.368K = 0.607 0.135
10.368 0.368

0.607 0.368 0.368
10.135
0.607 0.368 0.368 0.135

Example 18.6. An example of using
Bayesian Monte Carlo to estimate
the expected value and variance ofa function.
The figure compares theBayesian Monte Carlo method tothe sample mean for estimating the
expected value of a function. The
same randomly sampled z values
generated for each evaluated x
were input to each method.

1q = [0.577, 0.450, 0.450, 0.417, 0.417]

Ez∼ p [ f ] = 0.0

Varz∼ p [ f ] = 0.327
Below we plot the expected value as a function of x using the same approach with ten random samples of z at each point.
true
sample mean
Bayesian MCE[ f | x]

0.5

0

−0.5
−2

−1

0
x

1

2

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

18 . 5 . su mma ry 337

18.5

Summary• The expected value and variance of the objective function are useful when
optimizing problems involving uncertainty, but computing these quantities
reliably can be challenging.
• One of the simplest approaches is to estimate the moments using sampling ina process known as Monte Carlo integration.
• Other approaches, such as the Taylor approximation, use knowledge of theobjective function’s partial derivatives.
• Polynomial chaos is a powerful uncertainty propagation technique based onorthogonal polynomials.
• Bayesian Monte Carlo uses Gaussian processes to efficiently arrive at the moments with analytic results for Gaussian kernels.

18.6

Exercises

Exercise 18.1. Suppose we draw a sample from a univariate Gaussian distribution.
What is the probability that our sample falls within one standard deviation of themean (x ∈ [µ − σ, µ + σ])? What is the probability that our sample is less than
one standard deviation above the mean (x < µ + σ)?
Exercise 18.2. Let x (1) , x (2) , . . . , x (m) be a random sample of independent, identically distributed values of size m from a distribution with mean µ and variance ν.
Show that the variance of the sample mean Var(µ̂) is ν/m.
Exercise 18.3. Derive the recurrence relation equation (18.29) that is satisfied by
all orthogonal polynomials.
Exercise 18.4. Suppose we have fitted a polynomial chaos model of an objective function f (x, z) for a particular design point x using m evaluations withz(1) , . . . , z(m) . Derive an expression for estimating the partial derivative of thepolynomial chaos coefficients with respect to a design component xi .

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

338 c ha p te r 18 . u n c e rta i n ty p rop a g a ti on

Exercise 18.5. Consider an objective function f (x, z) with design variables x
and random variables z. As discussed in chapter 17, optimization under uncertainty often involves minimizing a linear combination of the estimated mean andvariance:f mod (x, z) = αµ̂(x) + (1 − α)ν̂(x)How can one use polynomial chaos to estimate the gradient of f mod withrespect to a design variable x?

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com