# How to use Bayesian updating to measure deviation from initial estimates and how delayed a Scrum project is?

Imagine a Product Owner estimates that the software development team using Scrum they are working with is going to be able to complete 80 user stories in 5 sprints:

• Sprint 1: 16 stories
• Sprint 2: 16 stories
• Sprint 3: 16 stories
• Sprint 4: 16 stories
• Sprint 5: 16 stories

Based on past team performance metrics, the team Product Owner has stated a belief that there is only a 5% chance the team will need more sprints to complete those 80 stories.

But after the first 3 sprints, things look like this:

• Sprint 1: 8 stories
• Sprint 2: 8 stories
• Sprint 3: 12 stories

Is there a way to use Bayes rule to update the original idea that there was only a 5% chance that the 80 stories will be finished on after the 5th sprint? What is the chance now? Would it had been possible to calculate it after the first sprint differed from the estimate?

Assuming sprints behave as they have behaved in the past... What is the way to calculate the number of sprints that are needed to finish the project so that there is only a 5% chance that more sprints are needed?

• Taking aside statistics, if you have an upfront plan for the next five "sprints" it's not Scrum, but Waterfall in disguise. What you are describing has nothing to do with Scrum.
– Tim
Apr 12 at 11:38
• @Tim as far as I know using expected velocity to estimate when a project is going to end is common practice in Scrum projects. I see no reason why using full probability curves instead of just mean velocity makes it “not scrum” Apr 13 at 23:22
• @Tim Scrum allows for “Product Owner to revise forecasted delivery timelines based on variations on velocity of the Development team” scrum.org/resources/blog/agile-metrics-velocity Apr 13 at 23:30
• @Tim he latest version of the Scrum guide only says “various practice exist to forecast progress” explicitly leaving the forecasting approach unspecified, with only rule being that “only what has already happened may be used for forward looking discussions” if that is what you mean, then I agree, goal here is precisely that. Apr 14 at 0:09
• I don't mean forecasting but planning everything upfront. What you are describing is Waterfall, where you first plan, make estimates, then execute.
– Tim
Apr 14 at 4:52

First, your question is very interesting from a curiosity perspective, but I doubt the answer will help the team get more consistent. It could show them they were wrong, but it does not mean they will learn to get better estimates. It might just add tension to the already late project.

Instead of bringing them measurement of how wrong they were, you could try to bring them a possible solution; Evidence Based Scheduling, by Joel Spolsky is an excellent article about how to incorporate probabilities and predictions into time estimates, and how to manage consistent but bad estimates.

Bayes Rule: can we update our belief?

Lets define two variables

• $D$ is the "Is the project Delayed?" variable, and the team stated $p(D) = 0.05$.
• $B$ is "What happened Before?", $B = (S_1 = 8 \cap S_2 = 8 \cap S_3 = 12)$, where $S_i$ is the number of stories done in week $i$.

Assuming the 5% estimate is "the truth", can we refine our information knowing the start of the project?

Bayes rules tells us that

$\underbrace{p(D|B)}_{\text{What you want to know}}p(B) = p(B|D)\underbrace{p(D)}_{\text{What you know}}$.

From the following set of equation,

• $p(D|B)p(B) = p(B|D)p(D)$
• $p(\neg D|B)p(B) = p(B|\neg D)p(\neg D)$, where $\neg D$ is "The project is not delayed".
• $p(D|B) + p(\neg D | B) = 1$, because there is no other possibility than being delayed or not.

We can get this expression,

$p(D|B) = \frac{p(B|D)p(D)}{p(B|D)p(D) + p(B|\neg D)p(\neg D)}$

Estimating $p(B|D)$: we know $p(D)$, so we have to find $p(B|D)$ and $p(B|\neg D)$, and we would need to estimate it from the data.

Why it is hard: without any assumption, we do not know the probability distribution of $S_4$ and $S_5$. We do not know how $S_{1,2,3}$ might influence $S_{4,5}$. The main different with the coin experiment (Wikipedia) you cite in your comment is the following; For the coin, each flip samples the same underlying distribution. The probability that the coin lands head does not depend on the last flip. Your data, on the other hand, is a time series. It is likely that the output of the team in one sprint influences it's output in another sprint. It is also likely that external factors that are not measured alter their output. In this sense, each coin flip is an experiment but measuring the output at a single sprint is not. Ideally, for this kind of task, you would want to have reports of past project, and the number of story done at each week and try to model the trend. If your objective is to predict $S_4$ and $S_5$, you do not have any exemple or any data, just assumptions.

But we can work with assumptions. The assumption you gave in your original question was that "Sprints behave the same as they did in the past". One interpretation would be: At any given week, the teams output follows the same probability distribution, let's say a Normal distribution $\mathcal{N}(\mu,\sigma^2)$, and each week is independent from all the others. Let's roll with it.

