# How to explain simply that the set of runs for Non Intrusive Polynomial Chaos cannot be used as a Monte Carlo sample

I had quite an annoying problem at work, a few days ago. I was doing a forward Uncertainty Quantification analysis using Non Intrusive Polynomial Chaos (NISP) (see for example here). Basically, you have some function $$Y=f(X_1,\ldots,X_n)$$ with $X_1,\ldots,X_n$ independent random variables of known distribution, and you want to summarize the distribution of $Y$. You don't have an analytical expression for $Y$, but for each realization of $\mathbf{X}=(X_1,\ldots,X_n)$ you can run a black box solver (for example a FEM or CFD code) and get the corresponding value of $Y$. This would seem to be a job for Monte Carlo, but there's a caveat: each run of the black box code requires a lot of time to converge, thus you cannot just fire up $10^6$ runs of your code and get a nice MC ensemble with a low sampling error. NISP approaches this issue by building an approximation for $f$ based on an expansion in a series of multivariate polynomials, orthonormal with respect to the input distribution. To do this, you run the code at some points $\mathbf{X_1},\mathbf{X_2},\ldots,\mathbf{X_m}$ , which form a design matrix (see figure 2 on page 9 of the linked paper). The design matrix is completely determined once you specify the input distribution, thus it's not a Monte Carlo ensemble, as you can see from its very regular structure.

Now, I performed the analysis with a matrix containing 32 points, i.e., I ran my black box code 32 times. The resulting approximation was unsatisfactory, and it was clear that we should increase the degree of the approximation, which would have required another 243 runs. My manager balked at the idea, and told me to instead use the 32 runs as a small Monte Carlo ensemble, to compute mean, standard deviation and a few percentiles of $Y$. I told her that would make no sense, because the samples for a Monte Carlo study must be generated with a pseudo-random number generator. She replied that any set of 32 points $\mathbf{X_1},\mathbf{X_2},\ldots,\mathbf{X_{32}}$ has the same "dignity" as any other equally sized set, because they all have the same probability of being extracted in random sampling. I replied: yes that's true, the probability is 0 for all samples, and knowing this doesn't help us in our analysis. The error estimate for Monte Carlo makes sense if the samples are random (and even in that case, $m=32$ would be quite a stretch for the application of the CLT). Using the 32 deterministic samples we can of course compute a sample mean, but we have no idea of the accuracy of the estimate. At this point she got upset and told me that I was a nitpicker. I'm sure this won't fare well on my evaluation. What could have I said, to explain in a simple way that her idea wasn't going to save us time, but would instead hurt the project goals? Did I make some mistakes in my explanation?

PS while writing this question, it came to my mind that I should have mentioned that without (pseudo)random sampling, we cannot even invoke the Law of Large Numbers, and thus state that our estimate should converge to the desired parameter in the limit $m\to\infty$. Other than this, I really have no idea how I could have made my point simpler.

EDIT: a user required a description of the experimental design. I think the idea is well described in the paper I linked, anyway I'll add an example. Suppose you have 3 random inputs, i.e., $\mathbf{X}=(X_1,X_2,X_3)$, and that each $X_i$ is uniformly distributed in $[-1,1]$. The experimental design of NISP, for a very coarse approximation, is then the tensor product of the 2-points Gauss-Legendre quadrature nodes. In other words, it's the 8x3 matrix

[[-0.57735, -0.57735, -0.57735];
[-0.57735, -0.57735,  0.57735];
[-0.57735,  0.57735, -0.57735];
[-0.57735,  0.57735,  0.57735];
[ 0.57735, -0.57735, -0.57735];
[ 0.57735, -0.57735,  0.57735];
[ 0.57735,  0.57735, -0.57735];
[ 0.57735,  0.57735,  0.57735]]


which means that 8 runs of the black box are required. In my case I had more input variables (5) and their distributions were not classical distributions, thus the experimental design was a bit more complicated, but the idea was exactly the same.