[PPT] - Quantification and Density Estimation Gadi Fibich, Tel Aviv PowerPoint Presentation

SLIDE 1

1

Adi Ditkowski
Amir Sagiv

A Spline-based approach to Uncertainty Quantification and Density Estimation

Gadi Fibich, Tel Aviv University

SLIDE 2

2

Motivation

SLIDE 3

3

Input pulse characteristics

Elliptic

typical shot average over 1000 shots

SLIDE 4

4

Input pulse characteristics

Elliptic
Varies from shot to shot
Always the case

average over 1000 shots typical shot

SLIDE 5

5

After 5 meters in air

Average

ver 1000

shots typical

SLIDE 6

6

After 5 meters in air

Average

ver 1000

shots typical

SLIDE 7

7 Average

ver 1000

shots

After 5 meters in air

typical

Want to predict

Average intensity after

5m

Probability for 0,1,2, or

3 filaments

SLIDE 8

𝜔 𝑨, 𝑦, 𝑧

utput

𝜔0(𝑦, 𝑧)

NLS initial condition

8

Mathematical model

SLIDE 9

𝜔 𝑨, 𝑦, 𝑧; 𝛽

random output

𝜔0(𝑦, 𝑧; 𝛽)

𝛽 - noise parameter

random initial condition

9

Shot to shot variation

NLS

SLIDE 10

𝜔 𝑨, 𝑦, 𝑧; 𝛽

random output

𝜔0(𝑦, 𝑧; 𝛽)

𝛽 - noise parameter

random initial condition

10

Shot to shot variation

NLS

Computational goals Moment estimation

e.g., average intensity 𝐅𝜷 𝜔 2

Density estimation

e.g., probability for 2 filaments

SLIDE 11

Nonlinear initial value problem ቊ 𝑣𝑢 𝑢, 𝒚 = 𝑅 𝒚, 𝑣 𝑣 𝑣 𝑢 = 0, 𝒚 = 𝑣0(𝒚)

“Quantity of interest” (model output) 𝑔 = 𝑔[𝑣]
e.g., 𝑔 = arg(𝑣 𝑢𝑗, 𝑦𝑗 ), f = ∫ u 2𝑒𝑦, …
𝑣, 𝑔[𝑣] not given explicitly, but can be evaluated numerically

11

General setting

SLIDE 12

General setting with randomness

Add randomness (in 𝑣0 and/or 𝑅) ቊ 𝑣𝑢 𝑢, 𝒚; 𝜷 = 𝑅 𝒚, 𝑣; 𝜷 𝑣 𝑣 𝑢 = 0, 𝒚; 𝜷 = 𝑣0(𝒚; 𝜷)

𝜷 distributed according to a known measure
“Quantity of interest” (model output) 𝑔 𝜷 ≔ 𝑔[𝑣 𝑢, 𝒚; 𝜷 ]
𝑔 𝜷 is a random variable

12

SLIDE 13

Add randomness (in 𝑣0 and/or 𝑅) ቊ 𝑣𝑢 𝑢, 𝒚; 𝜷 = 𝑅 𝒚, 𝑣; 𝜷 𝑣 𝑣 𝑢 = 0, 𝒚; 𝜷 = 𝑣0(𝒚; 𝜷)

𝜷 distributed according to a known measure
“Quantity of interest” (model output) 𝑔 𝜷 ≔ 𝑔[𝑣 𝑢, 𝒚; 𝜷 ]
𝑔 𝜷 is a random variable

13

Computational goals:

Moment estimation

𝐅𝜷 𝑔

Density estimation

Probability Density Function (PDF) of 𝑔 𝜷

General setting with randomness

SLIDE 14

Standard statistical methods

Moment estimation

Monte-Carlo 𝐅𝛽 𝑔 ≈

1 𝑂 σ𝑜=1 𝑂

𝑔

𝑜

…

14

Density (PDF) estimation

Histogram method
Kernel density estimators (KDE)
…

Step I – draw samples {𝛽1, … , 𝛽𝑂} Step II – compute 𝑔

1, … , 𝑔 𝑂 ,

𝑔

𝑜:=𝑔 𝛽𝑜

SLIDE 15

Standard statistical methods

Moment estimation

Monte-Carlo 𝐅𝛽 𝑔 ≈

1 𝑂 σ𝑜=1 𝑂

𝑔

𝑜

…

15

Density (PDF) estimation

Histogram method
Kernel density estimators (KDE)
…

Step I – draw samples {𝛽1, … , 𝛽𝑂} Step II – compute 𝑔

