Tweet
Tweet

Slide 1

Slide 1 text

Quasi-Monte Carlo Kernel Density Estimation Fred J. Hickernell, Illinois Tech, hickernell@iit.edu Yuhan Ding, Illinois Tech Brooke Feinberg, Scripps College Guillem Grao I Grasa, Deloitte (formerly Illinois Tech) Aiwen Li, University of Pennsylvania Aadit Jain, Rancho Bernardo High School 
 Larysa Matiukha, Illinois Tech Aleksei Sorokin, Illinois Tech Richard Varela, Sacramento State University August 21, 2024 Thanks to β€’ The organizers β€’ The US National Science Foundation #2316011.& #2244553 See this Jupyter notebook for computations

Slide 2

Slide 2 text

Problem If , where , then we estimate Y = f(X) X ∼ 𝒰 [0,1]d population mean ΞΌ := 𝔼 (Y) = integral ∫ [0,1]d f(x) dx by sample mean 1 n n βˆ‘ i=1 f(xi ) =: Μ‚ ΞΌn But what about approximating the probability density, , of the random variable ? Ο± Y

Slide 3

Slide 3 text

Problem If , where , then we estimate Y = f(X) X ∼ 𝒰 [0,1]d population mean ΞΌ := 𝔼 (Y) = integral ∫ [0,1]d f(x) dx by sample mean 1 n n βˆ‘ i=1 f(xi ) =: Μ‚ ΞΌn But what about approximating the probability density, , of the random variable ? Ο± Y Kernel density estimation (KDE) Ο±(y) = d dy β„™(Y ≀ y) β‰ˆ 1 n n βˆ‘ i=1 k(f(xi ), y) =: KDE(y; {xi }n i=1 , k) ∫ ℝ k(z, y) proportion of z assigned to y dy = 1 ⟹ ∫ ℝ KDE(y; {xi }n i=1 , k) dy = 1 kernel

Slide 4

Slide 4 text

Problem If , where , then we may approximate the probability density, , by a kernel density estimator (KDE): Y = f(X) X ∼ 𝒰 [0,1]d Ο± Ο±(y) = d dy β„™(Y ≀ y) β‰ˆ 1 n n βˆ‘ i=1 k(f(xi ), y) =: Μ‚ Ο±(y; {xi }n i=1 , k), e.g., k(z, y) = ˜ k((z βˆ’ y)/h)/h, ∫ ∞ βˆ’βˆž ˜ k(y) dy = 1

Slide 5

Slide 5 text

Problem If , where , then we may approximate the probability density, , by a kernel density estimator (KDE): Y = f(X) X ∼ 𝒰 [0,1]d Ο± Ο±(y) = d dy β„™(Y ≀ y) β‰ˆ 1 n n βˆ‘ i=1 k(f(xi ), y) =: Μ‚ Ο±(y; {xi }n i=1 , k)

Slide 6

Slide 6 text

Problem If , where , then we may approximate the probability density, , by a kernel density estimator (KDE): Y = f(X) X ∼ 𝒰 [0,1]d Ο± Ο±(y) = d dy β„™(Y ≀ y) β‰ˆ 1 n n βˆ‘ i=1 k(f(xi ), y) =: Μ‚ Ο±(y; {xi }n i=1 , k) β€’ How to choose the kernel, , including bandwidth, ? k h

Slide 7

Slide 7 text

Problem If , where , then we may approximate the probability density, , by a kernel density estimator (KDE): Y = f(X) X ∼ 𝒰 [0,1]d Ο± Ο±(y) = d dy β„™(Y ≀ y) β‰ˆ 1 n n βˆ‘ i=1 k(f(xi ), y) =: Μ‚ Ο±(y; {xi }n i=1 , k) β€’ How to choose the kernel, , including bandwidth, ? k h β€’ Are low discrepancy (LD) points, , better than independent and identically distributed (IID) points? {xi }n i=1

Slide 8

Slide 8 text

