Monte Carlo Algorithms Where the Integrand Size is Unknown Fred J. - - PowerPoint PPT Presentation

▶
monte carlo algorithms where the integrand size is unknown
SMART_READER_LITE
LIVE PREVIEW

Monte Carlo Algorithms Where the Integrand Size is Unknown Fred J. - - PowerPoint PPT Presentation

Jul 13, 2023 225 likes •640 views

Introduction My cubMC function Examples and Extensions Monte Carlo Algorithms Where the Integrand Size is Unknown Fred J. Hickernell Department of Applied Mathematics Illinois Institute of Technology Email: hickernell@iit.edu , Web:

slide-1
SLIDE 1

Introduction My cubMC function Examples and Extensions

Monte Carlo Algorithms Where the Integrand Size is Unknown

Fred J. Hickernell

Department of Applied Mathematics Illinois Institute of Technology Email: hickernell@iit.edu, Web: www.iit.edu/~hickernell Joint work with Lan Jiang, Yuewei Liu, and Art Owen Thanks to Ian Sloan, Frances Kuo, Josef Dick, and Gareth Peters for inviting me. This work is partially supported by NSF-DMS-1115392 and NSF-DMS-1135257 grants.

  • Feb. 16, 2012

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 1 / 28

slide-2
SLIDE 2

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-3
SLIDE 3

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi), where the Xi are i.i.d. ∼ ρ.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-4
SLIDE 4

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi), where the Xi are i.i.d. ∼ ρ. How large should I make n?

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-5
SLIDE 5

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi), where the Xi are i.i.d. ∼ ρ. How large should I make n? As large as your computational budget allows.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-6
SLIDE 6

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi), where the Xi are i.i.d. ∼ ρ. How large should I make n to obtain |µ − ˆ µ| ≤ ε? The Central Limit Theorem says n = 1.96σ ε 2 where σ2 is the variance of the integrand.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-7
SLIDE 7

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi), where the Xi are i.i.d. ∼ ρ. How large should I make n to obtain |µ − ˆ µ| ≤ ε? The Central Limit Theorem says n = 1.96σ ε 2 where σ2 is the variance of the integrand. How do I find σ2?

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-8
SLIDE 8

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi+nσ), where the Xi are i.i.d. ∼ ρ. How large should I make n to obtain |µ − ˆ µ| ≤ ε? The Central Limit Theorem says n = 1.96ˆ σ ε 2 How do I find σ2? Try the sample variance times a variance inflation factor: ˆ σ2 = C2 nσ − 1

nσ

  • i=1

[f(Xi) − ˆ µσ]2.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-9
SLIDE 9

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi+nσ), where the Xi are i.i.d. ∼ ρ. How large should I make n to obtain |µ − ˆ µ| ≤ ε? The Central Limit Theorem says n = 1.96ˆ σ ε 2 How do I find σ2? Try the sample variance times a variance inflation factor: ˆ σ2 = C2 nσ − 1

nσ

  • i=1

[f(Xi) − ˆ µσ]2. Does theory guarantee that this algo- rithm works (at least 95% of the time)?

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-10
SLIDE 10

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi+nσ), where the Xi are i.i.d. ∼ ρ. How large should I make n to obtain |µ − ˆ µ| ≤ ε? The Central Limit Theorem says n = 1.96ˆ σ ε 2 How do I find σ2? Try the sample variance times a variance inflation factor: ˆ σ2 = C2 nσ − 1

nσ

  • i=1

[f(Xi) − ˆ µσ]2. Does theory guarantee that this algo- rithm works (at least 95% of the time)? Yes! This algorithm, with minor modifi- cations, carries a limited warranty.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 2 / 28

slide-11
SLIDE 11

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Three Perspectives

µ = E[f(X)] =

  • Rd f(x) ρ(x) dx =?