1, … , 𝑔 𝑂 ,

𝑔

𝑜:=𝑔 𝛽𝑜

Constraint:

Computation of 𝑔 𝜷𝒌 is expensive (e.g., solving the (3+1)D NLS)
Can only use a small samples {𝑔 𝛽1 , … , 𝑔 𝛽𝑂 }

SLIDE 16

Standard statistical methods

16

Poor approximations for small N

Moment estimation

Monte-Carlo 𝐅𝛽 𝑔 ≈

1 𝑂 σ𝑜=1 𝑂

𝑔

𝑜

…

Density (PDF) estimation

Histogram method
Kernel density estimators (KDE)
…

e.g. Histogram method with N=10 samples Step I – draw samples {𝛽1, … , 𝛽𝑂} Step II – compute 𝑔

1, … , 𝑔 𝑂 ,

𝑔

𝑜:=𝑔 𝛽𝑜

SLIDE 17

Standard statistical methods

17

Poor approximations for small N

Moment estimation

Monte-Carlo 𝐅𝛽 𝑔 ≈

1 𝑂 σ𝑜=1 𝑂

𝑔

𝑜

…

Density (PDF) estimation

Histogram method
Kernel density estimators (KDE)
…

e.g. Histogram method with N=10 samples

17

Exact PDF Step I – draw samples {𝛽1, … , 𝛽𝑂} Step II – compute 𝑔

1, … , 𝑔 𝑂 ,

𝑔

𝑜:=𝑔 𝛽𝑜

SLIDE 18

Standard statistical methods

Moment estimation

Monte-Carlo 𝐅𝛽 𝑔 ≈

1 𝑂 σ𝑜=1 𝑂

𝑔

𝑜

…

18

Density (PDF) estimation

Histogram method
Kernel density estimators (KDE)
…

Given a sample {𝑔

1, … , 𝑔 𝑂} of 𝑔 𝛽

How to improve?

Above methods only use {𝑔

1, … , 𝑔 𝑂}

Uncertainty Quantification (UQ) approach: Utilize

1. The relation 𝑔 𝛽 2. Smoothness of 𝑔 𝛽

SLIDE 19

Approximation-based estimation

19

𝑔 𝛽 𝑔

𝑂 𝛽

Questions

How accurate are E𝛽 𝑔 − E𝛽 𝑔

𝑂 , ||𝑞 − 𝑞𝑂||?

Which approximation method should be used?

approximation moment, PDF moment, PDF

𝐅𝛽 𝑔 , 𝑞 𝐅𝛽 𝑔

𝑂 , 𝑞𝑂

𝑜 → ∞

can take a large sample, since computation of 𝑔

𝑂 𝛽 is cheap

𝑞 and 𝑞𝑂 are the PDFs of 𝑔 and 𝑞𝑂
𝑞 is the PDF of 𝑔

SLIDE 20

Approximation-based estimation

20

𝑔 𝛽 𝑔

𝑂 𝛽

Questions

How accurate are E𝛽 𝑔 − E𝛽 𝑔

𝑂 , ||𝑞 − 𝑞𝑂||?

Which approximation method should be used?

approximation moment, PDF moment, PDF

𝐅𝛽 𝑔 , 𝑞 𝐅𝛽 𝑔

𝑂 , 𝑞𝑂

𝑜 → ∞

can take a large sample, since computation of 𝑔

𝑂 𝛽 is cheap

cannot take a large sample

𝑞 and 𝑞𝑂 are the PDFs of 𝑔 and 𝑞𝑂
𝑞 is the PDF of 𝑔

SLIDE 21

Approximation-based estimation

21

𝑔 𝛽 𝑔

𝑂 𝛽

Questions

Which approximation should be used?
How small are E𝛽 𝑔 − E𝛽 𝑔

𝑂

and ||𝑞 − 𝑞𝑂|| ?

approximation moment, PDF moment, PDF

𝐅𝛽 𝑔 , 𝑞 𝐅𝛽 𝑔

𝑂 , 𝑞𝑂

𝑂 → ∞

can take a large sample, since computation of 𝑔

𝑂 𝛽 is cheap

cannot take a large sample

𝑞 is the PDF of 𝑔
𝑞𝑂 is the PDF of 𝑔

𝑂

SLIDE 22

Noise dimension

One-dimensional noise 𝛽 ∈ 𝑆
Random input power 𝜔0 = 1 + 𝛽 𝑓

− 𝑠2