Problem If , where , then we may approximate the probability density, , by a kernel density estimator (KDE): Y = f(X) X ∼ 𝒰 [0,1]d Ο± Ο±(y) = d dy β„™(Y ≀ y) β‰ˆ 1 n n βˆ‘ i=1 k(f(xi ), y) =: Μ‚ Ο±(y; {xi }n i=1 , k) β€’ How to choose the kernel, , including bandwidth, ? k h β€’ Are low discrepancy (LD) points, , better than independent and identically distributed (IID) points? {xi }n i=1 β€’ What is the convergence rate for LD points?

Slide 9

Slide 9 text

Problem If , where , then we may approximate the probability density, , by a kernel density estimator (KDE): Y = f(X) X ∼ 𝒰 [0,1]d Ο± Ο±(y) = d dy β„™(Y ≀ y) β‰ˆ 1 n n βˆ‘ i=1 k(f(xi ), y) =: Μ‚ Ο±(y; {xi }n i=1 , k) β€’ How to choose the kernel, , including bandwidth, ? k h β€’ Are low discrepancy (LD) points, , better than independent and identically distributed (IID) points? {xi }n i=1 β€’ What is the convergence rate for LD points? β€’ How to choose to reach the desired error tolerance? n

Slide 10

Slide 10 text

Where does this arise in practice? Y = f(X) = option payo ff underground water pressure with random rock porosity option price average water pressure = ΞΌ

Slide 11

Slide 11 text

Recent work β€’ [BALEOP21] demonstrate that the mean integrated squared error of KDEs using randomized LD sequences is theoretically and empirically superior to using IID sequences. They also discuss how to choose the bandwidth of the kernel. β€’ [LEPBA22] show how conditional Monte Carlo density estimators (CDEs) of using randomized LD sequences has a faster convergence rate than using IID sequences β€’ [LEP22] compare KDEs, CDEs, and likelihood ratio density estimators using randomized LD sequences β€’ The error of interpolated pointwise CDEs of is analyzed in the deterministic setting by [GKS23] Ο± Ο±

Slide 12

Slide 12 text

Why do we think that LD is better?

Slide 13

Slide 13 text

Why do we think that LD is better? Grids do not fi ll space well, IID are better, but LD is best

Slide 14

Slide 14 text

Big picture error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ |Ο±(y) βˆ’ ˜ Ο±(y; k)| + |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ˜ Ο±(y; k):= 𝔼 [k(Y, y)] = ∫ ∞ βˆ’βˆž k(z, y) Ο±(z) dz

Slide 15

Slide 15 text

Big picture error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ |Ο±(y) βˆ’ ˜ Ο±(y; k)| + |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ˜ Ο±(y; k):= 𝔼 [k(Y, y)] = ∫ ∞ βˆ’βˆž k(z, y) Ο±(z) dz If and k(z, y) = ˜ k((z βˆ’ y)/h)/h |Ο±(y) βˆ’ ˜ Ο±(y; k)| = π’ͺ (hp) and |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| = π’ͺ (nβˆ’qhβˆ’r) Then an optimal choice of is h hopt = π’ͺ (nβˆ’q/(p+r)) and |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| = π’ͺ (nβˆ’q/(1+r/p))

Slide 16

Slide 16 text

Big picture error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ |Ο±(y) βˆ’ ˜ Ο±(y; k)| + |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ˜ Ο±(y; k):= 𝔼 [k(Y, y)] = ∫ ∞ βˆ’βˆž k(z, y) Ο±(z) dz If and k(z, y) = ˜ k((z βˆ’ y)/h)/h |Ο±(y) βˆ’ ˜ Ο±(y; k)| = π’ͺ (hp) and |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| = π’ͺ (nβˆ’qhβˆ’r) Then an optimal choice of is h hopt = π’ͺ (nβˆ’q/(p+r)) and |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| = π’ͺ (nβˆ’q/(1+r/p)) For IID, q = r = 1/2 hIID opt = π’ͺ (nβˆ’1/(2p+1)) and |Ο±(y) βˆ’ KDE(y; {xIID i }n i=1 , k)| = π’ͺ (nβˆ’1/(2+1/p))

Slide 17

Slide 17 text