As a warm up, and to illustrate the point that three data points is not a lot, we will compute the probability of delay the data is telling us, under this assumption, $p_{\text{data}}(D)$, ignoring the predicted probability of delay $p_{\text{pred}}(D)$.

We can estimate the parameters (Wikipedia) of the distribution.

• The sample mean is $\hat{\mu} = \text{average}(8,8,12) = 9.\overline{3}$ and the sample variance is $s^2 = \frac{1}{2} \sum (S_i - \hat{\mu})^2 = 5.\overline{3}$.
• We can compute the 0.95 confidence interval (Wikipedia) of these values, with
• $\mu \in [\hat{\mu} - t_{n-1,0.975}\frac{1}{\sqrt{3}}s, \hat{\mu} + t_{n-1,0.975}\frac{1}{\sqrt{3}}s]$.
• $\sigma^2 \in \left[\frac{(n-1)s^2}{\chi^2_{n-1,0.975}}, \frac{(n-1)s^2}{\chi^2_{n-1,0.025}}\right]$
• the computation ([Mean], [Variance] on Wolfram) gives
• $\mu \in [3.596, 15.07]$.
• $\sigma^2 \in [1.14, 49.383]$

Even when assuming that every sprint comes from the same distribution and is independent from the other sprints, which is not a "trivial" assumption, the precision of our model is not that great since we have only three points.

Calculating $p_{\text{data}}(D)$:

• In the worst case scenario of our 95% confidence model (min $\mu$, min $\sigma$), $S_4$ and $S_5$ are generated from $\mathcal{N}(\mu = 3.596, \sigma^2 = 1.14)$. So $p_{\text{data}}(D) = p(S_4 + S_5 < 52)$, can be phrased as $p(X + Y < 52)$ with $X, Y \sim \mathcal{N}(\mu = 3.596, \sigma^2 = 1.14)$, or $p(X < 52)$, with $X \sim \mathcal{N}(\mu = 7.192,\sigma^2 = 2.28)$ since to add Gaussians we can add their mean and variance. This is certain (Wolfram)

• In the "best" case (max $\mu$, max $\sigma$), we have $p(X < 52), X \sim \mathcal{N}(\mu = 30.14, \sigma^2 = 98.766)$, or 0.986 (Wolfram)

Under the assumption that the number of story completed at each sprint are i.i.d. and follow a Normal distribution, using only the data we got for the first three weeks, we can say with 95% confidence that the probability of finishing the project in time is less that 1.4%.

This assumes that nothing changes between sprint, and will not change until the end of the project. So this estimate might be a little bit pessimistic. Also, we ignored that the team had an initial estimate of 5% of success.

Adding the team estimate If we assume that their estimate is true, that there was indeed a probability of 5% of delay at the beginning, and if we continue to assume that each sprint is independent from the other and follow the same Normal distribution, we have that

$p(S_1 + S_2 + ... < 80) = 0.05$, with $S_1,S_2,... \sim \mathcal{N}(\mu, \sigma^2)$ or $p(X < 80) = 0.05$, with $X \sim \mathcal{N}(5\mu, 5\sigma^2)$

This is our constraint, and we have to find the parameters $\mu,\sigma^2$ that best explain what we observed; $S_1 = 8, S_2 = 8, S_3 = 12$. We want to maximize the likelihood that the model we have generated our observations.

Using the cumulative distribution function of the Normal distribution (Wikipedia), we have that

$p(X < 80) = 0.05 = \frac{1}{2}\left[1 + \text{erf}\left(\frac{80 - 5\mu}{\sqrt{5}\sigma\sqrt{2}}\right)\right]$ $\Rightarrow$ $0.1 = 1 + \text{erf}\left(\frac{80 - 5\mu}{\sqrt{5}\sigma\sqrt{2}}\right)$ $\Rightarrow$ $\frac{80 -5 \mu}{\sqrt{5}\sigma\sqrt{2}} = \text{erf}^{-1}(-0.9)$ $\Rightarrow$ $\sigma = \frac{80 - 5\mu}{\sqrt{5}\text{erf}^{-1}(-0.9)\sqrt{2}}$. (Note that since $\text{erf}^{-1}(-0.9)$ is negative and $\sigma$ is positive, we need $\mu > 16$)

And we want to maximize $p(S_1 = 8)p(S_2 = 8)p(S_3 = 8)$, with $S_1,S_2,S_3 \sim \mathcal{N}(\mu,\sigma)$, with respect to the previous constraints. I wrote a little matlab script (code available at the end) to do so, since the beautifully named erf function does not have a closed formed, and we get the following result:

