# 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.

Well, what is the dimensionality of your input? Generally speaking, 32 model runs are not so much! The experimental design could indeed be a random sample, in which case you could use this for Monte Carlo simulation. Quasi-Monte Carlo would be pretty similar with non-random sequences. The error estimates will change then accordingly. There are many other ways to create experimental designs though. Do you have some details on how the inputs were created?

Best regards

• Hi, the dimensionality for this specific example was 5, but in the future I'll have to perform a similar study with a lot more variables (10 or more). The experimental design is not a random sample: explaining it here is not easy because the inputs in my case didn't have classical distributions, anyway I added an example which gives the idea. – DeltaIV Jan 23 '16 at 12:28
• Welcome to the site. Was this intended as an answer to the OP's question, or a comment requesting clarification from the OP? Please only use the "Your Answer" field to provide answers to the original question. You will be able to comment anywhere when your reputation is >50. Since you're new here, you may want to take our tour, which has information for new users. – gung - Reinstate Monica Jan 23 '16 at 13:24