Deterministic error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ |Ο±(y) βˆ’ ˜ Ο±(y; k)| + |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ˜ Ο±(y; k):= 𝔼 [k(Y, y)] = ∫ ∞ βˆ’βˆž k(z, y) Ο±(z) dz

Slide 18

Slide 18 text

Deterministic error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ |Ο±(y) βˆ’ ˜ Ο±(y; k)| + |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ˜ Ο±(y; k):= 𝔼 [k(Y, y)] = ∫ ∞ βˆ’βˆž k(z, y) Ο±(z) dz |Ο±(y) βˆ’ ˜ Ο±(y; k)| ≀ smooth(k, Ky ) density independent βˆ₯Ο±βˆ₯Ky ⏟ kernel independent smooth2(k, Ky ) = Ky (y, y) βˆ’ 2 ∫ ∞ βˆ’βˆž k(z, y) Ky (z, y) dz + ∫ ∞ βˆ’βˆž ∫ ∞ βˆ’βˆž k(z, y) Ky (z, t) k(t, y) dz dt

Slide 19

Slide 19 text

Deterministic error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ |Ο±(y) βˆ’ ˜ Ο±(y; k)| + |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ˜ Ο±(y; k):= 𝔼 [k(Y, y)] = ∫ ∞ βˆ’βˆž k(z, y) Ο±(z) dz |Ο±(y) βˆ’ ˜ Ο±(y; k)| ≀ smooth(k, Ky ) density independent βˆ₯Ο±βˆ₯Ky ⏟ kernel independent smooth2(k, Ky ) = Ky (y, y) βˆ’ 2 ∫ ∞ βˆ’βˆž k(z, y) Ky (z, y) dz + ∫ ∞ βˆ’βˆž ∫ ∞ βˆ’βˆž k(z, y) Ky (z, t) k(t, y) dz dt |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ discrepancy({x}n i=1 ) kernel independent βˆ₯k(f( β‹… ), y)βˆ₯Kx sample independent discrepancy2({x}n i=1 ) = ∫ [0,1]dΓ—[0,1]d Kx (x, t) dx dt βˆ’ 2 n n βˆ‘ i=1 ∫ [0,1]d Kx (x, xi ) dx + 1 n2 n βˆ‘ i,j=1 Kx (xi , xj )

Slide 20

Slide 20 text

Deterministic error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ |Ο±(y) βˆ’ ˜ Ο±(y; k)| + |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ˜ Ο±(y; k):= 𝔼 [k(Y, y)] = ∫ ∞ βˆ’βˆž k(z, y) Ο±(z) dz |Ο±(y) βˆ’ ˜ Ο±(y; k)| ≀ smooth(k, Ky ) density independent βˆ₯Ο±βˆ₯Ky ⏟ kernel independent smooth2(k, Ky ) = Ky (y, y) βˆ’ 2 ∫ ∞ βˆ’βˆž k(z, y) Ky (z, y) dz + ∫ ∞ βˆ’βˆž ∫ ∞ βˆ’βˆž k(z, y) Ky (z, t) k(t, y) dz dt |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ discrepancy({x}n i=1 ) kernel independent βˆ₯k(f( β‹… ), y)βˆ₯Kx sample independent discrepancy2({x}n i=1 ) = ∫ [0,1]dΓ—[0,1]d Kx (x, t) dx dt βˆ’ 2 n n βˆ‘ i=1 ∫ [0,1]d Kx (x, xi ) dx + 1 n2 n βˆ‘ i,j=1 Kx (xi , xj ) and are reproducing kernels Ky Kx

Slide 21

Slide 21 text