Random temperature
…
Multi-dimensional noise 𝛽 ∈ 𝑆𝑒
Random input power and incidence angle
…

22

SLIDE 23

Noise dimension

One-dimensional noise 𝛽 ∈ 𝑆
Random input power 𝜔0 = 1 + 𝛽 𝑓

− 𝑠2

Random temperature
…
Multi-dimensional noise 𝛽 ∈ 𝑆𝑒
Random input power and incidence angle
…

23

SLIDE 24

Standard uncertainty quantification approach:

Approximate 𝑔 using orthogonal polynomials 𝑟𝑜(𝛽)

𝑔

𝑂 𝛽 = ෍ 𝑜=0 𝑂−1

𝑟𝑜 , 𝑔 𝑟𝑜 𝛽

Spectral accuracy for moments

𝐅𝛽 𝑔 − 𝐅𝛽 𝑔

𝑂 = 𝑃 𝑓−𝛿𝑂 ,

𝑂 ≫ 1 if f is analytic

Generalized Polynomial Chaos (gPC)

24

SLIDE 25

Standard uncertainty quantification approach:

Approximate 𝑔 using orthogonal polynomials 𝑟𝑜(𝛽)

𝑔

𝑂 𝛽 = ෍ 𝑜=0 𝑂−1

𝑟𝑜 , 𝑔 𝑟𝑜 𝛽

Spectral accuracy for moments

𝐅𝛽 𝑔 − 𝐅𝛽 𝑔

𝑂 = 𝑃 𝑓−𝛿𝑂 ,

𝑂 ≫ 1 if f is analytic

Problem solved

Generalized Polynomial Chaos (gPC)

25

SLIDE 26

Standard uncertainty quantification approach:

Approximate 𝑔 using orthogonal polynomials 𝑟𝑜(𝛽)

𝑔

𝑂 𝛽 = ෍ 𝑜=0 𝑂−1

𝑟𝑜 , 𝑔 𝑟𝑜 𝛽

Spectral accuracy for moments

𝐅𝛽 𝑔 − 𝐅𝛽 𝑔

𝑂 = 𝑃 𝑓−𝛿𝑂 ,

𝑂 ≫ 1 if f is analytic

Problem solved

But,

Generalized Polynomial Chaos (gPC)

26

Moment estimation Spectral accuracy reached only for large N How to achieve ``good’’ accuracy with e.g. 𝑶 = 𝟐𝟏 samples?

Density estimation No theory for 𝑞 − 𝑞𝑂 Will it work in practice?

SLIDE 27

Example: Density estimation with gPC

27

𝛽

𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

SLIDE 28

28

𝛽

𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

PDF approximation, 𝑂 = 12

Example: Density estimation with gPC

SLIDE 29

29

𝛽

𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

PDF approximation, 𝑂 = 12

Example: Density estimation with gPC

Lemma 𝑞 𝑧 = ෍

𝑔 𝛽 =𝑧

1 𝑔′ 𝛽

SLIDE 30

30

𝛽

𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

PDF approximation, 𝑂 = 12

Example: Density estimation with gPC

𝛽 f′(α)

Lemma 𝑞 𝑧 = ෍

𝑔 𝛽 =𝑧

1 𝑔′ 𝛽

SLIDE 31

31

𝛽

𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

PDF approximation, 𝑂 = 12

Example: Density estimation with gPC

𝛽 f′(α)

Lemma 𝑞 𝑧 = ෍

𝑔 𝛽 =𝑧

1 𝑔′ 𝛽

SLIDE 32

32

𝛽

𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

PDF approximation, 𝑂 = 12

Example: Density estimation with gPC

𝛽 f′(α)

Conclusion Although gPC is spectrally accurate in 𝑀2, it produces “artificial” zero derivatives which ``contaminate’’ the PDF

Lemma 𝑞 𝑧 = ෍

𝑔 𝛽 =𝑧

1 𝑔′ 𝛽

SLIDE 33

𝑞𝑂 is the PDF of 𝑔

𝑂

Approximation-based estimation

33

𝑔 𝛽 𝑔

𝑂 𝛽

Question: Which approximation should be used?

moment, PDF moment, PDF

𝐅𝛽 𝑔 , 𝑞 𝐅𝛽 𝑔

𝑂 , 𝑞𝑂

𝑂 → ∞

can take a large sample, since computation of 𝑔

𝑂 𝛽 is cheap

cannot take a large sample

𝑞 is the PDF of 𝑔

Answer

For density approximation, require that 𝑔

