How to reconstruct gPCE response surface in case of numerically generated polynomials

[sitadrost]

EDIT: This question is an extended version of a question I asked last week. I've been experiencing some trouble with my GitHub account, couldn't find my question, and decided to re-post and extended version. Sorry about the duplicate, I'll try to remove the earlier version!

Following up on an earlier question I asked here, I was wondering how to reconstruct a gPCE response surface fitted by Dakota, in case numerically-generated orthogonal polynomials were used. Would that be possible based on the coefficients returned by Dakota?

To illustrate what I mean, say I have the coefficients below for response variable R:

Coefficients of Polynomial Chaos Expansion for R:
        coefficient   u1   u2   u3   u4   u5   u6   u7
        ----------- ---- ---- ---- ---- ---- ---- ----
   2.4062125336e+00 Num0 Num0 Num0 Num0 Num0 Num0   P0
  -2.6915310518e-02 Num1 Num0 Num0 Num0 Num0 Num0   P0
  -1.7194577501e-01 Num0 Num1 Num0 Num0 Num0 Num0   P0
  -2.4893045389e+01 Num0 Num0 Num1 Num0 Num0 Num0   P0
   2.4541212354e-04 Num0 Num0 Num0 Num1 Num0 Num0   P0
  -8.5396586596e+01 Num0 Num0 Num0 Num0 Num1 Num0   P0
   4.5426729065e-01 Num0 Num0 Num0 Num0 Num0 Num1   P0
   3.2663326749e+00 Num0 Num0 Num0 Num0 Num0 Num0   P1

Here, u7 has uniform distribution, so Legendre polynomials were used. Because I specified bounded normal distributions for u1 to u6 and didn't specify which scheme to use, Dakota automatically used the Extended option, with numerically generated orthogonal polynomials. Can I use this information to generate plots like the one attached/below (random example, colours represent value of response variable, white dots are points that were used for generating the gPCE, black dots are points that Dakota meant to use, but where the black box CFD simulation returned NaN)?

From experimenting with the rosen_uq_pce example (thanks for the tip @mseldre !), I managed to reconstruct the response surface using uniform distributions for both variables (Legendre polynomials), normal distributions (Hermite polynomials), but not in case of bounded normal distributions (numerically generated orthogonal polynomials).

The information in the manual is helpful, but describes several different methods:

Extended (default if no option is selected): The Extended option avoids the use of any nonlinear variable transformations by augmenting the Askey approach with numerically-generated orthogonal polynomials for non-Askey probability density functions. Extended polynomial selections with numerically-generated polynomials that are orthogonal to the prescribed probability density functions replace each of the sub-optimal Askey basis selections for bounded normal, lognormal, bounded lognormal, loguniform, triangular, gumbel, frechet, weibull, and bin-based histogram.

Numerically generated orthogonal polynomials
If all random inputs can be described using independent normal, uniform, exponential, beta, and gamma distributions, then Askey polynomials can be directly applied. If correlation or other distribution types are present, then additional techniques are required. One solution is to employ nonlinear variable transformations such that an Askey basis can be applied in the transformed space. This can be effective as shown in [EWC08], but convergence rates are typically degraded. In addition, correlation coefficients are warped by the nonlinear transformation [DKL86], and simple expressions for these transformed correlation values are not always readily available. An alternative is to numerically generate the orthogonal polynomials (using Gauss-Wigert [Sim78], discretized Stieltjes [Gau04], Chebyshev [Gau04], or Gramm-Schmidt [WB06] approaches) and then compute their Gauss points and weights (using the Golub-Welsch [GW69] tridiagonal eigensolution). These solutions are optimal for given random variable sets having arbitrary probability density functions and eliminate the need to induce additional nonlinearity through variable transformations, but performing this process for general joint density functions with correlation is a topic of ongoing research (refer to “Transformations to uncorrelated standard variables” for additional details).

Also, the answer to this question goes into more detail regarding to how the gPCE is generated, which is very interesting, but I'm not sufficiently familiar with the methods mentioned here and in the manual to understand whether it should be possible to reconstruct a response surface based on the coefficients that Dakota lists in the output file, and if so, how.

I'd be happy to do some background reading/studying if that helps, but for now I find the number of possibly relevant references rather overwhelming. I'd really appreciate it if someone could provide me with some pointers!