Algorithm Design Construct an automatic multivariate integrator analogous to MATLAB’s quad for univariate integrals. Information-Based Complexity Construct an algorithm, A, satisfying |µ − A(f)| ≤ ǫ (definitely, with high probability, or on average) with cost(ε, A, f) depending reasonably on ε and the unknown size(f). Statistics Find a nonparametric confidence interval of prescribed half-width ε for µ from a reasonable number of samples Yi = f(Xi).

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 3 / 28

slide-12
SLIDE 12

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

MATLAB’s Quadrature Routine quad Works Well, but It Can Be Fooled

0.2 0.4 0.6 0.8 1 0.4 0.6 0.8 1 1.2 1.4 x f(x) integral 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 x f(x) integral data 0.1 0.12 0.14 0.16 0.18 0.2 0.5 1 1.5 2 x f(x) integral

2 √π 1 e−x2 dx = 0.8427007929497149 1 f(x) dx = 1.5436 1 [1+cos(200πx)] dx = 1 quad → 0.8427007929497149 in 0.160521 seconds. but quad → 2 in 0.007092 seconds. but quad → 0.7636784919876782 in 0.205272 seconds.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 4 / 28

slide-13
SLIDE 13

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

Can We Have a

Guarantee Like This? For ni e integrands, f , quad will provide

b

a f(x) dx

with an error ≤ ε in a reasonable amount
  • f
time,
  • r
your money ba k. A ni e integrand, f , satisfies the following
  • nditions:

◮

. . . , i.e., quad won't be fooled,

◮

. . . , i.e., the number
  • f
fun tion values required is moderate. If f is not ni e (nasty), then this guarantee is void, and quad may return an in orre t answer.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 5 / 28

slide-14
SLIDE 14

Introduction My cubMC function Examples and Extensions How to Determine Sample Size MATLAB’s quad function

An

Impra ti al Guarantee For integrands, f , satisfying f ′′∞ ≤ M , a trapezoidal rule with n =
  • (b − a)3M/(12ε)
trapezoids will provide

b

a f(x) dx

with an absolute error ≤ ε .

To apply this guarantee, one must know M in advance, which is impractical. This is why quad (adaptive recursive Simpson’s rule) estimates the error and adaptively determines the number of function evaluations, n. If the algorithm works for f, it should normally work for cf.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 6 / 28

slide-15
SLIDE 15

Introduction My cubMC function Examples and Extensions Guarantee the Variance Estimate Guarantee the Integral (Mean) Guarantee the Cost (Sample Size, Time)

Recall Our Hypothetical Conversation

Practitioner You, the Expert I need to evaluate integrals µ =

  • Rd f(x) ρ(x) dx,

for many different f, where ρ is a given probability density function. Try a sample average, ˆ µ = 1 n

n

  • i=1

f(Xi+nσ), where the Xi are i.i.d. ∼ ρ. How large should I make n to obtain |µ − ˆ µ| ≤ ε? The Central Limit Theorem says n = 1.96ˆ σ ε 2 How do I find σ2? Try the sample variance times a variance inflation factor: ˆ σ2 = C2 nσ − 1

nσ

  • i=1

[f(Xi) − ˆ µσ]2. Does theory guarantee that this algo- rithm works (at least 95% of the time)? Yes! This algorithm, with minor modifi- cations, carries a limited warranty.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 7 / 28

slide-16
SLIDE 16

Introduction My cubMC function Examples and Extensions Guarantee the Variance Estimate Guarantee the Integral (Mean) Guarantee the Cost (Sample Size, Time)

Guarantee the Variance

The sample variance, v is an unbiased estimate of σ2 =

  • Rd[f(x) − µ]2 ρ(x) dx.

v = 1 nσ − 1

nσ

  • i=1

[f(Xi) − ˆ µnσ]2, ˆ µnσ = 1 nσ

nσ

  • i=1

f(Xi), X1, X2, . . . i.i.d. ∼ ρ E[v] = σ2, var(v) = σ4 nσ

  • κ − nσ − 3

nσ − 1

  • ,

κ :=

  • Rd[f(x) − µ]4 ρ(x) dx

