# Estimating the probability that a software change fixed a problem

At work we have a hardware device that is failing for some yet to be determined reason. I have been tasked to see if I can make this device not fail by making changes to its software driver. I have constructed a software test bench which iterates over the driver functions which I feel are most likely to cause the device to fail. So far I have forced 7 such failures and the iterations that the device failed on are as follows:

100 22 36 44 89 24 74

Mean = 55.57 Stdev = 31.81

Next, I made some software changes to the device driver and was able to run the device for 223 iterations without failure before I manually stopped the test.

I want to be able to go back to my boss and say "The fact that we were able to run the device for 223 iterations without failure means that my software change has a X% probability of fixing the problem." I would also be satisfied with the converse probability that the device will still fail with this fix.

If we were to assume that the iteration the device fails on is normally distributed, we can say that going 223 iterations without failure is 5.26 standard deviations from the mean which roughly has a 1-in-14 million chance of happening. However, because we only have a sample size of 7 (not including the 223), I'm fairly certain it would be unwise to assume normality.

This is where I think the Student's t-test comes into play. Using the t-test with 6 degrees of freedom, I've calculated that the actual population mean has a 99% probability of being less than 94.

So now my question to you guys is whether or not I am allowed to say with 99% certainty that hitting 223 iterations without failure is a 4.05 sigma event, i.e. $\frac{(223 - 94)}{31.81} = 4.05$ ? Am I allowed to use the 31.81 sample standard deviation in that calculation or is there some other test I should do to get a 99% confidence on what the maximum standard deviation is and then use that in my calculation for how many sigmas 223 really is away from the mean at the 99% confidence level?

Thanks!

UPDATE

The answers I received here are beyond any expectation I had. I truly appreciate the time and thought many of you have put into your answers. There is much for me to think about.

In response to whuber's concern that the data does not seem to follow an exponential distribution, I believe I have the answer for as to why. Some of these trials were run with what I thought would be a software fix but ultimately ended in failure. I would not be surprised if those trials were the 74 89 100 grouping that we see. Although I wasn't able to fix the problem it certainly seems like I was able to skew the data. I will check my notes to see if this is the case and my apologies for not remembering to include that piece of information earlier.

Lets assume the above is true and we were to remove 74 89 100 from the data set. If I were to re-run the device with the original driver and get additional failure data points with values 15 20 23, how would you then compute the exponentially distributed parametric prediction limit at the 95% confidence level? Would you feel that this prediction limit is still a better statistic than assuming independent Bernoulli trials to find the probability of no failure at 223 iterations?

Looking more closely at the wikipedia page on Prediction Limits I calculated the parametric prediction limits at the 99% confidence level assuming unknown population mean and unknown stdev on Excel as follows:

$\bar{X_n} = 55.57$
$S_n = 31.81$
$T_a = T.INV\Bigl(\frac{1+.99}{2},6\Bigr)$

$Lower Limit = 55.57 - 3.707*31.81*\sqrt{1+\frac{1}{7}} = -70.51$
$Upper Limit = 55.57 + 3.707*31.81*\sqrt{1+\frac{1}{7}} = 181.65$

Since my trial of 223 is outside the 99% confidence interval of [-70.51 , 181.65] can I assume with 99% probability that this is fixed assuming that the underlying distribution is the T-Distribution? I wanted to make sure my understanding was correct even though the underlying distribution is most likely exponential, not normal. I have no clue in the slightest how to adjust the equation for an underlying exponential distribution.

UPDATE 2

So I'm really intrigued with this 'R' software, I've never seen it before. Back when I took my stats class (several years ago) we used SAS. Anyway, with the cursory knowledge I gathered from Owe Jessen's example and a bit of help from Google, I think I came up with the following R code to produce the Prediction Limits with the hypothetical dataset assuming an Exponential Distribution

Let me know if I got this right:

fails <- c(22, 24, 36, 44, 15, 20, 23)
fails_xfm <- fails^(1/3)
Y_bar <- mean(fails_xfm)
Sy <- sd(fails_xfm)
df <- length(fails_xfm) - 1
no_fail <- 223

