Reporting uncertainty from a Monte Carlo simulation

Question

I am confused about the how to report uncertainty in Monte Carlo simulations. Take this simplified example:

Say I want to model a system with the following system equations:

X = A*B
Y = C*D
Z = X+Y

If all the inputs (a,b,c,d) have uncertainty around them as represented by probability density functions (e.g. triangular PDFs), I have been trying to use Monte Carlo analysis to propagate this uncertainty by randomly sampling from each pdf and multiplying them together as in the equations above. I then calculate the means, medians and 5th, 95th percentile. For example, in R:

library(triangle)

##  Number of iterations
n <- 10000

##  randomly sample inputs using triangular distribution
A <- rtriangle(n,1,40,10)
B <- rtriangle(n,5,10,8)
X <- A*B

C <- rtriangle(n,0,5,3)
D <- rtriangle(n,5,20,10)
Y <- C*D

Z <- X+Y

## histograms show outputs are not normally distributed
hist(X)
hist(Y)
hist(Z)

outputs <- data.frame(cbind(X,Y,Z))

## Calculate the mean, median and quantiles for X,Y,Z
means <- apply(outputs,2,mean)
median <- apply(outputs,2,median)
lower_quant <- apply(outputs,2,quantile,0.05)
upper_quant <- apply(outputs,2,quantile,0.95)

Am I correct in thinking that the percentiles I have calculated for each output (x,y,z) should be termed prediction intervals (i.e. 90% of all values fall in this range)?

How about the confidence intervals? From what I understand the confidence interval represents the uncertainty around the mean values. Therefore, does this mean in order to get confidence intervals I would have to repeat the above sampling process n times, and then calculate the confidence intervals based on the different means calculated (similar to what I think this post is suggesting I think)?

Is this necessary or even correct procedure, particularly when the distributions of the outputs are not normally distributed?

Many thanks in advance for your help.

Answer 1

This is mostly a stats question rather than a programming question, so I'll answer that.

You really have a probability model here for the 7 quantities A, B, C, D, X, Y and Z. If your system equations are exact, then the distribution of these is determined by the joint distribution of A, B, C and D, which you have assumed are independent with known PDFs. So a Monte Carlo method for studying the joint distribution of all 7 is simply to do what you did: simulate A, B, C, D and calculate the outputs from them.

The quantiles of the outputs won't give you what are usually called "prediction intervals", which usually involve a model containing unknown parameters. I'd simply call them "uncertainty intervals", but "prediction intervals" isn't really too misleading.

On the other hand, "confidence intervals" always involve unknown parameters of a distribution. Since your distribution is completely specified, you wouldn't normally calculate confidence intervals.

What you are doing is approximating the mean and quantiles of this specified distribution using Monte Carlo integration. It appears to me that you are doing this correctly. The fact that the distributions are not normal is irrelevant: you know exactly what they are, and are just using a particular numerical method to calculate them.

There is one place you might want confidence intervals: to calculate the uncertainty in your numerical approximation. The idea is to think of the mean of one of your outputs, say E(Z), as an unknown parameter estimated by mean(Z). Your n realizations of Z from the simulation could be used to estimate the uncertainty in it. Since your simulations are independent, a reasonable estimate of the standard error in mean(Z) is sd(Z)/sqrt(n). You can make n as large as you want to reduce this to a tolerable level. I'd call this "Monte Carlo uncertainty", but some people would just call mean(Z) +/- 1.96 sd(Z)/sqrt(n) an approximate 95% confidence interval for E(Z). Since your n is so large, mean(Z) will be very nearly normally distributed even if Z isn't, and you can use normal theory to get the interval.