σ4 Cantelli’s Inequality (Lin and Bai, 2010, 6.1e) guarantees that an inflated sample variance bounds the variance from above with uncertainty ˜ α, ˆ σ2 := C2v, Prob(ˆ σ2 ≥ σ2) ≥ 1 − ˜ α, C > 1 provided that the kurtosis of the integrand, κ, is not too large, i.e., κ ≤ nσ − 3 nσ − 1 + ˜ αnσ 1 − ˜ α 1 − 1 C2 2 =: κmax(˜ α, nσ, C).

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 8 / 28

slide-17
SLIDE 17

Introduction My cubMC function Examples and Extensions Guarantee the Variance Estimate Guarantee the Integral (Mean) Guarantee the Cost (Sample Size, Time)

Guarantee the Variance

ˆ σ2 = C2 nσ − 1

nσ

  • i=1

[f(Xi) − ˆ µnσ]2, Prob(ˆ σ2 ≥ σ2) ≥ 1 − ˜ α if κ ≤ κmax(˜ α, nσ, C)

10

−3

10

−2

10

−1

10 10 10

1

10

2

10

3

10

4

10

5

κmax ˜ α nσ = 30 nσ = 100 nσ = 1000 nσ = 10000 nσ = 100000

C = 1.5

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 9 / 28

slide-18
SLIDE 18

Introduction My cubMC function Examples and Extensions Guarantee the Variance Estimate Guarantee the Integral (Mean) Guarantee the Cost (Sample Size, Time)

Guarantee the Integral (Mean)

The Central Limit Theorem gives an asymptotic result for fixed z ≥ 0: Prob

  • Rd f(x) ρ(x) dx
  • µ

− 1 n

n

  • i=1

f(Xi+nσ)

  • ˆ

µ

  • ≤ zσ

√n

  • → 1 − 2Φ(−z)

as n → ∞ A non-uniform Berry-Esseen Inequality (Petrov, 1995, Theorem 5.16, p. 168) gives a hard upper bound: Prob

  • |µ − ˆ

µ| ≤ zσ √n

  • ≥ 1 − 2
  • Φ(−z) + 0.56κ3/4

√n (1 + |z|)−3

  • This guarantees that Prob [|µ − ˆ

µ| ≤ ε] ≥ 1 − ˜ α if the sample size is large enough: n ≥ NB(ε/σ, ˜ α, κ) := min

  • m ∈ N : Φ
  • −ε√m/σ
  • +

0.56κ3/4 √m (1 + ε√m/σ)

3 ≤ ˜

α 2

  • ≍ σ2

ε2 as ε σ → 0

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 10 / 28

slide-19
SLIDE 19

Introduction My cubMC function Examples and Extensions Guarantee the Variance Estimate Guarantee the Integral (Mean) Guarantee the Cost (Sample Size, Time)

cubMC

To evaluate µ =

  • Rd f(x) ρ(x) dx

given input f, ρ, ε, α, nσ, C, and Nmax:

◮ Compute ˜

α = 1 − √1 − α, and the maximum kurtosis allowed, κmax(˜ α, nσ, C).

◮ Overestimate the variance: ˆ

σ2 = C2 nσ − 1

nσ

  • i=1

[f(Xi) − ˆ µnσ]2.

◮ Choose the new sample size, n = min (max (nσ, NB(ε/ˆ

σ, ˜ α, κmax)) , Nmax), for the sample mean.

◮ Finally, compute the sample mean: ˆ

µ = 1 n

n

  • i=1

f(Xi+nσ). Then Prob [|µ − ˆ µ| ≤ ε] ≥ 1 − α provided κ ≤ κmax and n < Nmax.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 11 / 28

slide-20
SLIDE 20

Introduction My cubMC function Examples and Extensions Guarantee the Variance Estimate Guarantee the Integral (Mean) Guarantee the Cost (Sample Size, Time)

Guarantee the Time (Sample Size)