𝑂 ′ = 0 ⇔ 𝑔′ = 0

“Monotonicity-preserving” approximation

approximation

SLIDE 34

Ditkowski, Fibich, Sagiv, 18: Approximate f using a cubic spline over N grid points Thm (Ditkowski, Fibich, Sagiv, 18) : Let 𝑞 and 𝑞𝑂 denote the PDFs

f 𝑔(𝛽) and its cubic spline interpolant over 𝑂 points. Then

𝑞 − 𝑞𝑂 1 ≤ 𝐷𝑂−3, 𝐷 = 𝑃(1)

No equivalent thm for gPC

Thm: 𝐅𝛽 𝑔 − 𝐅𝛽 𝑔

𝑂

≤ 𝐷𝑂−4, 𝐷 = 𝑃(1)

Worse than gPC for large N
But, usually better than gPC for small N

Adopt a spline-based approach

34

SLIDE 35

35

𝑂

𝛽

Example – moment estimation

𝜏 𝑔 − 𝜏 𝑔

𝑂

Moment estimation 𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

𝑂

SLIDE 36

36

𝑂

𝛽

Example – moment estimation

𝜏 𝑔 − 𝜏 𝑔

𝑂

𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

Moment estimation Spline better than gPC for small N

SLIDE 37

𝜏 𝑔 − 𝜏 𝑔

𝑂

PDF estimation

37

𝛽

PDF approximation, 𝑂 = 12

Statistically optimal 𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

SLIDE 38

𝜏 𝑔 − 𝜏 𝑔

𝑂

PDF estimation

38

𝛽

PDF approximation, 𝑂 = 12

𝑔 = tanh 9𝛽 +

𝛽 2 , 𝛽 ∼ Uniform [−1, 1]

SLIDE 39

alpha_min = −1; alpha_max = 1 ; N = 18; f = @(x) tanh (9∗x) + x/2; samplingGrid = linspace(alpha_min, alpha_max, N) ; sample_s = f (samplingGrid); M = 2e6 ; denseGrid = linspace (alpha_min, alpha_max ,M) ; fNspline = spline ( samplingG rid, samples, denseGrid) Cf = 1 . 6 9 ; L =Cf∗Mˆ (1/3 ) ; [histogram, binsEdges ] = hist( fNspline ,L) ; binWidth = (max( binsEdges)−min (binsEdges)) /L; pdf = histogram / (sum(histogram) ∗binWidth ) ; plot(binsEdges, pdf )

Matlab code for PDF estimation

39

SLIDE 40

𝜔 𝑨, 𝑦, 𝑧; 𝛽

random output

𝜔0(𝑦, 𝑧; 𝛽)

𝛽 - noise parameter

random initial condition

40

Shot to shot variation

NLS

Loss of Phase Lemma (Sagiv, Ditkowski, Fibich, 2017) Let 𝑔 𝑨; 𝛽 ≔ arg 𝜔(𝑨, 𝑦 = 0, 𝑧 = 0; 𝛽)𝑛𝑝𝑒(2𝜌). Then lim

𝑨→∞ 𝑔 𝑨; 𝛽

∼ 𝑉(0,2𝜌)

SLIDE 41

Example: PDF of on-axis phase

Histogram 𝑨1 𝑨2 𝑨3

𝑨2 𝑨1 𝑦 𝑨3

Spline

𝑔 PDF PDF

Use N=10 NLS simulations 𝑔 𝛽 = arg 𝜔 𝑨, 0; 𝛽 mod 2𝜌

𝑔 𝑔 𝑔 𝑔 𝑔

SLIDE 42

Coupled NLS – loss of polarization angle

phase: 𝜒± 𝑢 = arg 𝐵± 𝑢, 𝑦 = 0

𝑛𝑝𝑒 (2𝜌)

polarization 𝜄 𝑢 = 𝜒+ 𝑢 − 𝜒− 𝑢 Random elliptical beam – 𝐵± 𝑢 = 0 = 1 + 𝛽 𝐷±𝑓−𝑦2, 𝛽 ∼ 𝑉(−0.1,0.1)

𝒋 𝝐 𝝐𝒖 𝑩± 𝒖, 𝒚 + 𝝐𝟑 𝝐𝒚𝟑 𝑩± + 𝟑 𝟒 𝑩±

𝟑 + 𝟑 𝑩∓ 𝟑 𝑩± = 𝟏

𝜄 𝜄 PDF, N=64 PDF, N=64

Patwardhan et al., 2018

SLIDE 43