percentile <- c(.9000, .9500, .9750, .9900, .9950, .9990, .9995, .9999)
quantile <- qt(percentile, df)

UPL <- (Y_bar + quantile*Sy*sqrt(1+(1/length(fails_xfm))))^3
plot(percentile,UPL)
abline(h=no_fail,col="red")
text(percentile, UPL, percentile, cex=0.9, pos=2, col="red")


-
The iteration numbers are surely just identifiers, and should not be used as numbers in computations. How many iterations did you do in total in your original trial? –  James Nov 3 '10 at 16:11
@James I could be wrong but if the iterations are simply identifiers then most people will naturally report them in an increasing order like so: 22 24 36 44 ... The fact that the numbers are reported in an apparently random sequence suggests that each number represents the failed iteration number of 7 separate tests. –  user28 Nov 3 '10 at 16:26
I take your point, but I'm still a little unsure of the value of the mean of these numbers. –  James Nov 3 '10 at 16:57
@Srikant: Yes, your interpretation is spot on. To James' question, I think it would be fair to say that with this given device and software test bench, the device fails on average at the 55th iteration. Does that not sit will with you? –  SiegeX Nov 3 '10 at 18:01
@SiegeX Your last approach has the right flavor, but please don't use the formulae on the Wikipedia page: they are only for normally distributed data. It is rare for failure time data to be normal. They typically are positively skewed. This implies the normal theory upper prediction limits can be (way) too low. BTW, a lower prediction limit is irrelevant here--although the fact that it is hugely negative is a clear indicator of how bad the normal theory methods are for these data. To adjust the equation, see the reference I provided in my answer. –  whuber Nov 3 '10 at 19:35
show 1 more comment

This question asks for a prediction limit. This tests whether a future statistic is "consistent" with previous data. (In this case, the future statistic is the post-fix value of 223.) It accounts for a chance mechanism or uncertainty in three ways:

1. The data themselves can vary by chance.

2. Because of this, any estimates made from the data are uncertain.

3. The future statistic can also vary by chance.

Estimating a probability distribution from the data handles (1). But if you simply compare the future value to predictions from that distribution you are ignoring (2) and (3). This will exaggerate the significance of any difference that you note. This is why it can be important to use a prediction limit method rather than some ad hoc method.

Failure times are often taken to be exponentially distributed (which is essentially a continuous version of a geometric distribution). The exponential is a special case of the Gamma distribution with "shape parameter" 1. Approximate prediction limit methods for gamma distributions have been worked out, as published by Krishnamoorthy, Mathew, and Mukherjee in a 2008 Technometrics article. The calculations are relatively simple. I won't discuss them here because there are more important issues to attend to first.

Before applying any parametric procedure you should check that the data at least approximately conform to the procedure's assumptions. In this case we can check whether the data look exponential (or geometric) by making an exponential probability plot. This procedure matches the sorted data values $k_1, k_2, \ldots, k_7$ = $22, 24, 36, 44, 74, 89, 100$ to percentage points of (any) exponential distribution, which can be computed as the negative logarithms of $1 - (1 - 1/2)/7, 1 - (2 - 1/2)/7, \ldots, 1 - (7 - 1/2)/7$. When I do that the plot looks decidedly curved, suggesting that these data are not drawn from an exponential (or geometric) distribution. With either of those distributions you should see a cluster of shorter failure times and a straggling tail of longer failure times. Here, the initial clustering is apparent at $22, 24, 26, 44$, but after a relatively long gap from $44$ to $74$ there is another cluster at $74, 89, 100$. This should cause us to mistrust the results of our parametric models.