Cantelli’s inequality also tells us that the estimated variance, ˆ σ2, will not overestimate the true variance, σ2, by much, and so the number of function values needed is not unnecessarily large: cost(ε, cubMC, σ) = sup

f:κ≤κmax

min

N {Prob[nσ + n ≤ N] ≥ 1 − β}

≤ nσ + max(nσ, NB(ε/(σγ), ˜ α, κ3/4

max)) ≍ σ2

ε2 , γ := C   1 +

  • ˜

α 1 − ˜ α 1 − β β 1 − 1 C2 2   

1/2

. Cost depends on σ2 = var(f), but the algorithm does not need to know σ2.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 12 / 28

slide-21
SLIDE 21

Introduction My cubMC function Examples and Extensions Guarantee the Variance Estimate Guarantee the Integral (Mean) Guarantee the Cost (Sample Size, Time)

Guarantee for cubMC For ni e integrands cubMC will provide the value
  • f µ =
  • Rd f(x) ρ(x) dx
with an absolute error
  • f

≤ ε ,

with probability 1 − α , in time ≍ (σ/ε)2 with probability 1 − β ,
  • r
your money ba k. A ni e integrand, f , satisfies the following
  • nditions:

◮

the kurtosis is not too large, i.e., κ ≤ κmax(˜

α, nσ, C)

, and

◮

the varian e is not
  • verwhelming,
i.e., σ2 ≤ cε2Nmax/d , where Nmax is the maximum number
  • f
s alar samples. If f is not ni e (nasty), cubMC may return the wrong answer.

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 13 / 28

slide-22
SLIDE 22

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Peak Function — cubMC

f(x) =    1 + σ

  • 1−p

p ,

0 ≤ x − z (mod 1) ≤ p, 1 − σ

  • p

1−p,

p < x − z (mod 1) ≤ 1, µ = 1, κ = 1 p(1 − p) − 3

0.2 0.4 0.6 0.8 1 2 4 6 8 10 12 x f(x) integral

z ∼ U(0, 1) p ∈

  • 10−5, 1/2
  • ,

σ ∈ [0.1, 10] log(p), log(σ) ∼ Uniform α = 5%, C = 1.5, nσ = 1024, ε = 0.001 κmax = 9.2, Nmax = 109 covered by guarantee kurtosis too large truncated sample

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 14 / 28

slide-23
SLIDE 23

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Peak Function — cubMC vs. quad & quadgk

           MATLAB’s quad ε = 0.001 fast tolerance rarely met no guarantee

  • my cubMC

ε = 0.001 covered by guarantee kurtosis too large truncated sample            MATLAB’s quadgk ε = 0.001 fast tolerance rarely met no guarantee

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 15 / 28

slide-24
SLIDE 24

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Peak Function — cubMC i.i.d. vs. Sobol’, Plus More Robustness

       i.i.d. sampling covered by guarantee kurtosis too large truncated sample    Sobol’ sampling no guarantee yet faster nσ = 1024, κmax = 9.2 nσ = 131072, κmax = 1050 ε = 0.001

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 16 / 28

slide-25
SLIDE 25

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Peak Function for d = 3

       i.i.d. sampling covered by guarantee kurtosis too large truncated sample        Sobol’ sampling quasi-standard error no guarantee yet faster nσ = 1024, κmax = 9.2 nσ = 131072, κmax = 1050 ε = 0.001

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 17 / 28

slide-26
SLIDE 26

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Asian Geometric Mean Call, d = 1, 2, 4, . . . , 64

  • i.i.d. sampling
  • Sobol’ sampling

nσ = 1024, κmax = 9.2 nσ = 131072, κmax = 1050 ε = 0.02

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 18 / 28

slide-27
SLIDE 27

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Why an Adaptive Algorithm?

cubMC non-adaptive algorithm

10

3

10

7

10

11

10

3

10

7

10

11

G u a r a n t e e d T

  • B

i g N

  • t

G u a r a n t e e d N a s t y I m p

  • s

s i b l e f − µ2

2/ε2 = σ 2/ε2 ∝ time

f − µ2

