# Refining Monte Carlo predictions using observed measurements

I'm trying to build a monte-carlo simulation that can revise it's distribution of outcomes of a project based on observed measurements after the project has started.

I have a few questions about the best way to do this. I'm not a statistician, so please correct me if I am doing something wildly wrong.

For example, let's say I've observed that task x has been selected by person y (whose original 90% CI estimate for the task was [l,h]) , and that y has logged w hours of work to the task.

I can use that data to re-simulate the project under new constraints and compute a new, more accurate, distribution of outcomes.

For example, if w > l, then I know that the lower bound for the time to complete x is now w, not l, and can adjust the distribution used accordingly. However, w is not a 5% lower bound. It's a 0% lower bound (i.e. the limit), so using [w, h] as a 90% CI didn't quite seem correct. As a result I was thinking I could just pick some arbitrarily small number for p(w), say 0.0001, and continue using .05 for p(h) and then generate a new distribution for [w, h] (of course, I would just use the number of deviations for h and w rather than the probabilities).

Is that sound?

What's not immediately clear is what I would do in the case where w > h. I have calibrated estimates with a 90% CI, so I should expect to see this 5% of the time. If I ask: "what do I know in that case", I come up with the following:

1. I know w and I know my arbitrarily low p(w)
2. From the original confidence interval (which assumes a normal distribution), I can determinep(w + sigma).
3. So, I could produce: [w, w + sigma] as an interval, using p(w) and p(w + sigma), and then derive a normal distribution from that (again, just using the z-values).

Is that sound as well?

• Are you familiar with Bayesian inference? – Sycorax says Reinstate Monica Aug 2 '14 at 20:41
• No. I had to look it up. Are you suggestion that I wouldn't need to re-run the simulation, but instead just adjust the output distribution using the formula here: en.wikipedia.org/wiki/Bayesian_inference? – Scott Wisniewski Aug 2 '14 at 21:05
• Well, what's remarkable to me (as a statistician) is that you've already framed your problem in nearly Bayesian terms, wherein you have a basic idea of your quantity of interest and want to iteratively improve upon that idea as you receive more data. Whether you have to reevaluate your simulation or simply do a little algebra depends on the model you decide is "best" for your problem. – Sycorax says Reinstate Monica Aug 2 '14 at 21:10
• I'm also saying that an investment of time to understand Bayesian inference would provide you with a very flexible toolbox for understanding this and many, many other problems. – Sycorax says Reinstate Monica Aug 2 '14 at 21:11
• Thanks for the pointers. I think I may have figured out how to apply Bayesian inference to my case (see my 'answer' below). – Scott Wisniewski Aug 3 '14 at 0:01

Here's my attempt at a (self) answer, based on user777's hints.

At first I was confused by the P(E|H) term in the Bayesian inference formula, thinking P(E|H) is 0 whenever H > E, and I only care about cases where H > E when I'm computing P(H|E).

However, then I reformulated the problem slightly and said "what is the probability that it will take at least H given that it takes at least E". The probably of at least E given at least H is one, so that yields:

$\int {p(H|E)}{dH} = \frac{\int{p(H)}{dH}}{\int^E_0{p(x)}{dx}}$

I believe that means that if I compute a random number using the original normal distribution that I should just be able to multiply it by $\int^E_0{p(x)}{dx}$ to get my new random number.

Rationale: The normal random number should be $inv(\int{p(x)}{dx})(rand())$. To get a random number corresponding to $p(H|E)$ I need $inv\bigl({\int {p(H|E)}{dH}}\bigr)(rand())$, which should be equal to $\int^E_0{p(x)}{dx}*inv\bigl({\int{p(H)dH}}\bigr)$.

That should be exactly equal to normal_rand(sigma, mu)*cdf(w).

Update:

This needs to be truncated somehow, as the df is only valid when x > w. I'm not quite sure how to do that.

Update 2:

I ended up using the info here: http://en.wikipedia.org/wiki/Truncated_distribution to produce a truncated normal distribution. I think I can use that, without the Bayesian update.

• Since the quantity of interest is additional time required to a task, we might expect that the lower bound of this variable is zero (we can't go back in time). The gamma distribution is often used to accomplish this, of which the exponential distribution is a special case. But there are many other distributions with support over positive reals. Conversely, the normal distribution has support over the real line, so any choice of parameters will always have positive probability below zero. – Sycorax says Reinstate Monica Aug 3 '14 at 6:58