## 1 Introduction  
Studying physical systems often involves assessing the influence of random parameters upon a system's behavior and outputs. Such uncertainty quantification (UQ) studies are becoming increasingly popular in a variety of technical domains. Of particular importance to industrial scientists and engineers are so-called non-intrusive UQ methods, where the term \"non-intrusive\" signifies that the often complex and typical proprietary simulation software applied to model physical systems are incorporated into UQ studies without any modifications. Sampling methods such as (quasi-) Monte Carlo (MC) [1] or latin hypercube sampling (LHS) [2] are the most common approaches applied for non-intrusive UQ, as they are easy to implement, straightforwardly parallelizable, and can be applied to problems with a large number of parameters under relatively mild assumptions regarding the regularity of the underlying mathematical model. However, their slow convergence rates render them unattractive or even intractable in cases where high accuracy is demanded and single simulations are highly time-consuming. Computationally efficient alternatives are spectral UQ methods [3, 4], which approximate the functional dependence of a system's quantity of interest (QoI) on its input parameters. The most popular black-box methods of this category employ global polynomial approximations by aid of either Lagrange interpolation schemes [5, 6] or generalized polynomial chaos (gPC) [7], where the latter approach is typically based on least squares (LS) regression [8; 9; 10] or pseudo-spectral projection methods [11; 12; 13], or, less commonly, on interpolation [14]. Assuming a smooth input-output dependence, spectral UQ methods provide fast convergence rates, in some cases even of exponential order. However, their performance decreases rapidly for an increasing number of input parameters due to the curse of dimensionality [15]. To overcome this bottleneck, state-of-the-art approaches rely on algorithms for the adaptive construction of sparse approximations [8; 9; 10; 16; 17; 18]. Adaptive sparse approximation algorithms for Lagrange interpolation methods typically employ nested sequences of univariate interpolation nodes. While not strictly necessary [19], nested node sequences are helpful for the efficient construction of sparse grid interpolation algorithms. Nested node sequences are readily available for the well studied cases of uniformly or normally distributed random variables (RVs) [20; 21], but not for more general probability measures. In principle, one could derive nested interpolation (or quadrature) rules tailored to an arbitrary probability density function (PDF) [22; 23], however, this is a rather cumbersome task for which dedicated analyses are necessary each and every time a new PDF is considered. Diversely, the main requirement in gPC approximations, either full or sparse, is to find a complete set of polynomials that are orthogonal to each other with respect to the input parameter PDFs. In the case of arbitrary PDFs, such polynomials can be numerically constructed [24; 25; 26]. As a result, despite the fact that according to a number of studies Lagrange interpolation methods enjoy an advantage over gPC-based methods with respect to the error-cost ratio [27; 28; 29], they are abandoned whenever the input parameter PDFs are not uniform, or normal, or closely related distributions, e.g., log-uniform or log-normal. Consequently, the literature lacks comparisons between Lagrange interpolations and gPC-based approximations for other probability measures. This work investigates the use of Lagrange interpolation methods based on the recently introduced weighted Leja node sequences [30] for the approximation of systems with arbitrarily distributed random parameters. Weighted Leja sequences provide interpolation and quadrature nodes tailored to an arbitrary PDF in a seamless way. Moreover, being nested by definition, they allow for the efficient construction of sparse grids, thus enabling high-dimensional UQ studies.