Deterministic error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ |Ο±(y) βˆ’ ˜ Ο±(y; k)| + |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ˜ Ο±(y; k):= 𝔼 [k(Y, y)] = ∫ ∞ βˆ’βˆž k(z, y) Ο±(z) dz |Ο±(y) βˆ’ ˜ Ο±(y; k)| ≀ smooth(k, Ky ) density independent βˆ₯Ο±βˆ₯Ky ⏟ kernel independent smooth2(k, Ky ) = Ky (y, y) βˆ’ 2 ∫ ∞ βˆ’βˆž k(z, y) Ky (z, y) dz + ∫ ∞ βˆ’βˆž ∫ ∞ βˆ’βˆž k(z, y) Ky (z, t) k(t, y) dz dt |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ discrepancy({x}n i=1 ) kernel independent βˆ₯k(f( β‹… ), y)βˆ₯Kx sample independent discrepancy2({x}n i=1 ) = ∫ [0,1]dΓ—[0,1]d Kx (x, t) dx dt βˆ’ 2 n n βˆ‘ i=1 ∫ [0,1]d Kx (x, xi ) dx + 1 n2 n βˆ‘ i,j=1 Kx (xi , xj ) and are reproducing kernels Ky Kx smooth(k, Ky ) ↓ ⟺ βˆ₯k(f( β‹… ), y)βˆ₯Kx ↑

Slide 22

Slide 22 text

Deterministic error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ smooth(k, Ky ) βˆ₯Ο±βˆ₯Ky + discrepancy({xi }n i=1 ) π’ͺ (nβˆ’1+Ξ΄) βˆ₯K(f( β‹… ), y)βˆ₯Kx

Slide 23

Slide 23 text

Deterministic error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ smooth(k, Ky ) βˆ₯Ο±βˆ₯Ky + discrepancy({xi }n i=1 ) π’ͺ (nβˆ’1+Ξ΄) βˆ₯K(f( β‹… ), y)βˆ₯Kx For k(z, y) = ˜ k((z βˆ’ y)/h)/h), ˜ k even smooth(k, Ky ) = Ky (y, y) βˆ’ 2 ∫ ∞ βˆ’βˆž k(z, y) Ky (z, y) dz + ∫ ∞ βˆ’βˆž ∫ ∞ βˆ’βˆž k(z, y) Ky (z, t) k(t, y) dz dt = π’ͺ (hr+2) if ∫ ∞ βˆ’βˆž ˜ k(w) w2l dw = Ξ΄0,l , l = 0,…, r

Slide 24

Slide 24 text

Deterministic error analysis |Ο±(y) βˆ’ KDE(y; {xi }n i=1 , k)| ≀ smooth(k, Ky ) βˆ₯Ο±βˆ₯Ky + discrepancy({xi }n i=1 ) π’ͺ (nβˆ’1+Ξ΄) βˆ₯K(f( β‹… ), y)βˆ₯Kx r = 0 r = 1 For k(z, y) = ˜ k((z βˆ’ y)/h)/h), ˜ k even smooth(k, Ky ) = π’ͺ (hr+2) if ∫ ∞ βˆ’βˆž ˜ k(w) w2l dw = Ξ΄0,l , l = 0,…, r

Slide 25

Slide 25 text

Deterministic error analysis |˜ Ο±(y; k) βˆ’ KDE(y; {xi }n i=1 , k)| = ∫ ∞ βˆ’βˆž g(x; y) dx βˆ’ 1 n n βˆ‘ i=1 g(xi ; y) , g(x; y) := k(f(x), y) ≀ discrepancy({xi }n i=1 ) π’ͺ (nβˆ’1+Ξ΄) βˆ₯g( β‹… ; y)βˆ₯Kx E.g., discrepancy2({xi }n i=1 ) = d ∏ β„“=1 ( 1 + Ξ³2 β„“ 12) βˆ’ 2 n n βˆ‘ i=1 d ∏ β„“=1 [ 1 + Ξ³2 β„“ 2 ( xiβ„“ βˆ’ 1/2 βˆ’ xiβ„“ βˆ’ 1/2 2 )] + 1 n2 n βˆ‘ i,j=1 d ∏ β„“=1 [ 1 + Ξ³2 β„“ 2 ( xiβ„“ βˆ’ 1/2 + xjβ„“ βˆ’ 1/2 βˆ’ xiβ„“ βˆ’ xjβ„“ )] variation2(g) = ∫ [0,1] βˆ‚g(x1 , 1/2) Ξ³1 βˆ‚x1 2 dx1 + β‹― + ∫ [0,1]2 βˆ‚2g(x1 , x2 , 1/2) Ξ³1 Ξ³2 βˆ‚x1 βˆ‚x2 2 dx1 dx2 +β‹― + ∫ [0,1]d βˆ‚dg(x) Ξ³1 β‹―Ξ³d βˆ‚x1 β‹―βˆ‚xd 2 dx

