Piecewise polynomial approximation of probability density functions with application to uncertainty quantification for stochastic PDEs

TL;DR

This paper introduces a scalable, linear polynomial-based method for approximating probability density functions from samples, applicable to uncertainty quantification in stochastic PDEs, avoiding complex computations and bandwidth tuning.

Contribution

It presents a novel piecewise polynomial density estimator that is scalable, does not require solving linear systems, and simplifies bandwidth selection, with applications to stochastic PDEs.

Findings

The estimator is consistent and enforces unit integral.

It does not require positivity constraints.

Validated through numerical examples with known and unknown PDFs.

Abstract

The probability density function (PDF) associated with a given set of samples is approximated by a piecewise-linear polynomial constructed with respect to a binning of the sample space. The kernel functions are a compactly supported basis for the space of such polynomials, i.e. finite element hat functions, that are centered at the bin nodes rather than at the samples, as is the case for the standard kernel density estimation approach. This feature naturally provides an approximation that is scalable with respect to the sample size. On the other hand, unlike other strategies that use a finite element approach, the proposed approximation does not require the solution of a linear system. In addition, a simple rule that relates the bin size to the sample size eliminates the need for bandwidth selection procedures. The proposed density estimator has unitary integral, does not require a constraint to enforce positivity, and is consistent. The proposed approach is validated through numerical examples in which samples are drawn from known PDFs. The approach is also used to determine approximations of (unknown) PDFs associated with outputs of interest that depend on the solution of a stochastic partial differential equation.

Tables

Table 1. Table 1: Computational time (in seconds) for determining the approximation ( 5 ) of the truncated standard Gaussian PDF ( 17 ).

				Computation time in seconds
$k$	$N_{\delta }$	$\delta =\frac{11}{N_{\delta }}$	$M=2^{4k}$	$N_{bins}$	$N_{\Gamma }=1$	$N_{bins}$	$N_{\Gamma }=2$	$N_{bins}$	$N_{\Gamma }=3$
2	4	2.75	256	4	$7.710\text{e}-04$	16	$3.777\text{e}-02$	64	$1.190\text{e}-01$
3	8	1.375	4096	8	$1.356\text{e}-02$	64	$5.114\text{e}-01$	512	$1.955\text{e}-00$
4	16	0.6875	65536	16	$2.796\text{e}-01$	256	$8.106\text{e}-00$	4096	$3.545\text{e+01}$
5	32	0.34375	1048576	32	$5.870\text{e}-00$	1024	$1.281\text{e+02}$	32768	$1.013\text{e+03}$

Equations

$$\displaystyle{f}^{histo}_{\delta,M}(Y)\big{|}_{\mathcal{B}_{\ell}}=\dfrac{1}{\delta M}\sum_{m=1}^{M}{\mathcal{X}}_{\mathcal{B}_{\ell}}(Y_{m}),\quad\ell=1,\ldots,N_{bins},$$

$$\displaystyle{f}^{kde}(Y)=\dfrac{1}{bM}\sum_{m=1}^{M}{\mathcal{K}}\Big{(}\dfrac{Y-Y_{m}}{b}\Big{)},$$

$${\mathcal{V}}_{\delta}=\big{\{}v\in C(\Gamma)\,\,:\,\,v|_{\mathcal{B}_{\ell}}\in{\mathcal{P}}_{1}(\mathcal{B}_{\ell})\quad\mbox{for $\ell=1,\ldots,N_{bins}$}\big{\}},$$

$$\phi_{j}({\bm{Y}})=\big{\{}\phi_{j}({\bm{Y}})\in{\mathcal{V}}_{\delta}\,:\,\phi_{j}({\widehat{\bm{Y}}}_{j^{\prime}})=\delta_{jj^{\prime}}\quad\mbox{for $j^{\prime}=1,\ldots,N_{nodes}$}\big{\}},$$

$$\sum _{j=1}^{N_{nodes}}{\varphi }_j(\mathbf{Y})=1\forall \mathbf{Y}\in \Gamma .$$

$$C_j={\int }_{\Gamma }{\varphi }_j(\mathbf{Y})d\mathbf{Y}={\int }_{S_j(\mathbf{Y})}{\varphi }_j(\mathbf{Y})d\mathbf{Y}=\sum _{{\mathcal{B}}_{\ell }\in S_j(\mathbf{Y})}{\int }_{{\mathcal{B}}_{\ell }}{\varphi }_j(\mathbf{Y})d\mathbf{Y}.$$

$$\framebox{$\displaystyle\begin{aligned} &f({\bm{Y}})\approx f_{\delta,M}({\bm{Y}})=\sum_{j=1}^{N_{nodes}}F_{j}\phi_{j}({\bm{Y}})\in{\mathcal{V}}_{\delta},\\[4.30554pt] &\qquad\qquad\mbox{where}\quad F_{j}=\dfrac{1}{MC_{j}}\,\sum_{{\bm{Y}}_{m}\in S_{j}({\bm{Y}})}\phi_{j}({\bm{Y}}_{m})\\[4.30554pt] &\qquad\qquad\qquad\mbox{with}\quad C_{j}=\int_{S_{j}({\bm{Y}})}\phi_{j}({\bm{Y}})d{\bm{Y}},\,\,j=1,\ldots,N_{nodes}.\end{aligned}$}$$

$${\int }_{\Gamma }{f}_{\delta ,M}(\mathbf{Y})d\mathbf{Y}$$

$$=\sum _{j=1}^{N_{nodes}}\frac{1}{M{C}_j}\sum _{m=1}^M{\varphi }_j({\mathbf{Y}}_m){\int }_{\Gamma }{\varphi }_j(\mathbf{Y})d\mathbf{Y}$$

$$=\frac{1}{M}\sum _{i=1}^{N_{nodes}}\sum _{m=1}^M{\varphi }_j({\mathbf{Y}}_m)=\frac{1}{M}\sum _{m=1}^{M}1=1.$$

$$\displaystyle\Big{|}f(\widehat{{\bm{Y}}}_{j})-\mathbb{E}[F_{j}]\Big{|}\leq C\delta^{\alpha},$$

$$\displaystyle\phi_{j}(\mathbf{Y}_{m})=\phi\Big{(}\dfrac{\widehat{{\bm{Y}}}_{j}-\mathbf{Y}_{m}}{\delta}\Big{)},\qquad F_{j}=\dfrac{1}{\delta^{N_{\Gamma}}\,M}\sum_{m=1}^{M}\phi\Big{(}\dfrac{\widehat{{\bm{Y}}}_{j}-\mathbf{Y}_{m}}{\delta}\Big{)}.$$

$$\displaystyle\Big{|}f(\widehat{{\bm{Y}}}_{j})-\mathbb{E}[F_{j}]\Big{|}\leq$$

$$\displaystyle\mbox{+}\dfrac{\delta^{2}}{2}\Big{|}\int_{[-1,1]^{N_{\Gamma}}\cap\Gamma}{\mathbf{Y}^{{}^{\prime}}}^{T}\dfrac{\partial^{2}f(\widehat{{\bm{Y}}}_{j})}{\partial\mathbf{Y}^{2}}\mathbf{Y}^{{}^{\prime}}d\mathbf{Y}^{{}^{\prime}}\Big{|}+\mathcal{O}(\delta^{2}).$$