One approach in this situation is to use a nonparametric prediction limit. That's a dead simple procedure in this case: if the post-fix value is the largest of all the values, that should be evidence that the fix actually lengthened the failure times. If all eight values (the seven pre-fix data and the one post-fix value) come from the same distribution and are independent, there is only a $1/8$ chance that the eighth value will be the largest. Therefore, we can say with $1 - 1/8 = 87.5$% confidence that the fix has improved the failure times. This procedure also correctly handles the censoring in the last value, which really records a failure time of some unknown value greater than 233. (If a parametric prediction limit happens to exceed 233--and I suspect [based on experience and on the result of @Owe Jessen's bootstrap] it would be close if we were to calculate it with 95% confidence--we would determine that the number 233 is not inconsistent with the other data, but that would leave unanswered the question concerning the true time to failure, for which 233 is only an underestimate.)

Based on @csgillespie's calculations, which--as I argued above--likely overestimate the confidence as $98.3$%, we nevertheless have found a window in which the actual confidence is likely to lie: it's at least $87.5$% and somewhat less than $98.3$% (assuming we have any faith in the geometric distribution model).

I will conclude by sharing my greatest concern: the question as stated could easily be misinterpreted as an appeal to use statistics to make an impression or sanctify a conclusion, rather than provide genuinely useful information about uncertainty. If there are additional reasons to suppose that the fix has worked, then the best course is to invoke them and don't bother with statistics. Make the case on its technical merits. If, on the other hand, there is little assurance that the fix was effective--we just don't know for sure--and the objective here is to decide whether the data warrant proceeding as if it did work, then a prudent decision maker will likely prefer the conservative confidence level afforded by the non-parametric procedure.

Edit

For (hypothetical) data {22, 24, 36, 44, 15, 20, 23} the exponential probability plot is not terrifically non-linear:

(If this looks non-linear to you, generate probability plots for a few hundred realizations of seven draws from an Exponential[25] distribution to see how much they will wiggle by chance alone.)

Therefore with this modified dataset you can feel more comfortable using the equations in Krishnamoorthy et al. (op. cit.) to compute a prediction limit. However, the harmonic mean of 25.08 and relatively small SD (around 10) indicate the prediction limit for any typical confidence level (e.g., 95% or 99%) will be much less than 223. The principle in play here is that one uses statistics for insight and to make difficult decisions. Statistical procedures are of little (additional) help when the results are obvious.

-
+1 very nice answer. –  csgillespie Nov 3 '10 at 16:49
thank you so much for this thoughtful answer. I have included an Update section in my question that attempts to address your concern as well as ask an additional question. –  SiegeX Nov 3 '10 at 17:54
I've added an UPDATE 2 section that attempts to use the equations in the Normal Based Methods in a Gamma Distribution paper by Krishnamoorthy, Mathew, and Mukherjee. Do you agree with my results? If so, I intend to get more failure data by removing the fix then testing the new data set against an exponential probability plot to ensure they fit the exponential distribution. I'll then update my R code with the new dataset to get my new upper limits. –  SiegeX Nov 4 '10 at 1:51
@SiegeX Good job! Using equation (2) of the paper with your data and your confidence levels, but using a different computing platform as a check, I obtain UPLs of 42.5763,50.1065,58.4814,71.5132,83.4132,121.271,143.914,220.331. That looks very close to what you have plotted. –  whuber Nov 4 '10 at 2:13

There are a few ways of doing this problem. The way I would tackle this problem is as follows.

The data you have comes from a geometric distribution. That is, the number of Bernoulli trials before a failure. The geometric distribution has one parameter p, which is the probability of failure at each point. For your data set, we estimate p as follows:

$$\hat p^{-1} = \frac{100 + 22 + 36 + 44 + 89 + 24 + 74}{7} = 55.57$$

So $\hat p = 1/55.57 = 0.018$. From the CDF, the probability of having a run of 223 iterations and observing a failure is:

$$1-(1-\hat p)^{223} = 0.983$$

So the probability of running 223 iterations and not having a failure is

$$1- 0.983 = 0.017$$

So it seems likely (but not overwhelming so) that you have fixed the problem. If you have a run of about 300 iterations than the probability goes down to 0.004

Some notes