Coupled NLS – loss of polarization angle

𝑂 𝑞 − 𝑞𝑂 1

𝒋 𝝐 𝝐𝒖 𝑩± 𝒖, 𝒚 + 𝝐𝟑 𝝐𝒚𝟑 𝑩± + 𝟑 𝟒 𝑩±

𝟑 + 𝟑 𝑩∓ 𝟑 𝑩± = 𝟏

𝜄 𝜄 PDF, N=64 PDF, N=64

phase: 𝜒± 𝑢 = arg 𝐵± 𝑢, 𝑦 = 0

𝑛𝑝𝑒 (2𝜌)

polarization 𝜄 𝑢 = 𝜒+ 𝑢 − 𝜒− 𝑢 Random elliptical beam – 𝐵± 𝑢 = 0 = 1 + 𝛽 𝐷±𝑓−𝑦2, 𝛽 ∼ 𝑉(−0.1,0.1)

𝒋 𝝐 𝝐𝒖 𝑩± 𝒖, 𝒚 + 𝝐𝟑 𝝐𝒚𝟑 𝑩± + 𝟑 𝟒 𝑩±

𝟑 + 𝟑 𝑩∓ 𝟑 𝑩± = 𝟏

Patwardhan et al., 2018

~𝑂−3.7

(𝑢ℎ𝑓𝑝𝑠𝑧: 𝑂−3)

SLIDE 44

Burgers equation – shock location

Initial condition: u0 x = 𝛽 sin(𝑦) Shock location at 𝑢 → ∞

𝛽 = − 𝑑𝑝𝑡 𝑌𝑡

Distribution of random initial amplitude –

𝛽 𝜉 = ൞ −1 + 1 + 4𝜉2 2𝜉 𝜉 ≠ 0 𝜉 = 0 𝜉 ∼ 𝑂(0, 𝜏) 44

𝒗𝒖 𝒖, 𝒚 + 𝟐 𝟑 𝒗𝟑

𝒚 = 𝟐

𝟑 𝒕𝒋𝒐 𝒚

𝒚

𝑌𝑡 𝑌𝑡 PDF, N=7 PDF, N=7

Chen, Gottlieb, Hesthaven, JCP 2005

SLIDE 45

Burgers equation – shock location

Initial condition: u0 x = 𝛽 sin(𝑦) Shock location at 𝑢 → ∞

𝛽 = − 𝑑𝑝𝑡 𝑌𝑡

Distribution of random initial amplitude –

𝛽 𝜉 = ൞ −1 + 1 + 4𝜉2 2𝜉 𝜉 ≠ 0 𝜉 = 0 𝜉 ∼ 𝑂(0, 𝜏) 45

𝑂 𝑞 − 𝑞𝑂 1

𝒗𝒖 𝒖, 𝒚 + 𝟐 𝟑 𝒗𝟑

𝒚 = 𝟐

𝟑 𝒕𝒋𝒐 𝒚

𝒚

𝑌𝑡 𝑌𝑡 PDF, N=7 PDF, N=7

~𝑂−3.1

Chen, Gottlieb, Hesthaven, JCP 2005

(𝑢ℎ𝑓𝑝𝑠𝑧: 𝑂−3)

SLIDE 46

Noise dimension

One-dimensional noise 𝛽 ∈ 𝑆
Random input power 𝜔0 = 1 + 𝛽 𝑓

− 𝑠2

Random temperature
…
Multi-dimensional noise 𝛽 ∈ 𝑆𝑒
Random input power and incidence angle
…

46

SLIDE 47

𝑞𝑂 is the PDF of 𝑔

𝑂

Approximation-based estimation

47

𝑔 𝛽 𝑔

𝑂 𝛽

Questions

Which approximation should be used?
How small are E𝛽 𝑔 − E𝛽 𝑔

𝑂

and ||𝑞 − 𝑞𝑂|| ?

moment, PDF moment, PDF

𝐅𝛽 𝑔 , 𝑞 𝐅𝛽 𝑔

𝑂 , 𝑞𝑂

𝑂 → ∞

can take a large sample, since computation of 𝑔

𝑂 𝛽 is cheap

cannot take a large sample

𝑞 is the PDF of 𝑔

approximation

SLIDE 48

Ditkowski, Fibich, Sagiv, 18: If 𝛽 is d-dimensional, approximate f with a tensor product cubic spline over N grid points Thm (Ditkowski, Fibich, Sagiv, 18): 𝑞 − 𝑞𝑂 1 = 𝑃(𝑂−3/𝒆)