Slide 26

Slide 26 text

Deterministic error analysis g(x; y):= k(f(x), y) = ˜ k((f(x) βˆ’ y)/h)/h E.g., variation2(g) = ∫ [0,1] βˆ‚g(x1 , 1/2) Ξ³1 βˆ‚x1 2 dx1 + β‹― + ∫ [0,1]2 βˆ‚2g(x1 , x2 , 1/2) Ξ³1 Ξ³2 βˆ‚x1 βˆ‚x2 2 dx1 dx2 +β‹― + ∫ [0,1]d βˆ‚dg(x) Ξ³1 β‹―Ξ³d βˆ‚x1 β‹―βˆ‚xd 2 dx βˆ‚g(x 𝔲 , 1/2) βˆ‚x 𝔲 = βˆ‘ P∈Π( 𝔲 ) 1 h|P|+1 ˜ k(|P|) ( f(x 𝔲 , 1/2) βˆ’ y h )∏ 𝔳 ∈P βˆ‚| 𝔳 | (f(x 𝔳 , 1/2)) βˆ‚x 𝔳 Ξ ( 𝔲 ) := set of partitions of 𝔲 Even if is additive, i.e., , f f(x) = f1 (x1 ) + β‹― + fd (xd ) variation(˜ k(f(x βˆ’ y)/h)/h) = π’ͺ (hβˆ’d)?

Slide 27

Slide 27 text

Ex. sum of uniforms, squared exponential K Y = X1 + β‹― + Xd , X ∼ 𝒰 [0,1]d

Slide 28

Slide 28 text

Ex. sum of uniforms, squared exponential K Y = X1 + β‹― + Xd , X ∼ 𝒰 [0,1]d

Slide 29

Slide 29 text

Ex. sum of uniforms, squared exponential K Y = X1 + β‹― + Xd , X ∼ 𝒰 [0,1]d errorIID = π’ͺ (nβˆ’2/5) errorLD = π’ͺ (nβˆ’1/(1+d/2))

Slide 30

Slide 30 text

Ex. sum of uniforms, squared exponential K Y = X1 + β‹― + Xd , X ∼ 𝒰 [0,1]d hIID = π’ͺ (nβˆ’1/5) hLD = π’ͺ (nβˆ’1/(2+d))

Slide 31

Slide 31 text

Ex. weighted sum of Gaussians, squared exponential K Y = Ξ³1 X1 + β‹― + Ξ³d Xd , Xβ„“ IID ∼ 𝒩 (0,1), Ξ³β„“ = 1/β„“

Slide 32

Slide 32 text

Ex. weighted sum of Gaussians, squared exponential K Y = Ξ³1 X1 + β‹― + Ξ³d Xd , Xβ„“ IID ∼ 𝒩 (0,1), Ξ³β„“ = 1/β„“ errorIID = π’ͺ (nβˆ’2/5) errorLD = π’ͺ (nβˆ’1/(1+d/2))

Slide 33

Slide 33 text

Ex. weighted sum of Gaussians, Gauss kernel Y = Ξ³1 X1 + β‹― + Ξ³d Xd , Xβ„“ IID ∼ 𝒩 (0,1), Ξ³β„“ = 1/β„“ hIID = π’ͺ (nβˆ’1/5) hLD = π’ͺ (nβˆ’1/(2+d))

Slide 34

Slide 34 text

Ex. sum of uniforms, squared exp x quad K Y = X1 + β‹― + Xd , X ∼ 𝒰 [0,1]d

Slide 35

Slide 35 text

Ex. sum of uniforms, squared exp x quad K Y = X1 + β‹― + Xd , X ∼ 𝒰 [0,1]d errorIID = π’ͺ (nβˆ’3/7) errorLD = π’ͺ (nβˆ’1/(1+d/3))