1. A bernoulli trial is just tossing a coin, i.e. there are only two outcomes.
2. The geometric distribution is usually phrased in terms of success (rather than failure). For you a success is when the machine breaks!
-
Surely p^ is 7/100 (or however many trails were used in the original test)? –  James Nov 3 '10 at 16:09
No I don't think so. I took the question to mean that he ran device until he reached a failure. So it failed on 100 iterations, 22 iterations, ... –  csgillespie Nov 3 '10 at 16:14
Then shouldn't each test have its own failure probability? –  James Nov 3 '10 at 16:56
Each test does have it's own failure probability - it's 1. We've got the data to prove it ;) Anyway I'm assuming that each test is a replicate from the Geometric distribution. –  csgillespie Nov 3 '10 at 16:59
Assuming a geometric distribution is assuming that the failures are independent. In other words, on any trial the device fails with probability p, independent of the other trials (like tossing a coin). This is a reasonable starting point, but if the failure is due to something like a memory leak, then the longer the device runs, the more likely it is to fail, and the assumption of independent failures would be invalid. –  PeterR Nov 3 '10 at 19:06

I think you could torture your data a bit with bootstrapping. Following cgillspies calculations with the geometric distribution, I played around a bit and came up with the following R-code - any corrections greatly appreciated:

fails <- c(100, 22, 36, 44, 89, 24, 74) # Observed data
N <- 100000 # Number of replications
Ncol <- length(fails) # Number of columns in the data-matrix
boot.m <- matrix(sample(fails,N*Ncol,replac=T),ncol=Ncol) # The bootstrap data matrix
# it draws a vector of Ncol results from the original data, and replicates this N-times
p.hat <- function(x){p.hat = 1/(sum(x)/length(x))} # Function to calculate the
# probability of failure
p.vec <- apply(boot.m,1,p.hat) # calculates the probabilities for each of the
# replications
quant.p <- quantile(p.vec,probs=0.01) # calculates the 1%-quantile of the probs.
hist(p.vec) # draws a histogram of the probabilities
abline(v=quant.p,col="red") # adds a line where quant.p is
no.fail <- 223 # Repetitions without a fail after the repair
(prob.fail <- 1 - pgeom(no.fail,prob=quant.p)) # Prob of no fail after 223 reps with
# failure prob qant.p


The idea was to get a worst-case value for the probability, and then use it to calculate the probability of observing no fail after 223 iterations, given the prior failure probability. The worst case of course being a low failure probability to begin with, which would raise the likelihood of observing no failure after 223 iterations without fixing the problem. The result was 6.37% - as I understand it, you would have had a 6%-probability of not observing a failure after 223 trials if the problem still exists.

Of course, you could generate samples of trials and calculate the probability from that:

boot.fails <- rbinom(N,size=no.fail, prob=quant.p) # repeats draws with succes-rate
# quant.prob N times.
mean(boot.fails==0) # Ratio of no successes


with the result of 6.51%.

-
What is this doing exactly? Im not familiar with R-code but if I had to guess it is using the original 7 figures to extrapolate out the exponential series with 100,000 data points? –  SiegeX Nov 3 '10 at 18:26
I edited the code with comments (which in R are after a #) –  Owe Jessen Nov 3 '10 at 18:59
PS: Of course, this is not really the worst case (which would be a draw of 7x100), but then you wouldn't have to use bootstrap, because you would ignore the other observations. –  Owe Jessen Nov 3 '10 at 19:15
Thank you for the comments, I understand it much better now. However, could you put into layman's term what "quantile(p.vec,probs=0.01) #calculates the 1%-quantile of the probs" is essentially telling us? I see that we get a value of 0.0122 for the 1% quantile. How do I interpret that? –  SiegeX Nov 3 '10 at 22:46
The quantile is a rank-based value. One you probably know is the median: You sort your observations according to size, and the median is the observation in the middle (say, you have the income of 100 people, then the median income is the 50th-highest income). Now generalizing, the 1%-quantile is the probability at 1% of the observations. The advantage of a quantile is that it is non-parametric, you don't have to make assumptions about the underlying distribution, and that it is less biased than the quantile-function of a symmetrical distribution when you have skewed data. –  Owe Jessen Nov 4 '10 at 9:16