So the model $\mathcal{N}(\mu = 23.22, \sigma^2 = 9.82^2)$ is the best model to explain the observed data $S_1,S_2,S_3$ while allowing the prediction that the prediction of a 5% chance of delay being true, under the assumptions we made. So, what is the updated probability $p(D|B)$?

$p(D|B) = p(S_4 + S_5 < 52) = p(X < 52)$, with $X \sim \mathcal{N}(2\mu,2\sigma^2)$, and Wolfram says the probability of the project being delayed, assuming the sprints are i.i.d., that the initial estimate was correct and what they delivered in the first three sprints is 99.2%.

Oh, by the way, this new probability is not saying that their initial estimate was wrong, on the contrary, it assumes that it is right! It is an update of it, knowing that the first three sprints went $8,8,12$.

Comment questions:

What confuses me is the "I don't think we can do 26,26 at $S_4,S_5$". It seems intuitively obvious that a human would think it is unlikely after 8,8,12, but under the probability distribution configured with the assumptions you described, it does not seem unlikely.

In both scenarios I described, (Finding the distribution from the data, or from the initial probability estimate + the data), it is very unlikely that the team will do 26,26. But it is not impossible.

• In the [Data-Only] case, we estimate the distribution of the stories/week, and the "best" model for the first three sprints is $\hat{\mu} = 9.\overline{3}$, $\hat{sigma}^2 = 3.555$. Under this model, is it not possible for 26,26 to happen, but the model might be wrong. We have only 3 datapoints, after all. So we compute the confidence interval on $\mu$ and $\sigma^2$ and get $\mu \in [3.596,15.07]$, $\sigma^2 \in [1.14, 49.383]$.

In the most extreme case, $\mathcal{N}(\mu = 15.07,\sigma^2 = 49.383)$, it is possible to do two sprints that give 52 stories, with probability 0.986 of delay, so a probability of 1.4% of success, at best. But we do not know what is the true probability, since there is an incertainty in the model. We can just say that the probability of having no extension is between 0% and 1.4%, given the assumptions we made.

• In the [Data + Estimation] case, we assume that the distribution initially allowed a 95% chance of success; the initial distribution must allow that to happen. We find the distribution that best explain the observed data in the first three sprints, and it seems to be $\mathcal{N}(\mu = 23.33, \sigma^2 = 9.82^2)$. The probability of success is only 0.8%. In other words, the team, by giving a probability of failure, gave an idea of what the distribution is, and if they had hit the mean at each time they would have easily made it. However, since they hit "way under" in the first three sprints, the probability that they will hit "way above" is low.

You wrote "using only the data we got for the first three weeks, we can say with 95% confidence that the probability of finishing the project in time is less that 1.4%" How come, if there is only 1.4% chance that that it is finished in time, there is only 5% chance that is delayed ? Isn't that a contradiction?

I think there is a confusion with the 95% confidence interval of the model and the probability of delay; We build a model of the distribution underlying the sprints output, and we compute the confidence interval on the parameters of the model ($\mu,\sigma^2$). At that point, we do not know the probability that the project is delayed, but we have a probability distribution of the parameters of the model. In order to be completely thorough, we should compute the probability of the value for each parameter, and compute the probability of delay under this model, and sum all of it to get the estimate of the probability. But to "ballpark" the probability of delay, we take the two extreme models, and we get that the probability of delay should be between 98.6% and 100%, if we take the 95% confidence interval on the model we computed.

Conclusion:

• The data you have is more complex than a coin flip. Ideally, you would need more of it to be able to do something useful
• Without more data, assumptions can make something, but if the assumptions are wrong (and they probably are), the results are garbage.
• We've seen how to make an estimate out of data (and assumptions)
• We've seen how to update an estimate, assumed true, with new data (and assumptions)
• But essentially, with three data points, it is not much better than "Hey guys, we went 8,8,12 the last three sprints... I don't think we can do 26,26 for the next two, so let's mail the client".
• If you really want to master time management in software dev. teams using probability and statistics, read something about that

Matlab code:

mu_range = [16,64];
steps = 2^10;

mus = linspace(mu_range(1), mu_range(2), steps);

sigmas = zeros(size(steps));
likelihoods = zeros(size(steps));
for step = 1:steps
mu = mus(step);
sigma = abs((80 - 5*mu)/(sqrt(5)*sqrt(2)*erfinv(-0.9)));

sigmas(step) = sigma;

% Computing the likelihood
likelihood_8 = 1/sqrt(2*pi*sigma^2)*exp(-(8-mu)^2/(2*sigma^2));
likelihood_12 = 1/sqrt(2*pi*sigma^2)*exp(-(12-mu)^2/(2*sigma^2));
likelihoods(step) = likelihood_8^2 * likelihood_12;
end

[maxLikelihood, bestStep] = max(likelihoods);

bestMu = mus(bestStep);
bestSigma = sigmas(bestStep);