Slide 36

Slide 36 text

Ex. sum of uniforms, squared exponential K Y = X1 + β‹― + Xd , X ∼ 𝒰 [0,1]d errorIID = π’ͺ (nβˆ’2/5) errorLD = π’ͺ (nβˆ’1/(1+d/2))

Slide 37

Slide 37 text

Ex. weighted sum of Gaussians, squared exp x quad K Y = Ξ³1 X1 + β‹― + Ξ³d Xd , Xβ„“ IID ∼ 𝒩 (0,1), Ξ³β„“ = 1/β„“ errorIID = π’ͺ (nβˆ’3/7) errorLD = π’ͺ (nβˆ’1/(1+d/3))

Slide 38

Slide 38 text

Ex. weighted sum of Gaussians, squared exponential K Y = Ξ³1 X1 + β‹― + Ξ³d Xd , Xβ„“ IID ∼ 𝒩 (0,1), Ξ³β„“ = 1/β„“ errorIID = π’ͺ (nβˆ’2/5) errorLD = π’ͺ (nβˆ’1/(1+d/2))

Slide 39

Slide 39 text

Conditional Monte Carlo for Density Estimation (CDE) Estimate the probability density, , for , where Ο± Y = f(X) X ∼ 𝒰 [0,1]d

Slide 40

Slide 40 text

Conditional Monte Carlo for Density Estimation (CDE) Estimate the probability density, , for , where Ο± Y = f(X) X ∼ 𝒰 [0,1]d y = f(x) ⟺ x1 = u(y; x2:d ) If one can identify , such that u

Slide 41

Slide 41 text

Conditional Monte Carlo for Density Estimation (CDE) Estimate the probability density, , for , where Ο± Y = f(X) X ∼ 𝒰 [0,1]d y = f(x) ⟺ x1 = u(y; x2:d ) Then Ο±(y) = ∫ [0,1]dβˆ’1 uβ€² οΏΌ (y; x) dx If one can identify , such that u

Slide 42

Slide 42 text

Ex. sum of uniforms, KDE vs. CDE Y = X1 + β‹― + Xd , X ∼ 𝒰 [0,1]d

Slide 43

Slide 43 text

Ex. weighted sum of Gaussians, KDE vs CDE Y = Ξ³1 X1 + β‹― + Ξ³d Xd , Xβ„“ IID ∼ 𝒩 (0,1), Ξ³β„“ = 1/β„“

Slide 44

Slide 44 text

(Preliminary) conclusions β€’ KDE allows density estimation for black box β€’ Theory for LD sequences for KDE matches practice sometimes β€’ Performance of LD sequences for KDE compared to IID is disappointing for larger ; can this be overcome? β€’ Don’t know how to β€’ Adjust to accommodate known boundaries in the sample space of β€’ Choose bandwidth from data β€’ Choose to get desired accuracy β€’ CDE, when available, outperforms KDE f d k Y n

Slide 45

Slide 45 text

References [BALEOP21] A. Ben Abdellah, P. L'Ecuyer, A. B. Owen, and F. Puchhammer, Density estimation by randomized quasi-Monte Carlo, SIAM/ASA J. Uncertain. Quantif. 9 (2021), pp. 280-301. [GKS23] A. D. Gilbert, F. Y. Kuo, and I. H. Sloan, Analysis of preintegration followed by quasi-Monte Carlo integration for distribution functions and densities, SIAM J. Numer. Anal. 61 (2023), 135–166. [LEP22] P. L'Ecuyer and F. Puchhammer, Density estimation by Monte Carlo and quasi-Monte Carlo, Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Oxford, England, August 2020 (A. Keller, ed.), Springer Proceedings in Mathematics and Statistics, Springer, Cham, 2022, pp. 3–21. LEPBA22] P. L'Ecuyer, F. Puchhammer, and A. Ben Abdellah, Monte Carlo and quasi-Monte Carlo density estimation via conditioning, INFORMS J. Comput. 34 (2022), no. 3, 1729–1748.