Optimal statistical method (KDE) converges as 𝑂−2/5
Hence, our method is faster for d ≤ 7
For higher dimensions, can use mth-order splines:

Thm (Ditkowski, Fibich, Sagiv, 18): 𝑞 − 𝑞𝑂 1 = 𝑃(𝑂−𝑛/𝒆)

Tensor product spline

48

Lemma 𝑞 𝑧 = 1 𝜈(Ω) න

𝑔−1(𝑧)

1 𝛼𝑔 𝑒𝜏

SLIDE 49

2D example

49

𝛽1 𝛽2

𝑔 𝛽1, 𝛽2 = tanh 6𝛽1𝛽2 + 𝛽1 2 + 𝛽1 + 𝛽2 2 , 𝛽1, 𝛽2 ∼ Uni −1,1 , 𝑗. 𝑗. 𝑒.

𝑔

𝑂 = 82

SLIDE 50

2-dimensional example

50

𝛽1 𝛽2 𝑔

𝑂

𝑂−2.1

𝑔 𝛽1, 𝛽2 = tanh 6𝛽1𝛽2 + 𝛽1 2 + 𝛽1 + 𝛽2 2 , 𝛽1, 𝛽2 ∼ Uni −1,1 , 𝑗. 𝑗. 𝑒. 𝑞 − 𝑞𝑂 1 (𝑢ℎ𝑓𝑝𝑠𝑧: 𝑂−3/2)

𝑂 = 82

SLIDE 51

func = @(x,y) atan(1.3*y.^2 - .4*x.*y+1.1*x.^2)+(x+y) MCpointsPerBlock = 2e3; MCsqrtNumBlock = 10; binsNum = 750; N=10; [xs,ys] = ndgrid(linspace(xmin,xmax,N),linspace(xmin,xmax,N)); [pdf_spline,y_spline] = pdfSampleSquare(@(t,v) interpn(xs,ys,func(xs,ys),t,v,'cubic'),... MCsqrtNumBlock,MCpointsPerBlock,binsNum); function [pdf,binsEdges] = pdfSampleSquare(func,sqrtNumBlocks,sqrtSamplesBlock, numBins) xmin =-1; xmax = 1; smallGridSize = 3e6; blockLength = (xmax-xmin)/sqrtNumBlocks; nsmall =1e3; [x_calibrate,y_calibrate] = ndgrid(linspace(xmin,xmax,1e3),linspace(xmin,xmax,1e3)); [~,binsEdges] = hist(func(x_calibrate,y_calibrate),numBins); binWidth = (max(binsEdges)-min(binsEdges))/numBins; histogram = zeros(1,numBins); for k=1:sqrtNumBlocks for m=1:sqrtNumBlocks [xg,yg] = ndgrid(linspace(xmin+(k-1)*blockLength,xmin+(k)*blockLength, sqrtSamplesBlock),... linspace(xmin+(m-1)*blockLength,xmin+(m)*blockLength,sqrtSamplesBlock)); funcBlock =func(xg,yg); [hist_temp] = hist(funcBlock(:),binsEdges); histogram = histogram+hist_temp; end end pdf = histogram/(sum(histogram)*binWidth);

Matlab code for 2D PDF estimation

51

SLIDE 52

3 dimensional example

52

𝑔 𝛽1, 𝛽2, 𝛽3 = tanh 2𝛽1 + 3𝛽2 + 3𝛽3 + 𝛽1 + 𝛽2 + 𝛽3 3 , 𝛽1, 𝛽2, 𝛽3 ∼ Uni −1,1 , 𝑗. 𝑗. 𝑒.

𝑔

PDF

𝑂

𝑂−1.09

𝑞 − 𝑞𝑂 1

𝑂 = 83

(𝑢ℎ𝑓𝑝𝑠𝑧: 𝑂−1)

SLIDE 53

Conclusions

53

References

A. Sagiv, A. Ditkowski, G. Fibich

A spline-based approach to uncertainty quantification and density estimation ArXiv 1803.10991

A. Sagiv, A. Ditkowski, G. Fibich

Loss of phase and universality of stochastic interactions between laser beams Optics Express 25: 24387-24399, 2017

New method for computing PDF and moments of nl PDEs

with randomness

Outperforms standard statistical methods and gPC
Guaranteed to converge for PDF approximation
Non-intrusive: can use any deterministic numerical solver
Achieves good accuracy using small samples
Extends to multi-dimensional noise
Can also handle non-smooth ``quantity of interest’’