$$\underset{\delta \to 0}{\lim }\underset{M\to \infty }{\lim }\Vert f-f_{\delta ,M}{\Vert }_{L^2(\Gamma )}=0.$$

$$\underset{M\to \infty }{\lim }\Vert f-f_{\delta ,M}{\Vert }_{L^2(\Gamma )}\le C{\delta }^2,$$

$$\Vert f-f_{\delta ,M}{\Vert }_{L^2(\Gamma )}$$

$$\le C_1{\delta }^2\text{mbox}+\Vert {\mathcal{I}}_{\delta }f-f_{\delta ,M}{\Vert }_{L^2(\Gamma )},$$

$$\Vert {\mathcal{I}}_{\delta }f-f_{\delta ,M}{\Vert }_{L^2(\Gamma )}$$

$$\displaystyle\leq\sqrt{\int_{\Gamma}\Big{[}\Big{(}\sum\limits_{j=1}^{N_{nodes}}|f(\widehat{{\bm{Y}}}_{j})-\mathbb{E}[F_{j}]|\phi_{j}\Big{)}\mbox{+}\Big{(}\sum\limits_{j=1}^{N_{nodes}}|\mathbb{E}[F_{j}]-F_{j}|\phi_{j}\Big{)}\Big{]}^{2}}$$

$$\displaystyle\leq\sqrt{\int_{\Gamma}\Big{[}C_{2}\delta^{\alpha}\mbox{+}\Big{(}\sum\limits_{j=1}^{N_{nodes}}|\mathbb{E}[F_{j}]-F_{j}|\phi_{j}\Big{)}\Big{]}^{2}}.$$

$$|\mathbb{E}[F_j]-F_j|$$

$$\displaystyle=\dfrac{1}{\delta^{N_{\Gamma}}}\Big{|}\mathbb{E}\Big{[}\phi_{j}\Big{]}-\dfrac{1}{M}\sum\limits_{m=1}^{M}\phi_{j}(\mathbf{Y}_{m})\Big{|}\leq\dfrac{\sigma(\phi_{j})}{\delta^{N_{\Gamma}}\,\sqrt{M}},$$

$$\Vert f-f_{\delta ,M}{\Vert }_{L^2(\Gamma )}\le C_1{\delta }^2\text{mbox}+C_2{\delta }^{\alpha }\text{mbox}+\frac{C_3}{\sqrt{{\delta }^{N_{\Gamma }}}\sqrt{M}},$$

$${\mathcal{E}}_{f^{approx}}=\Big{(}\dfrac{1}{M}\sum_{m=1}^{M}\big{(}f({\bm{Y}}_{m})-{f}^{approx}({\bm{Y}}_{m})\big{)}^{2}\Big{)}^{1/2}.$$

$$M=(b-a{)}^{2r}{\delta }^{-2r}={N_{\delta }}^{2r}.$$

$$N_{\delta }=2^{(3-r)k}$$

$$\displaystyle f({\bm{Y}})=\prod\limits_{n=1}^{N_{\Gamma}}\dfrac{1}{\sqrt{2\pi}C_{G}}\exp{\Big{(}-\dfrac{Y_{n}^{2}}{2}\Big{)}}\qquad\mbox{for ${\bm{Y}}\in\Gamma=[-5.5,5.5]^{N_{\Gamma}}$}$$

$$\displaystyle\qquad\qquad\mbox{with}\quad C_{G}=\frac{1}{2}\big{(}\mbox{erf}(5.5/\sqrt{2})-\mbox{erf}(-5.5/\sqrt{2})\big{)}.$$

$$M=10^7\text{mbox}and\delta =11/2^k\text{mbox}fork=3,4,5,6,$$

$$\delta =11/2^8\text{mbox}andM=10^k\text{mbox}fork=3,4,5,6,$$