4/ε2 = √κσ 2/ε2

2 min 0.01 sec 10

3

10

7

10

11

10

3

10

7

10

11

G u a r a n t e e d T

  • B

i g I m p

  • s

s i b l e f − µ2

2/ε2 = σ 2/ε2 ∝ time

f − µ2

4/ε2 = √κσ 2/ε2

2 min 0.01 sec ◮ Cost depends on size of integrand ◮ Algorithm parameters determine

robustness to nasty integrands

◮ Tiny integrands handled regardless ◮ Huge integrands cannot be handled ◮ Cost is fixed and high if you want it

reach tolerance for lots of integrands

◮ Huge integrands cannot be handled

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 19 / 28

slide-28
SLIDE 28

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Why Make cubMC Dependent on Kurtosis?

σ = difficulty κ = nastiness The kurtosis of a random variable or function, κ = E[(Y − µ)4] {E[(Y − µ)2]}2

  • σ4
  • r

κ :=

  • Rd[f(x) − µ]4 ρ(x) dx
  • Rd[f(x) − µ]2 ρ(x) dx

2

  • σ4

is difficult to estimate. Why should cubMC’s guarantee depend on bounded κ?

◮ Practically, we need κ bounded to justify the estimates of σ2. ◮ Bounded κ yields sets of probability distributions or functions that are non-convex.

◮ Nonparametric confidence intervals are impossible for convex sets of distributions

(Bahadur and Savage, 1956, Corollary 2).

How we break convexity ◮ Adaptive information does not help for convex sets of integrands in the worst case and

probabilistic settings (Traub et al, 1988, Chapter 4, Theorem 5.2.1; Chapter 8, Corollary 5.3.1).

How we break convexity hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 20 / 28

slide-29
SLIDE 29

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Quasi-Standard Error (Internal Replications) for Quasi-Random Sequences

Let X1, X2, . . . be a (random or deterministic) sequence, let r be fixed, and let ˆ µm = 1 2m

2m

  • i=1

f(Xi) = 1 2r

2r

  • j=1

ˆ µm,j, ˆ µm,j = 1 2m−r

2m−r

  • i=1

f(X(j−1)2m−r+i) The quasi-standard error (Owen, 1997) measures the variation of among the means of parts of the whole sample qsem =

  • 1

2r(2r − 1)

2r

  • j=1

(ˆ µm,j − ˆ µm)2 Given error tolerance, ε, and parameters r ∈ N, m1 ≥ r, and C > 1, for m = m1, m1 + 1, . . .,

◮ Compute f(X2m−1+1), . . . , f(X2m), and ˆ

µm,1, . . . , ˆ µm,2r, ˆ µm.

◮ If C qsem ≤ ε, then stop. Else continue.

When this works hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 21 / 28

slide-30
SLIDE 30

Introduction My cubMC function Examples and Extensions Numerical examples Discussion Summary & Further Work

Further Work

Quasi-Monte Carlo Sampling — What is a good measure of an integrand being nasty or nice? Variance Reduction Techniques — Can we preserve the guarantee? Different Error Criteria — Worst case? Randomized? Lower Bounds on Cost — The typical fooling functions are nasty (high kurtosis). Does assuming moderate kurtosis make the problem easier? Relative Error Tolerances — Both the variance and the mean are needed to determine the eventual sample size. Unbounded or Infinite d — Can automatic integrators for finite d be used in multilevel methods to improve efficiency?

Look here

Other Problems — Are there any guarantees for MATLAB’s quad, or any other univariate adaptive quadrature routine that estimates error? What about guarantees for function approximation?

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 22 / 28

slide-31
SLIDE 31

References Extra Material

References I