figure();
[ax,p1,p2] = plotyy(mus, likelihoods, mus, sigmas);
legend('Likelihood','Value of \sigma');
title(['Maximum likelihood at \mu =', num2str(bestMu), ' & \sigma = ', num2str(bestSigma)]);
xlabel('Value of \mu');
ylabel(ax(1), 'Likelihood');
ylabel(ax(2), '\sigma');

• Thanks! But I don't quite get how it is possible to check if a coin is fair ( en.m.wikipedia.org/wiki/… ) but what I am asking lacks enough data... I have my assumptions and 3 experiments... Can you please elaborate on the difference ? Jan 21, 2016 at 16:31
• @Luxspes : You are right, I was not very clear on that. I updated my answer to adress your point, and showed that you can indeed do it. What I meant was something along the lines of: your data is not as simple as coin tossing, in order to learn something useful from it you would need a lot of it, and without more we would have to make wild assumption, and that it wouldn't be that useful. Jan 22, 2016 at 1:27
• Althought what confuses me on the end is the "I don't think we can do 26,26". Seems like even while for a human is intuitively obvious that 26, 26 is unlikely after 8,8,12, that is not the case if the probability distribution is configured with the assumptions you described... Am I right? Jan 22, 2016 at 2:12
• Also, you wrote: "Under the assumption that the number of story completed at each sprint are i.i.d. and follow a Normal distribution, using only the data we got for the first three weeks, we can say with 95% confidence that the probability of finishing the project in time is less that 1.4%" How come, if there is only 1.4% chance that that it is finished in time, there is only 5% chance that is delayed ? Isn't that a contradiction @Winks? Jan 22, 2016 at 2:22
• @Luxspes I added an answer to both of your question, I hope I understood what was unclear. Jan 22, 2016 at 12:24

There's a simple solution assuming the conjugate Poisson-Gamma model. Sprints have fixed lengths, you are counting the number of done stories per those intervals, so this can be modeled using the Poisson distribution. Gamma is a conjugate prior for Poisson, so if your counts $$X_1, X_2, \dots, X_n$$ from $$n$$ sprints follow the Poisson distribution with unknown rate $$\lambda$$

$$X_i \sim \mathcal{P}(\lambda)$$

then you can assume the Gamma prior for $$\lambda$$ with hyperparameters $$\alpha$$ for the count assumed a priori $$16 \times 5$$ and $$\beta$$ for the assumed number of sprints $$5$$. In such a case, the posterior distribution for $$\lambda$$ follows the Gamma distribution

\begin{align} \alpha' &= \alpha + \sum_{i=1}^n X_i \\ \beta' &= \beta+n \\ \lambda &\sim \mathcal{G}\big(\alpha',~\beta'\big) \end{align}\tag{1}

with expected value

$$E[\lambda] = \frac{\alpha'}{\beta'} = \frac{16 \times 5 + 8 + 8 + 12}{5 + 3} = 13.5$$

But that's "on average", while you are interested in the posterior predictive distribution of the $$X_i$$ counts, that is negative binomial with parameters

$$\tilde X | X_1,X_2,\dots,X_n \sim \mathcal{NB}\big(\alpha', \tfrac{\beta'}{1+\beta'} \big)$$

It has the same mean as the posterior, but a bigger variance. To learn the upper and lower bounds, just check the quantile function.

> qnbinom(c(0.05, 0.5, 0.95), (16 * 5 + 8 + 8 + 12), (5 + 3)/(1 + 5 + 3))
[1]  7 13 20


There's 5% chance that the team would be doing not more than 5 stories per spring and 95% chance that it would be no more than 20 stories, with the median and average around 13 stories. If you need to do 80 stories, this gives between 4 and 11 sprints:

> 80 / qnbinom(c(0.05, 0.5, 0.95), (16 * 5 + 8 + 8 + 12), (5 + 3)/(1 + 5 + 3))
[1] 11.428571  6.153846  4.000000


Given the simple update rule (1) between previous estimates $$\alpha, \beta$$, and the newly observed counts, this can be updated each time the new data arrives.

Finally, with statistical software, you can easily run a Monte Carlo simulation to find out what is the chance that you'll manage to the 80 stories in 5 sprints, by simulating the five draws from the negative binomial distribution.

> mean(replicate(100000, sum(rnbinom(5, (16 * 5 + 8 + 8 + 12), (5 + 3)/(1 + 5 + 3))) >= 80))
[1] 0.0868


and it's around 9%.

As you can see, most of the math is rather straightforward and things like expected values can be calculated using pen and paper, but for the quartiles (the worst 5% case) you would need statistical software.

• Very interesting answer, I am even considering making it the accepted one as it aligns better with the original intent of my question (although @Winks answer indeed opened my eyes to entirely different interpretation of the situation) Apr 15 at 0:09