$$f_{[-1.5,1.5]}(Y)=\left\{\begin{aligned} 1&\quad\mbox{for $Y\in[-1,1]$}\\ 0&\quad\mbox{for $Y\in[-1.5,-1]$ and $Y\in[1,1.5]$.}\end{aligned}\right.$$

$$f(Y)=\dfrac{1}{3C_{L}}\exp{\Big{(}\dfrac{-|Y|}{1.5}\Big{)}}$$

$$f({\bm{Y}})=\dfrac{1}{\sqrt{8\pi}C^{\prime}_{G}}\exp{\Big{(}-\frac{Y_{1}^{2}}{8}\Big{)}}\dfrac{1}{\sqrt{2\pi}C_{G}}\exp{\Big{(}-\frac{Y_{2}^{2}}{2}\Big{)}},$$

Full text

11footnotetext: Funding: This work was supported by US Air Force Office of Scientific Research grant FA9550-15-1-0001 and by the Sandia National Laboratories contract 1985151.

1 Introduction

The problem of estimating a probability density function (PDF) associated with a given set of samples is of major relevance in a variety of mathematical and statistical applications; see, e.g., [1, 5, 6, 9, 11, 14, 20, 27, 31]. Histograms are perhaps the most popular means used in practice for this purpose. A histogram is a piecewise-constant approximation of an unknown PDF that is based on the subdivision of the sample space into subdomains that are commonly referred as bins. For simplicity, consider the case in which all bins have equal volume. The value of the histogram on any bin is given by the number of samples that lie within that bin along with a global scaling applied to ensure that the histogram has unitary integral. As an example, assume that a bounded one-dimensional sample domain $\Gamma$ has been discretized into a set $\mathcal{B}^{\delta}=\{\mathcal{B}_{\ell}\}_{\ell=1}^{N_{bins}}$ of $N_{bins}$ non-overlapping, covering, bins (i.e., intervals) of length $\delta$. Assume also that one has in hand $M$ samples $\{Y_{m}\}_{m=1}^{M}$ of an unknown PDF $f(Y)$. For every $\mathcal{B}_{\ell}\in{\mathcal{B}^{\delta}}$, the value of the histogram function ${f}^{histo}_{\delta,M}(Y)$ on $\mathcal{B}_{\ell}$ is given by

[TABLE]

where ${\mathcal{X}}_{\mathcal{B}_{\ell}}(Y)$ denotes the indicator function for $\mathcal{B}_{\ell}$. The above formula easily generalizes to higher-dimensional sample domains. Although the histogram is probably the easiest, with regards to implementation, means for estimating a PDF, it suffers from some limitations. For instance, being only a piecewise-constant approximation, it is discontinuous across bin boundaries and also the use of piecewise-constant approximations severely limits the discretization accuracy one can achieve.

A popular alternative method capable of overcoming the differentiability issue is kernel density estimation (KDE) for which the PDF is approximated by a sum of kernel functions ${\mathcal{K}}(\cdot)$ centered at the samples so that a desired smoothness of the approximation can be obtained [20]. Considering again a one-dimensional sample domain $\Gamma$, the KDE approximation is defined as

[TABLE]

where $b$, which is referred to as the bandwidth, usually governs the decay rate of the kernel as $|Y-Y_{m}|$ increases. KDE approximations have been shown to be effective in a variety of applications; see, e.g., [14, 22, 29, 30]. A limitation of KDE is that the choice of the bandwidth strongly affects the accuracy [19, 28] of the approximation. More important, the naive KDE method of (2) does not scale well with the dimension $M$ of the sample data set, i.e., for a given $Y$, the evaluation of the KDE approximation (2) requires $M$ kernel evaluations so that clearly evaluating (2) becomes more expensive as the value of $M$ grows. A way to overcome this issue is performing an appropriate binning of the sample set, as for the histogram, and appropriately transferring the information from the samples to the grid points. The kernel functions are then centered at the grid points and scalability with respect to the sample size can be achieved [12, 17]. Another method related to KDE is the one presented in [18], which is based on spline smoothing and on a finite element discretization of the estimator. This method is related to our approach because the PDF is also approximated by a finite element function. Our method presents several advantages compared to that of [18]. First, to determine the coefficients, the solution of a linear system is not required. Second, our method only involves the binning size as a smoothing parameter, whereas the one in [18] also requires the treatment of an additional smoothing parameter $\lambda$. Third, no special constraints have to be introduced to ensure positivity of the PDF approximation. In turn, the method in [18] has features that our approach does not provide such as the ability to match the sample moments up to a certain degree and the possibility of allowing aggregate data. We note that in [25], sparse-grid basis functions are substituted for the standard finite element basis functions in the method of [18], allowing for consideration of larger $N_{\Gamma}$.

In our approach, the unknown PDF is approximated by a piecewise-polynomial function, specifically a piecewise-linear polynomial, that is defined, as are histograms, with respect to a subdivision of the sample domain $\Gamma$ into bins. The piecewise-linear approximation we propose is a finite element function obtained as a linear combination of hat functions. The procedure to determine the coefficients of the linear combination is purely algebraic and can be efficiently carried out. The approach is consistent in the sense that the approximate PDF converges to the exact PDF with respect to the $L^{2}(\Gamma)$ norm as the number of samples tends to infinity and the volume of the bins tend to zero. Moreover, the only smoothing parameter involved is the bin size that can be related heuristically to the sample size by a simple rule.

The paper is structured as follows. In Section 2, the mathematical foundation of our approach is laid out; there, it is shown analytically that the proposed approximation satisfies several requirements needed for it to be considered as a PDF estimator. A numerical investigation of the accuracy and computational costs incurred by our approach is then provided. First, in Section 3, our method is tested and validated using sample sets associated with different types of known PDFs so that exact errors can be determined. Then, in Section 4, the method is applied to the estimation of unknown PDFs of outputs of interest associated with the solution of a stochastic partial differential equation. Finally, concluding remarks and future work are discussed in Section 5.

2 Piecewise-linear polynomial approximations of PDFs

Let ${\bm{Y}}$ denote a multivariate random variable belonging to a closed, bounded parameter domain $\Gamma\subset\mathbb{R}^{N_{\Gamma}}$ which, for simplicity, we assume is a polytope in $\mathbb{R}^{N_{\Gamma}}$. The probability density function (PDF) $f({\bm{Y}})$ corresponding to ${\bm{Y}}$ is not known. However, we assume that we have in hand a data set of $M$ samples ${\bm{Y}}_{m}\in\Gamma$, $m=1,\ldots,M$, of ${\bm{Y}}\in\Gamma$. The goal is to construct an estimator for $f({\bm{Y}})$ using the given data set $\{{\bm{Y}}_{m}\}_{m=1}^{M}$ of samples. To this end, we use an approximating space that is popular in the finite element community [7, 10].

Let ${\mathcal{B}}^{\delta}=\{\mathcal{B}_{\ell}\}_{\ell=1}^{N_{bins}}$ denote a covering, non-overlapping subdivision of the sample domain $\Gamma$ into ${N_{bins}}$ bins.222In the partial differential equation (PDE) setting, what we refer to as bins are often referred to as grid cells or finite elements or finite volumes. We instead refer to the subdomains $\{\mathcal{B}_{\ell}\}_{\ell=1}^{N_{bins}}$ as bins because that is the notation in common use for histograms which we use to compare to our approach. Furthermore, in Section 4, we also use finite element grids for spatial discretization of partial differential equations, so that using the notation “bins” for parameter domain subdivisions helps us differentiate between subdivisions of parameter and spatial domains. For the same reason, we use $\delta$ instead of $h$ to parametrize parameter bin sizes because $h$ is in common use to parametrize spatial grid sizes. Here, $\delta$ parametrizes the subdivision and may be taken as, e.g, the largest diameter of any of the bins $\{\mathcal{B}_{\ell}\}$. The bins are chosen to be hyper-quadrilaterals; for example, if $N_{\Gamma}=2$, they would be quadrilaterals. It is also assumed that the faces of the bins are either also complete faces of abutting bins or part of the boundary of $\Gamma$. From a practical point of view, our considerations are limited to relatively small ${N_{\Gamma}}$ because ${N_{bins}}={\mathcal{O}}(1/\delta^{N_{\Gamma}})$. Detailed discussion about the subdivisions we use can be found in, e.g., [7, 10].

Let $\{\widehat{{\bm{Y}}}_{j}\}_{j=1}^{{N_{nodes}}}$ denote the set of nodes, i.e., vertices, of ${\mathcal{B}}^{\delta}$ with ${N_{nodes}}$ denoting the number of nodes of ${\mathcal{B}}^{\delta}$. Note that we also have that ${N_{nodes}}={\mathcal{O}}(1/\delta^{N_{\Gamma}})$. Based on the subdivision ${\mathcal{B}}^{\delta}$, we define the space of continuous piecewise polynomials

[TABLE]

where, for hyper-quadrilateral elements, ${\mathcal{P}}_{1}(\cdot)$ denotes the space of $N_{\Gamma}$-linear polynomials, e.g., bilinear and trilinear polynomials in two and three dimensions, respectively.

A basis $\{\phi_{j}({\bm{Y}})\}_{j=1}^{N_{nodes}}$ for ${\mathcal{V}}_{\delta}$ is given by, for $j=1,\ldots,N_{nodes}$,

[TABLE]

where $\delta_{jj^{\prime}}$ denotes the Kronecker delta function. In detail, we have that $\{\phi_{j}({\bm{Y}})\}_{j=1}^{N_{nodes}}$ denotes the continuous piecewise-linear or piecewise ${N_{\Gamma}}$-linear Lagrangian FEM basis corresponding to ${\mathcal{B}}^{\delta}$, i.e, we have that, for $j=1,\ldots,{N_{nodes}}$,

• –

for hyper-quadrilateral bins, $\phi_{j}({\bm{Y}})$ is an ${N_{\Gamma}}$-linear function on each bin $\mathcal{B}_{\ell}$, $\ell=1,\ldots,{N_{bins}}$, e.g., for ${N_{\Gamma}}=\{1,2,3\}$, a linear, bilinear, or trilinear function, respectively;

• –

$\phi_{j}({\bm{Y}})$ is continuous on $\Gamma$;

• –

$\phi_{j}(\widehat{{\bm{Y}}}_{j})=1$ at the $j$-th node $\widehat{{\bm{Y}}}_{j}$ of the subdivision ${\mathcal{B}}^{\delta}$; and

• –

if $j^{\prime}\neq j$, $\phi_{j}(\widehat{{\bm{Y}}}_{j^{\prime}})=0$ at the $j^{\prime}$-th node $\widehat{{\bm{Y}}}_{j^{\prime}}$ of the subdivision ${\mathcal{B}}^{\delta}$.

For $j=1,\ldots,{N_{nodes}}$, let $S_{j}({\bm{Y}})=\mbox{\em support}\,\{\phi_{j}({\bm{Y}})\}\subset\Gamma$ and let $V_{j}=\mbox{\em volume}\,\{S_{j}({\bm{Y}})\}$; note that $S_{j}({\bm{Y}})$ consists of the union of the bins $\mathcal{B}_{\ell}\in{\mathcal{B}}^{\delta}$ having the node $\widehat{{\bm{Y}}}_{j}$ as one of its vertices. Thus, the basis functions have compact support with respect to $\Gamma$. An illustration of the basis functions in one dimension is given in Figure 1. We further let $N_{bins,j}$, for $j=1,\ldots,N_{nodes}$, denote the number of bins in $S_{j}({\bm{Y}})$, i.e., the number of bins that share the vertex $\widehat{{\bm{Y}}}_{j}$.

Note that the approximating space ${\mathcal{V}}_{\delta}$ and the basis $\{\phi_{j}({\bm{Y}})\}_{j=1}^{N_{nodes}}$ are in common use for the finite element discretization or partial differential equations. Details about the geometric subdivision ${\mathcal{B}}^{\delta}$, the approximation space ${\mathcal{V}}_{\delta}$, and the basis functions $\{\phi_{j}({\bm{Y}})\}$ and their properties may be found in, e.g., [7, 10].

Below we make use of two properties of the basis $\{\phi_{j}(\mathbf{Y})\}_{j=1}^{N_{nodes}}$. First, we have the well-known relation

[TABLE]

We also have that

[TABLE]

Note that, in general, $C_{j}$ is proportional to $V_{j}$.

2.1 The piecewise-linear approximation of a PDF

Given the $M$ samples values $\{{\bm{Y}}_{m}\}_{m=1}^{M}$ in $\Gamma$, we define the approximation $f_{\delta,M}({\bm{Y}})$ $\in{\mathcal{V}}_{\delta}$ of the unknown PDF $f({\bm{Y}})$ given by

[TABLE]

Note that only the samples ${\bm{Y}}_{m}\in S_{j}({\bm{Y}})$, i.e., only the samples in the support $S_{j}({\bm{Y}})$ of the basis function $\phi_{j}({\bm{Y}})$, are used to determine $F_{j}$. We observe that the proposed estimator can be regarded as a kernel density estimator with linear binning [17], where the binning kernel and the kernel associated with a given grid point are equal to the same hat function. With this choice, the smoothing parameter of the kernel becomes the binning parameter $\delta$, so no additional tuning of the bandwidth is necessary.

Of course, the approximate PDF (5) should be a PDF in its own right. That it is indeed a PDF is shown in the following lemma.

Lemma 1

$f_{\delta,M}(\mathbf{Y})\geq 0$* for all $\mathbf{Y}\in\Gamma$ and $\int_{\Gamma}f_{\delta,M}(\mathbf{Y})d\mathbf{Y}=1$.*

Proof 1

Clearly $f_{\delta,M}({\bm{Y}})$ is non-negative because it is a linear combination of non-negative functions with non-negative coefficients.

[TABLE]

The third and fourth equalities hold because of (4) and (3), respectively.

The next lemma is useful to prove the convergence of our approximation.

Lemma 2

Let $f\in C^{2}(\Gamma)$ with $f|_{\partial\Gamma}=0$ and let $\mathbb{E}[F_{j}]$ denote the expectation of $F_{j}$ with respect to $f$. Then

[TABLE]

where the constant $C$ does not depend on either $\delta$ or $M$, and $\alpha$ is a positive integer. If $[-1,1]^{N_{\Gamma}}\cap\Gamma=[-c_{1},c_{1}]^{{N_{\Gamma}}}$ for some positive constant $c_{1}$, then $\alpha=2$. Otherwise $\alpha=1$.

Proof 2

Let $\chi_{[-1,1]}$ be the characteristic function of $[-1,1]$ and define $\phi(Y):=(1-|Y|)\chi_{[-1,1]}$. Let $\phi(\mathbf{Y}):=\prod_{n=1}^{N_{\Gamma}}\phi(Y_{n})$ be defined in the usual tensor product fashion. Assuming $f|_{\partial\Gamma}=0$, we have $C_{j}=\delta^{N_{\Gamma}}$ for all $j$. Then

[TABLE]

$F_{j}$* is the value of a naive kernel density estimator of $f$ evaluated at $\widehat{{\bm{Y}}}_{j}$, with the function $\phi$ as a kernel. Using a standard argument for the bias of kernel density estimators we have that*

[TABLE]

The above inequality proves the result for $\alpha=1$. Thanks to the symmetry of $\phi$, if $[-1,1]^{N_{\Gamma}}\cap\Gamma=[-c_{1},c_{1}]^{N_{\Gamma}}$ for some positive constant $c_{1}$, then $\int_{[-1,1]^{N_{\Gamma}}\cap\Gamma}\phi(\mathbf{Y}^{{}^{\prime}})\mathbf{Y}^{{}^{\prime}}d\mathbf{Y}^{{}^{\prime}}=0$, hence the result follows with $\alpha=2$.

The next theorem shows that the approximate PDF obtained with our method converges to the exact PDF with respect to the $L^{2}(\Gamma)$ norm.

Theorem 1

Let $\Gamma$ be a polytope in $\mathbb{R}^{N_{\Gamma}}$ and $f\in C^{2}(\Gamma)$ with $f|_{\partial\Gamma}=0.$ If $f_{\delta,M}$ is the approximation of $f$ given in (5), then:

[TABLE]

Moreover, if $[-1,1]^{N_{\Gamma}}\cap\Gamma=[-c_{1},c_{1}]^{N_{\Gamma}}$ for some positive constant $c_{1}$, then

[TABLE]

where $C$ is a constant that does not depend on $\delta$ or $M$.

Proof 3

Let $\mathcal{I}_{\delta}f=\sum\limits_{j=1}^{N_{nodes}}f(\widehat{{\bm{Y}}}_{j})\phi_{j}$ be the finite element nodal interpolant of $f$, then

[TABLE]

where $C_{1}$ is a constant that does not depend on $\delta$ [7]. Considering the second term in the above inequality, we have

[TABLE]

The last inequality is obtained using Lemma 2 and (3). Considering that \mathbb{E}[F_{j}]=\mathbb{E}\Big{[}\dfrac{\phi_{j}}{\delta^{N_{\Gamma}}}\Big{]}=\dfrac{1}{\delta^{N_{\Gamma}}}\mathbb{E}\Big{[}\phi_{j}\Big{]}, we have

[TABLE]

where $\sigma(\phi_{j})=\sqrt{\mathbb{E}[\phi_{j}^{2}]-\mathbb{E}[\phi_{j}]^{2}}\leq C_{3}\sqrt{\delta^{N_{\Gamma}}}$ for all $j$, with $C_{3}$ independent of both $\delta$ and $M$. Hence

[TABLE]

so the first result is obtained. If $[-1,1]^{N_{\Gamma}}\cap\Gamma=[-c_{1},c_{1}]^{N_{\Gamma}}$ for some positive constant $c_{1}$, then $\alpha=2$ in (13), so the second result also follows taking the limit as $M\rightarrow\infty$.

We note that the numerical examples considered below show that convergence can be obtained even for cases where the PDF is not in $C^{2}(\Gamma)$, even when the PDF is not differentiable or even continuous.

2.2 Numerical illustrations

In Section 3, we validate our approach by approximating known joint PDFs $f({\bm{Y}})$. Of course, in comparing approximations to an exact known PDF, we pretend that we have no or very little knowledge about the latter except that we have available $M$ samples of the PDF $\{f({\bm{Y}}_{m})\}_{m=1}^{M}$ at $M$ points ${\bm{Y}}_{m}\in\Gamma$, $m=1,\ldots,M$. For the rest of the paper, whenever $M_{1}<M_{2}$, then $M_{1}\subset M_{2}$, meaning that smaller sample sets are obtained as subsets of a larger sample set. Comparing with known PDFs allows us to precisely determine errors in the approximation of the PDF determined using our method. Then, in Section 4, we use our method to approximate the PDFs of outputs of interest associated with the solution of a stochastic partial differential equation; in that case, the PDF is not known. All computations were performed on a Dell Inspiron 15, 5000 series laptop with the CPU {Intel(R) Core(TM) i3-4030U CPU1.90GHz, 1895 MHz} and 8 GB of RAM.

Note that in all the numerical examples, $\Gamma$ denotes a sampling domain, i.e., all samples $\{\bm{Y}_{m}\}_{m=1}^{M}$ lie within $\Gamma$. For most cases, $\Gamma$ is also the support domain for the PDF. However, we also consider the case in which the support of the PDF is not known beforehand so that the sampling domain $\Gamma$ is merely assumed to contain, but not be the same as, the support domain.

For simplicity, the sample space is assumed to be a bounded ${N_{\Gamma}}$-dimensional box $\Gamma=[a,b]^{N_{\Gamma}}$ with $a<b$. We subdivide the parameter domain $\Gamma$ into a congruent set of bins ${\mathcal{B}}^{\delta}=\{\mathcal{B}_{\ell}\}_{k=1}^{N_{bins}}$ consisting of ${N_{\Gamma}}$-dimensional hypercubes of side $\delta=(b-a)/N_{\delta}$, where $N_{\delta}$ denotes the number of intervals in the subdivision ${\mathcal{B}}^{\delta}$ of $\Gamma$ along each of the ${N_{\Gamma}}$ coordinate directions. We then have that the number of bins is given by $N_{bins}=N_{\delta}^{N_{\Gamma}}$ and the number of nodes is given by ${N_{nodes}}=(N_{\delta}$+$1)^{N_{\Gamma}}$. For simplicity, we assume throughout that the components of the random variable ${\bm{Y}}$ are independently distributed so that the joint PDFs are given as the product of univariate PDFs; our method can also be applied in a straightforward way to cases in which the components of ${\bm{Y}}$ are correlated.

3 Validation through comparisons with known PDFs

In this section, we assume that we have available $M$ samples $\{{\bm{Y}}_{m}\}_{m=1}^{M}$ drawn from a known PDF $f({\bm{Y}})$. The error incurred by any approximation ${f}_{approx}({\bm{Y}})$ of the exact PDF $f({\bm{Y}})$ is measured by

[TABLE]

In particular, we use this error measure for our approximation ${f}_{\delta,M}({\bm{Y}})$ defined in (5).

The accuracy of approximations of a PDF, be they by histograms or by our method, depends on both $M$ (the number of samples available) and $\delta$ (the length of the bin edges). Thus, $M$ and $\delta$ should be related to each other in such a way that errors in (14) due to sampling and bin size are commensurate with each other. Thus, if the bin size errors in (14) are of ${\mathcal{O}}(\delta^{r})$ and the sampling error is of ${\mathcal{O}}(M^{-1/2})$, we set

[TABLE]

Thus, once $a$ and $b$ are specified, one can choose $N_{\delta}$ (or equivalently $\delta$) and the value of $M$ is determined by (15) or vice versa. Clearly, $M$ increases as $\delta$ decreases.

For most of the convergence rate illustrations given below, we

[TABLE]

Note that neither (15) or (16) depend on the dimension $N_{\Gamma}$ of the parameter domain but, of course, $N_{bins}=N_{\delta}^{N_{\Gamma}}$ and $N_{nodes}=(N_{\delta}$ + $1)^{N_{\Gamma}}$ do. If the variance of the PDF is large, one may want to increase the size of $M$ by multiplying the term $\delta^{-2r}$ in (15) by the variance.

The computation of the coefficients $F_{j}$ defined in (5) may be costly if $M$ is large and consequently $\delta$ is small. To improve the computational efficiency, we evaluate a basis function $\phi_{j}({\bm{Y}})$ at a sample point ${\bm{Y}}_{m}$ only if the point is within the support of $\phi_{j}({\bm{Y}})$. However, the determination of the bin $\mathcal{B}_{\ell}\in{\mathcal{B}}^{\delta}$ such that ${\bm{Y}}_{m}\in\mathcal{B}_{\ell}$ may be expensive in case of large $M$ and $N_{\Gamma}>2$. For this task, we employ the efficient point locating algorithm described in [8].

3.1 A smooth PDF with known support

For the first example we consider, we ignore the fact that we know the exact PDF we are trying to approximate, but assume we know the support of the PDF so the sampling domain $\Gamma$ is also the support domain. We also use this example to illustrate that the number of Monte Carlo samples needed is independent of the dimension $N_{\Gamma}$ of the parameter space.

We set $\Gamma=[-5.5,5.5]^{N_{\Gamma}}$ so that $(b-a)=11$ and assume that the components of the random vector ${\bm{Y}}=(Y_{1},\ldots,Y_{N_{\Gamma}})$ are independently and identically distributed according to a truncated standard Gaussian PDF so that the joint PDF is given by

[TABLE]

The scaling factor $C_{G}$ is introduced to insure that we indeed have a PDF, i.e., that the integral of $f({\bm{Y}})$ over $\Gamma$ is unity. Note that because the standard deviation of the underlying standard Gaussian PDF is unity, the values of the truncated Gaussian distribution (17) near the faces of the box $\Gamma=[-5.5,5.5]^{N_{\Gamma}}$ are very small so that the results of this example are given to a precision such that they would not change if one considers instead the (non-truncated) standard Gaussian distribution. Also note that because the second moment of the standard Gaussian distribution is unity, the absolute error (14) is also very close to the error relative to the given PDF.

Before we use the formula (16) to relate $M$ and $\delta$, we first separately examine the convergence rates with respect to $\delta$ and $M$. To this end, to illustrate the convergence with respect to $\delta$, we set

[TABLE]

so that the error due to sampling is relatively negligible. For illustrating the convergence with respect to $M$, we set

[TABLE]

so that the error due to the bin size is relatively negligible. The plots for $N_{\Gamma}=1$ in Figure 2 illustrate a second-order convergence rate with respect to $\delta$ and a half-order convergence rate with respect to $1/M$.

We now turn to relating $M$ and $\delta$ using the formula (16). We consider the multivariate truncated standard Gaussian PDF (17) for $N_{\Gamma}=1,2,3$. Plots of the error vs. both $\delta$ and $M$ are given in Figure 3 from which we observe, in all cases, the second-order convergence rate with respect $\delta$ and the half-order convergence rate with respect to $1/M$. We also observe that the errors and the number of samples used are largely independent of the value of $N_{\Gamma}$. A visual comparison of the exact truncated standard Gaussian distribution (17) and its approximation (5) for the bivariate case is given in Figure 4.

Computational costs are reported in Table 1 in which, for each $N_{\Gamma}=1,2,3$, we choose $k=2,3,4,5$ in (16) to determine $\delta$ and $M$. Reading vertically for each $N_{\Gamma}$, we see the increase in computational costs due to the decrease in $\delta$ and the related increase in $M$, although the method scales linearly with respect to the sample size $M$. Reading horizontally so that $\delta$ and $M$ are fixed, the increase in costs is due to the increasing number of bins and nodes as $N_{\Gamma}$ increases. We note that our method is amenable to highly scalable parallelization not only as $\delta$ decreases and $M$ increases, but also as $N_{\Gamma}$ increases so that, through parallelization, our method may prove to be useful in dimensions higher than those considered here. When developing a parallel implementation of our method, using a point locating algorithm such as that of [8] to locate a sample on a finite element grid shared by several processors would be crucial to realizing the gains in efficiency due to parallelization.

3.2 A smooth PDF with unknown support

Still considering a known PDF, we now consider a case for which we not only pretend we do not know the PDF, but also we do not know its support. Specifically, we consider the uniform distribution $f(Y)=0.5$ on $[-1,1]$. A univariate distribution suffices for the discussions of this case; multivariate distributions can be handled by the obvious extensions of what is said here about the univariate case. We assume that we know that the support of the known PDF lies within a larger interval $\Gamma$. Of course, we may be mistaken about this so that once we examine the sample set $\{{\bm{Y}}_{m}\}_{m=1}^{M}$, we may observe that some of the samples fall outside of $\Gamma$. In this case we can enlarge the interval $\Gamma$ until we observe that the interval spanned by smallest to largest sample values is contained within the new $\Gamma$.

We first simply assume that we have determined, either through external knowledge or by the process just described, that the support of the PDF we are considering lies somewhere within the interval $\Gamma=[-1.5,1.5]$. Not knowing the true support, we not only sample in the larger interval $\Gamma$ (so that here we have $(b-a)=3$ and $\delta=3/N_{\delta}$), but we also build the approximate PDF with respect to $\Gamma$. We remark that a uniform distribution provides a stern test when the support of the distribution is not known because that distribution is as large at the boundaries of its support as it is in the interior. Distributions that are small near the boundaries of their support, e.g., the truncated Gaussian distribution of Section 3.1, would yield considerably smaller errors and better convergence rates compared to what are obtained for the uniform distribution. Choosing $k=2,3,4,5$ and $r=2$ in (16), we obtain the errors plotted in Figure 5. Clearly, the convergence rates are nowhere near optimal. Of course, the reason for this is that by building the approximation with respect to $\Gamma$, we are not approximating the uniform distribution on $[-1,1]$, but instead we are approximating the discontinuous distribution

[TABLE]

For comparison purposes we provide, in Figure 6, results for the case where we use the support interval $[-1,1]$ for both sampling and for approximation construction. Because now the PDF is smooth, in fact constant, throughout the interval in which the approximation is constructed, we obtain optimal convergence rates.

One can improve on the results of Figure 5, even if one does not know the support of the PDF one is trying to approximate, by taking advantage of the fact that the samples obtained necessarily have to belong to the support of the PDF and therefore provide an estimate for that support. For instance, for the example we are considering, one could proceed as follows.

For a chosen $M$, sample $\{Y_{m}\}_{m=1}^{M}$ over $[-1.5,1.5]$. 2. 2.

Determine the minimum and maximum values $Y_{min}$ and $Y_{max}$, respectively, of the sample set $\{Y_{m}\}_{m=1}^{M}$. 3. 3.

Choose the number of bins $N_{bins}$ and set $\delta=(Y_{min}-Y_{max})/N_{bins}$. 4. 4.

Build the approximation over the interval $[Y_{min},Y_{max}]$ with a bin size $\delta$.

It is reasonable to expect that as $M$ increases, the interval $[Y_{min},Y_{max}]$ becomes a better approximation to the true support interval $[-1,1]$. Figure 7 illustrates the convergence of $[Y_{min},Y_{max}]$ to $[-1,1]$. Note that because $[Y_{min},Y_{max}]\subset[-1,1]$, the exact PDF is continuous within $[Y_{min},Y_{max}]$. Thus, it is also reasonable to expect that because the approximate PDF is built with respect an interval which is contained within the support of the exact PDF, that there will be an improvement in the accuracy of that approximation compared to that reported in Figure 5 and, in particular, that as one increases $N_{bins}$ so that $\delta$ decreases and $M$ increases, better rates of convergence will be obtained. Figure 8 corresponds to the application of this procedure and shows the substantially smaller errors and substantially higher convergence rates compared that reported in Figure 5.

A visual comparisons of the approximations obtained using the smallest $\delta$/largest $M$ pairing corresponding to Figures 5, 6, and 8 are given in Figure 9. The defects resulting from the use of the interval $[-1.5,1.5]$ for constructing the approximation of a uniform PDF that has support on the interval $[-1,1]$ are clearly evident. On the other had, using the support interval approximation process outlined above results in a visually identical approximation as that obtained using the correct support interval $[-1,1]$. Note that for the smallest value of $\delta$, we have that $Y_{min}$ approximates $-1$ and $Y_{max}$ approximates $1$ to seven decimal places.

3.3 A non-smooth PDF

We next consider the approximation of a non-smooth PDF. Specifically, we consider the centered truncated Laplace distribution

[TABLE]

over $\Gamma=[-5.5,5.5]$, where $C_{L}=1-\exp(-5.5/1.5)$ is a scaling factor that ensures a unitary integral of the PDF over $\Gamma$. Here, the support domain and sampling domain are the same. This distribution is merely continuous. i.e., its derivative is discontinuous at $Y=0$, so one cannot expect optimally accurate approximations. However, as illustrated in Figure 10, it seems the approximation does converge, but at a lower rate with respect to $\delta$ and at the optimal rate with respect to $M$. The latter is not surprising because Monte Carlo sampling is largely impervious to the smoothness or lack thereof of the function being approximated.

Whenever there is any information about the smoothness of the PDF, one can choose an appropriate value of $r$ in (16). Alternately, possibly through a preliminary investigation, one can estimate the convergence rate of the approximation (5). In the case of the Laplace distribution which is continuous but not continuously differentiable, one cannot expect a convergence rate greater than one. Selecting $r=1$ in (16) to relate $M$ and $\delta$, we obtain the results given in Figure 11 which depicts rates somewhat worse that we should perhaps expect.

The Laplace distribution, although not globally $C^{2}$, is piecewise smooth, with failure of smoothness only occurring at the symmetry point of the distribution. For example, for the particular case of the centered distribution (20), the distribution is smooth for $Y>0$ and $Y<0$. Thus, in general, one could build two separate, optimally accurate approximations, one for the right of the symmetry point and the other for the left of that point. Of course, doing so requires knowledge of where that point is located. If this information is not available, then one can estimate the location of that point by a process analogous to what we described in Section 3.2 for distributions whose support is not known a priori. Such a process can be extended to distributions with multiple points at which smoothness is compromised.

3.4 Bivariate mixed PDF

We now consider a bivariate PDF in which the random variables $Y_{1}$ and $Y_{2}$ are independently distributed according to different PDFs. Specifically, we have that $Y_{1}$ is distributed according to a truncated Gaussian distribution with zero mean and standard deviation $2$, whereas $Y_{2}$ is distributed according to a truncated standard Gaussian. We choose $\Gamma=[-5.5,5.5]^{2}$ so that the joint PDF is given by

[TABLE]

where $C_{G}$ is as in (17) and $C^{\prime}_{G}=1/2(\mbox{erf}(2.75/\sqrt{2})-\mbox{erf}(-2.75/\sqrt{2}))$. Results for this case are shown in Figure 12, where we observe optimal convergence rates with respect to both $\delta$ and $M$. Visual evidence of the accuracy of our approach is given in Figure 13 that shows the approximation of the exact PDF (21) and zoom-ins of the approximate and approximate PDFs. Computational times are very similar to those for $N_{\Gamma}=2$ in Table 1 so that they are not provided here.

4 Application to an unknown PDFs associated with a stochastic PDE

In this section, we consider the construction of approximations of the PDF of outputs of interest that depend on the solution of a stochastic PDE. In general, such PDFs are unknown a priori.

The boundary value problem considered is the stochastic Poisson problem

[TABLE]

where $D\subset\mathbb{R}^{d}$ denotes a spatial domain with boundary $\partial D$ and $\Gamma_{input}\subset\mathbb{R}^{N_{\Gamma_{input}}}$ is the sample space for the input random vector variable ${\bm{Z}}$ which we assume is distributed according to a known input joint PDF $f_{input}({\bm{Z}})$.

For the coefficient function $\kappa(\textbf{x},{\bm{Z}})$, we assume that there exists a positive lower bound $\kappa_{\min}>0$ almost surely on $\Gamma_{input}$ for all $\textbf{x}\in D$. We also assume that $\kappa(\textbf{x},{\bm{Z}})$ is measurable with respect to ${\bm{Z}}$. It is then known that the system (22) is well posed almost surely for ${\bm{Z}}\in\Gamma_{input}$; see, e.g., [2, 16, 23, 24] for details.

Stochastic Galerkin approximation of the solution of the PDE

We assume that we have in hand an approximation $u_{approx}(\textbf{x},{\bm{Z}})$ of the solution $u(\textbf{x},{\bm{Z}})$ of the system (22). Specifically, spatial approximation is effected via a piecewise-quadratic finite element method [7, 10]. Because this aspect of our algorithm is completely standard, we do not give further details about how we effect spatial approximation. For stochastic approximation, i.e., for approximation with respect to the parameter domain $\Gamma_{input}$, we employ a spectral method. Specifically, we approximate using global orthogonal polynomials, where orthogonality is with respect to $\Gamma_{input}$ and the known PDF $f_{input}({\bm{Z}})$. Thus, if $\{\Phi_{j}({\textbf{x}})\}$ denotes a basis for the finite element space used for spatial approximation and $\{\Psi_{i}({\bm{Z}})\}$ denotes a basis for the spectral space used for approximation with respect to ${\bm{Z}}$, with the stochastic Galerkin method (SGM) we obtain an approximation of the form

[TABLE]

Note that once the SGM approximation (23) is constructed, it may be evaluated at any point ${\textbf{x}}\in D$ and for any parameter vector ${\mathbf{Z}}\in\Gamma_{input}$. Having chosen the types of spatial and stochastic approximations we use, the approximate solution $u_{approx}(\textbf{x},{\bm{Z}})$, i.e., the set of coefficients $\{U_{i,j}\}$, is determined by a stochastic Galerkin projection, i.e., we determine an approximation to the Galerkin projection of the exact solution with respect to both the spatial domain $D$ and the parameter domain $\Gamma_{input}$. This approach is well documented so that we do not dwell on it any further; one may consult, e.g., [3, 4, 9, 15, 16], for details. Note that once the surrogate (23) for the solution of the PDE is built, it can be used through direct evaluation to cheaply determine an approximation of the solution of the PDE for any ${\textbf{x}}\in\ D$ and any ${\bm{Z}}\in\Gamma_{input}$ instead of having to do a new approximate PDE solve for any new choice of x and ${\bm{Z}}$.

The error in $u_{approx}(\textbf{x},{\bm{Z}})$ depends on the grid-size parameter $h$ used for spatial approximation and the degree of the orthogonal polynomials used for approximation with respect to the input parameter vector ${\bm{Z}}$. In practice, these parameters should be chosen so that the two errors introduced are commensurate. However, here, because our focus is on stochastic approximation and because throughout we use the same finite element method for spatial approximation, we use a small enough spatial grid size so that the error due to spatial approximation is, for all practical purposes, negligible compared to the errors due to stochastic approximation.

Outputs of interest depending on the solution of the PDE

In our context, outputs of interest are spatially-independent functionals of the solution $u(\textbf{x},{\bm{Z}})$ of the system (22). Here, we consider the two specific functionals

[TABLE]

or

[TABLE]

i.e., the average of the approximate solution and one involving the integral of the square of $u_{i}$, respectively, over the spatial domain $D$. The output of interest $Y$ is a random variable that depends on the random input vector ${\bm{Z}}$. Note that although we consider scalar outputs of interest, the extension to vector-valued outputs of interest is straightforward. Note that throughout, all outputs of interest are standardized, namely, they are translated by their mean and scaled by their standard deviation. Hence, we seek approximations of the PDFs of standardized outputs of interest.

It is important to keep in mind that we are dealing with two random variables. First, we have the input random variable ${\mathbf{Z}}$ having a known PDF $f_{input}({\mathbf{Z}})$ supported over the known parameter domain $\Gamma_{input}$. Second, we have the output random variable $Y$ having an unknown PDF $f_{output}(Y)$ supported over an unknown parameter domain $\Gamma_{output}$. Although we do not know the output PDF, in fact that is what we want to construct so that further samples of the output of interest $Y$ can be obtained by simple direct sampling of the PDF $f_{output}(Y)$ for $Y$.

Thus, the task at hand is

given the known PDF $f_{input}({\bm{Z}})$ of the random input ${\bm{Z}}\in\Gamma_{input}$, determine an approximation of the unknown PDF $f_{output}(Y)$ of an output of interest $Y({\bm{Z}})$.

To deal with this task, one simply follows the recipe:

construct the stochastic Galerkin approximation $u_{approx}(\textbf{x},{\bm{Z}})$ given in (23) of the solution $u(\textbf{x},{\bm{Z}})$ of the PDE (22); 2. 2.

choose $M$ samples $\{{\bm{Z}}_{m}\}_{m=1}^{M}$ of the input random vector ${\bm{Z}}$ according to the known given input PDF $f_{input}({\bm{Z}})$; 3. 3.

determines $M$ samples of the approximate solution $\{u_{approx}(\textbf{x},{\bm{Z}}_{m})\}_{m=1}^{M}$ of the PDE (22) by evaluating (23) at each of the samples ${\bm{Z}}_{m}$ chosen in step 2; 4. 4.

use the approximate solution samples obtained in step 3 to determine $M$ samples $\{Y_{m}=Y({\bm{Z}}_{m})\}_{m=1}^{M}$ of an output of interest from, e.g., (24) or (25); 5. 5.

use the output of interest samples $\{Y_{m}\}_{m=1}^{M}$ obtained in step 4 to determine, using (5), an approximation to the output PDF f_{output}(Y)=f_{output}\big{(}Y({\bm{Z}})\big{)}.

Of course, because the exact PDF $f_{output}(Y)$ is not known, we cannot use (14) to compute errors. Thus, as a surrogate for the exact PDF, we use a histogram approximation $f^{histo}_{{\widehat{M}},{\widehat{\delta}}}(Y)$ obtained with a large number of bins (and therefore a very small ${\widehat{\delta}}$ and a large number of samples ${\widehat{M}}$), where “large” is relative to what is used in obtaining approximations using, e.g., (5). Thus, we now use

[TABLE]

as a measure of the error in any approximation ${f}_{approx}(Y)$ of ${f}_{output}(Y)$ that involves $M\ll{\widehat{M}}$ samples and a bin width $\delta\gg{\widehat{\delta}}$.

Illustrative results for a specific choice for the coefficient of the PDE

For the coefficient function in the PDE (22), we choose

[TABLE]

where $\{\lambda_{n},\psi_{n}({\textbf{x}})\}$ are the eigenpairs of a given covariance function, with the eigenvalues arranged in non-increasing order and the random variables $\{Z_{n}\}_{n=1}^{N_{input}}$ are independent and identically distributed standard Gaussian variables. One recognizes that $\gamma(\textbf{x},{\bm{Z}})$ is a truncated Karhunen-Loève (KL) expansion corresponding to a correlated Gaussian random field with mean $\mu$ [13, 21, 26].

Here, we consider the specific covariance function

[TABLE]

where $\sigma_{\gamma}^{2}$ denotes a variance and $0<L\leq\mbox{diam}(D)$ a correlation length. The eigenpairs satisfy the generalized eigenvalue problem

[TABLE]

The eigenpairs $\{\lambda_{n},\psi_{n}\}$ are approximately determined by means of a Galerkin projection as in [9, 16].

The $N_{input}$ components of the input stochastic variable ${\bm{Z}}$ are independent and are all distributed according to a standard Gaussian PDF. Thus, in the spectral method discretization with respect to the input stochastic variable ${\bm{Z}}$, we use tensor products of Hermite polynomials as the orthonormal basis $\{\Psi_{j}({\mathbf{Z}})\}$. Due to the orthogonality of these polynomials with respect to the Gaussian PDF, this choice results in very substantial savings in both the assembly and solution aspects of the discretized SGM system. Details can be found in, e.g., [9, 16].

For the numerical tests, we choose the spatial domain $D$ to be the unit square $[0,1]^{2}$, the correlation length $L=0.1$, $\kappa_{\min}=0.01$, and $\mu=0$.

Output of interest (24).* Here, we consider the approximation (5) of the PDF $f_{output}(Y)$ of the standardized output of interest $Y$ given by (24).*

We consider $\Gamma_{output}=[-5,3]$, and in (27) and (28) we set $N_{input}=2$, $\sigma_{\gamma}=1.4$. Because ${\bm{Z}}=(Z_{1},Z_{2})$ with $Z_{1}$ and $Z_{2}$ being standard Gaussian variables, $\Gamma_{input}=(-\infty,\infty)^{N_{input}}$. A plot of the output of interest (24) as a function of the input variables ${\bm{Z}}=(Z_{1},Z_{2})$ is given, for ${\bm{Z}}\in[-3.5.3.5]^{2}$, in Figure 14 (left). Plots of the approximate output PDF $f_{\delta,M}(Y)$ determined using (5) is given in Figure 15. For comparison purposes, a plot of the histogram approximation $f^{histo}_{\widehat{\delta},\widehat{M}}(Y)$ determined using (1) is also provided in that figure, but with larger sample size $M$ and smaller bin size $\delta$ compared to those used for $f_{\delta,M}(Y)$.

We next examine the convergence behavior of the approximation (5) of the output PDF. Because the exact PDF is unknown, we measure the error using (26) with the histogram surrogate $f^{histo}_{\widehat{\delta},\widehat{M}}$ obtained with a bin size of $\widehat{\delta}=8/2^{8}$ and $\widehat{M}=10^{8}$ samples. To study the order of convergence with respect to $\delta$, we choose

[TABLE]

whereas for the convergence with respect to $M$, we choose

[TABLE]

Note that for these values, we have $M<{\widehat{M}}$ and $\delta>{\widehat{\delta}}$. In Figure 16, we observe that the convergence rates with respect to $\delta$ and $M$ are approximately $1.75$, and $0.45$, respectively. Given that, from examining Figure 15, the output PDF seems to be $C^{2}(\Gamma_{output})$, these rates are lower than the values $2$ and $0.5$, respectively, that one might expect. Likely causes of these lower rates are that errors are determined by comparing to a histogram approximation and not to the an exact PDF and also because the values of $\widehat{M}$ and $\widehat{\delta}$ used for the histogram surrogate are “close” to the corresponding values used to determine the approximation (5).

Further results about errors and convergence rates are given in Figure 17 for which (16) with $r=2$ is used to relate $\delta$ and $M$. The histogram used for comparison to estimate errors is obtained with $\widehat{\delta}=0.125$ and $\widehat{M}=16^{6}$ samples. The results in this figure are consistent with those of Figure 16.

Output of interest (25). We next consider the approximation (5) of the PDF $f_{output}(Y)$ of the standardized output of interest $Y$ given by (25). The various inputs are the same as those used for the output of interest (24) except that now $\Gamma_{output}=[-3.5,4]$ and $\sigma_{\gamma}=2$. A plot of this output of interest as a function of the input random variables ${\bm{Z}}\in[-3.5.3.5]^{2}$ is given in the right plot of Figure 14. Plots of the approximate output PDF $f_{\delta,M}(Y)$ determined using (5) is given in Figure 18. For comparison purposes, plots of the histogram approximation $f^{histo}_{\widehat{\delta},\widehat{M}}(Y)$ determined using (1) are also provided in that figure, but with larger sample size $M$ and smaller bin size $\delta$ compared to those used for $f_{\delta,M}(Y)$. Figure (19) provides plots of the errors in the approximation (5) determined through comparisons with histogram approximations determined with $\widehat{\delta}=0.1171875$ and $\widehat{M}=16^{6}$. Convergence rates of $1.75$ and $0.4$ are observed with respect to $\delta$ and $M$, respectively.

5 Concluding remarks

A piecewise-linear density estimation method (5) for the approximation of the PDF associated with a given set of samples is presented. The approximation is naturally scalable with respect to the sample size, is intrinsically positive, and has a unitary integral. It is also consistent, meaning that it converges to the exact PDF in the $L^{2}$ norm if the sample size goes to infinity and the bin size goes to zero. The construction of the approximation does not require the solution of a linear system and is fast even for a large number of samples. The computational time has been shown to scale linearly with the sample size. Moreover, the binning size is the only smoothing parameter involved, and it can be related to the sample size by a simple rule.

Future work would involve strategies to extend the method to higher dimensions and decrease the computational time. For instance, the finite element basis functions may be replaced with sparse-grid basis functions in the way done in [25], in order to achieve scalability also with respect to the sample dimension. Computational speed-ups could be obtained by considering parallelization. Another means to speed up the computation would be to use adaptive refinement to coarsen the binning subdivision near the tails of the distribution. With this technique, fewer bins would have to be checked by the point locating algorithm of [8] that we employ to locate a given point in the binning subdivision shared by several processors. In addition, the Monte Carlo sampling used in our method can be replaced by, e.g., quasi-Monte Carlo or sparse-grid sampling, resulting in efficiency gains.