Bahadur RR, Savage LJ (1956) The nonexistence of certain statistical procedures in nonparametric problems. Ann Math Stat 27:1115–1122 Halton JH (2005) Quasi-probability: Why quasi-Monte-Carlo methods are statistically valid and how their errors can be estimated statistically. Monte Carlo Methods and Appl 11:203–350 Lin Z, Bai Z (2010) Probability Inequalities. Science Press and Springer-Verlag, Beijing and Berlin Owen AB (1997) Scrambled net variance for integrals of smooth functions. Ann Stat 25:1541–1562 Owen AB (2006) On the Warnock-Halton quasi-standard error. Monte Carlo Methods and Appl 12:47–54 Petrov VV (1995) Limit Theorems of Probability Theory:Sequences of Independent Random Variables. Clarendon Press, Oxford Snyder WC (2000) Accuracy estimation for quasi-Monte Carlo simulations. Math Comput Simulation 54:131–143 Traub JF, Wasilkowski GW, Wo´ zniakowski H (1988) Information-Based Complexity. Academic Press, Boston

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 23 / 28

slide-32
SLIDE 32

References Extra Material Breaking Convexity When Quasi-Standard Error Works

A Set of Distributions with Bounded Kurtosis is Non-Convex

Prob(Y1 = y) = 0.5, y = ±1 Prob(Y2 = y) = 0.5, y = ±2 f3 = 1 2f1 + 1 2f2 Prob(Y3 = y) = 0.25, y = ±1, ±2 κ3 = 34 25 > 1 = κ1 = κ2

−3 −2 −1 1 2 3 0.5 1 f1 κ=1 −3 −2 −1 1 2 3 0.5 1 f2 κ=1 −3 −2 −1 1 2 3 0.5 1 f3 κ=34/25

Return hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 24 / 28

slide-33
SLIDE 33

References Extra Material Breaking Convexity When Quasi-Standard Error Works

A Set of Integrands with Bounded Kurtosis is Non-Convex

f1(x) =

  • −1,

0 ≤ x < 1/2 1, 1/2 ≤ x ≤ 1 f2(x) =      1, 0 ≤ x < 1/4 −1, 1/4 ≤ x ≤ 3/4 1, 3/4 ≤ x ≤ 1 f3(x) = 1

2f1(x) + 1 2f2(x)

=          0, 0 ≤ x < 1/4, −1, 1/4 ≤ x ≤ 1/2, 0, 1/2 ≤ x ≤ 3/4, 1, 3/4 ≤ x ≤ 1, κ3 = 2 > 1 = κ1 = κ2

0.2 0.4 0.6 0.8 1 −1 1 f1 κ=1 0.2 0.4 0.6 0.8 1 −1 1 f2 κ=1 0.2 0.4 0.6 0.8 1 −1 1 f3 κ=2

Return hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 25 / 28

slide-34
SLIDE 34

References Extra Material Breaking Convexity When Quasi-Standard Error Works

When Does Quasi-Standard Error Work?

Quasi-standard error has been proposed by Warnock and studied by Owen (1997); Snyder (2000); Halton (2005); Owen (2006). Suppose that X1, X2, . . . is a scrambled digital (t, d)-sequence in base 2, Pm = {X1, . . . , X2m}, and the integrand can be expanded in a Walsh series: f(x) =

  • k∈Nd

ˆ f(k)eπ√−1k⊗x, (SBf)(x) :=

  • k∈B

ˆ f(k)eπ√−1k⊗x (filtered f) µ =

  • [0,1)d f(x) dx = S{0}f,

ˆ µm = 1 2m

2m

  • i=1

f(Xi) = (SP ⊥

mf)(X1),

where ⊗ is a bitwise dot product modulo 2, and P ⊥

m = {k ∈ Nd 0 : k ⊗ x = 0 ∀x ∈ Pm} is

the dual net (wavenumbers aliased with 0). The quasi-standard error may be expressed as qsem =

  • 1

2r − 1

  • l∈(P ⊥

m−r/P ⊥ m)\{0}

(SP ⊥

m⊕lf)2(X1)

which is a surrogate for |µ − ˆ µm| =

  • (SP ⊥

m\{0}f)(X1)

  • .

hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 26 / 28

slide-35
SLIDE 35

References Extra Material Breaking Convexity When Quasi-Standard Error Works

When Does Quasi-Standard Error Work?

f(x) =

  • k∈Nd

ˆ f(k)eπ√−1k⊗x |µ − ˆ µm| =

  • (SP ⊥

m\{0}f)(X1)

  • P ⊥

6 \ {0} 1 10 100 1 10 100 k1 k2

qsem =

  • 1

7

  • l∈(P ⊥

m−3/P ⊥ m)\{0}

(SP ⊥

m⊕lf)2(X1)

  • P ⊥

3 \ P ⊥ 6

Return hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 27 / 28

slide-36
SLIDE 36

References Extra Material Breaking Convexity When Quasi-Standard Error Works

When Does Quasi-Standard Error Work?

f(x) =

  • k∈Nd

ˆ f(k)eπ√−1k⊗x |µ − ˆ µm| =

  • (SP ⊥

m\{0}f)(X1)

  • P ⊥

7 \ {0} 1 10 100 1 10 100 k1 k2

qsem =

  • 1

7

  • l∈(P ⊥

m−3/P ⊥ m)\{0}

(SP ⊥

m⊕lf)2(X1)

  • P ⊥

4 \ P ⊥ 7

Return hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 27 / 28

slide-37
SLIDE 37

References Extra Material Breaking Convexity When Quasi-Standard Error Works

When Does Quasi-Standard Error Work?

f(x) =

  • k∈Nd

ˆ f(k)eπ√−1k⊗x |µ − ˆ µm| =

  • (SP ⊥

m\{0}f)(X1)

  • P ⊥

8 \ {0} 1 10 100 1 10 100 k1 k2

qsem =

  • 1

7

  • l∈(P ⊥

m−3/P ⊥ m)\{0}

(SP ⊥

m⊕lf)2(X1)

  • P ⊥

5 \ P ⊥ 8

Return hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 27 / 28

slide-38
SLIDE 38

References Extra Material Breaking Convexity When Quasi-Standard Error Works

When Does Quasi-Standard Error Work?

f(x) =

  • k∈Nd

ˆ f(k)eπ√−1k⊗x |µ − ˆ µm| =

  • (SP ⊥

m\{0}f)(X1)

  • P ⊥

9 \ {0} 1 10 100 1 10 100 k1 k2

qsem =

  • 1

7

  • l∈(P ⊥

m−3/P ⊥ m)\{0}

(SP ⊥

m⊕lf)2(X1)

  • P ⊥

6 \ P ⊥ 9

Return hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 27 / 28

slide-39
SLIDE 39

References Extra Material Breaking Convexity When Quasi-Standard Error Works

When Does Quasi-Standard Error Work?

f(x) =

  • k∈Nd

ˆ f(k)eπ√−1k⊗x |µ − ˆ µm| =

  • (SP ⊥

m\{0}f)(X1)

  • P ⊥

10 \ {0} 1 10 100 1 10 100 k1 k2

qsem =

  • 1

7

  • l∈(P ⊥

m−3/P ⊥ m)\{0}

(SP ⊥

m⊕lf)2(X1)

  • P ⊥

7 \ P ⊥ 10

Return hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 27 / 28

slide-40
SLIDE 40

References Extra Material Breaking Convexity When Quasi-Standard Error Works

Asian Geometric Mean Call Execution Times

10 10

1

10

2

10

−3

10

−2

10

−1

10 10

1

10

2

d Time (seconds) 10 10

1

10

2

10

−3

10

−2

10

−1

10 10

1

10

2

d Time (seconds)

  • i.i.d. sampling

10 10

1

10

2

10

−3

10

−2

10

−1

10 10

1

10

2

d Time (seconds) 10 10

1

10

2

10

−3

10

−2

10

−1

10 10

1

10

2

d Time (seconds)

  • Sobol’ sampling

nσ = 1024, κmax = 9.2 nσ = 131072, κmax = 1050

Return hickernell@iit.edu Adaptive Monte Carlo MCQMC 2012 28 / 28