Quick Takes - Denver Edition

Snow!

Watch This to be Lulled to Sleep by the Dulcet Sounds of my Voice

...as I talk about one of the most boring subjects known to man.

For those of you with persistent insomnia, a macabre fascination with what I've been doing at work over the last few years, or my adoring fan who just have to know what I'm up to (hi mom!), read on.

Firstly, point your web browsing thingamabobber to (I'll give this in the phonetic): India Sierra Oscar Golf Romeo Alpha Papa Hotel (dash) Sierra Oscar Foxtrot Tango Whiskey Alpha Romeo Echo (dot com). If you entered it correctly, you will be greeted with a smattering of words—or as I like to call it, a smwordening. You may close the smwordening (or if you're really interested, you can click to learn more).

No, you will not understand what half these terms mean.

Once you're into the website proper, click the "Downloads" link on the right side of the menu that runs across the top of the page. To continue, you will need to fill out a download form with your name, company, address, and the time that you leave for work. If you don't want to fill out this form with your own personal information, just do what I do:

Ask for Jenny.

Or you can just put that you're Joe's mom and you want to see what he's been hard at work on.

Where it asks for your email address, be sure to use the email address of some poor sap who you don't know, but who has a convenient-to-type address. Or just something random and funny. I recommend using joecool(at)gmail(dot)com.

Once inside the download area, select the first "Download" button.

Then wait. And wait. (Our FTP server is kinda slow.)

After the hamsters have personally delivered the requisite bits to your ISP, you can run the installer. Psych! No, you can't. You need the password to install! Alas and alack, that is not something I can post here, even in the phonetic. To get the password for the demo installer, you will need to contact me via email. I can be reached at something_random_and_funny@gmail.com, or some.poor.sap@yahoo.com.

Once I've given you the password—which isn't "password", "swordfish", or "xyzzy"—you can install Reliability Workbench. It requires the .NET framework 4.0, which will be installed if you don't already have it. When it's done installing, run the program. The "demonstration mode" dialog will appear. Click "Continue in demonstration mode".

OK, now you're ready to see what I've been doing for the last few centimeters of hairline. If you've followed me up to this point, you should see something that looks like this:

Rollover to see the UI terminology.

Choose "Reliability Growth" from the module selector.

The Reliability Growth module is my baby. Even before Damien was born, I was hard at work nurturing it and developing it into something worthwhile. My boss basically said. "We need a Reliability Growth module. Here's a document that explains what that is. Here's the framework for the program. Go." I did all the coding, UI design, testing, and documentation.

But what is Reliability Growth, you ask? Excellent question! Reliability Growth is an analysis method used to determine if the reliability of a system is increasing, decreasing, or staying about the same over time. It can be used to examine the effectiveness of design modifications or changes in a maintenance policy. For instance, by recommending oil changes every 3000 miles vs. every 5000 miles, does the engine reliability of a fleet of vehicles increase? By using a different weight of oil, is their an appreciable impact on performance?

There are three basic Reliability Growth methods that our software uses: Power Law, Power Law for grouped failures, and the Crow Discrete method.

The first method, the Power Law model, is suitable for continuously-operating systems with immediately-revealed failures. Data consists of the time to each failure. The software will analyze the times to failure to look for an increasing or decreasing trend. It also provides a statistical test for growth, as well as a goodness-of-fit test that verifies the predicted growth curve matches the data.

To begin, select the Add, Growth Set menu option, or corresponding toolbar button. The Growth Set Properties dialog will appear. Choose the "Items" tab and enter the following data:

Click OK when done.

Let's pretend that we tracked the failures of our disgronificator over time. We want to know if our new maintenance policy is having a statistically significant impact on its reliability. We recorded the time at which each failure occurred, and so we enter them into the Growth module.

If it is not already shown, select "Plot" from the Right Window Mode selector. The program will show a plot of the cumulative number of failures vs. time.

If it is not already selected, choose the "Power Law" option from the Method selection drop down list on the toolbar, and select the Analysis, Calculate Method Parameters Automatically menu option or corresponding toolbar button (big red Eta on the toolbar).

The program will automatically calculate values λ and β such that:

E[N(T)] = λT^β

Where E[N(T)] is the expected number of failures to have elapsed by time T. These parameters describe the best-fit growth curve, and are displayed to the right of the plot. Also, the curve described by λ and β is displayed on the plot.

λ is the scale parameter. Ignore it. β is the one we want. This is the shape parameters. β > 1 indicates negative reliability growth (or the reliability is getting worse), β = 1 means zero growth (neither improving nor worsening reliability), β < 1 indicates positive growth.

Next is the U statistic. I won't go into details on how it is calculated, but you can look in the Help file if you're interested. Under the hypothesis that there is zero growth, U will be a standard normally-distributed random variable. By comparing U to the critical values of the standard normal distribution, the Growth module can determine if there is evidence of positive or negative growth. A two-sided test at the α significance level has critical values u_1–α/2 and –u_1–α/2. If U falls between these two values, we can say that U is a standard normal random variable at the α significance level, and that there is no evidence of reliability growth. If U < –u_1-α/2 or U > u_1-α/2 then there is evidence at the α significance level for positive or negative reliability growth, respectively.

Don't pretend you understood any of that, unless your screen name starts with "M" and ends in "r. Pickles".

You can change the significance level in the Growth Set Properties dialog. Double-click "GS1" in the project tree on the left; select "GS1" and choose the Edit, Properties menu option or corresponding toolbar button; or select "GS1", right-click it, and choose "Edit Properties" from the context menu. Either way. I'm not picky.

You can choose the significance level on the Analysis Options tab.

The program will show the calculated U statistic, the standard normal critical value for the chosen significance level, and a line indicating whether there is positive, negative, or no growth.

The program will also perform a goodness-of-fit test to determine if the power law model adequately fits the supplied data. This is done by calculating a Cramér-von Mises statistic and comparing it to the Cramér-von Mises critical value for the 10% significance level. If the statistic is less than the critical value, then the test is passed and the curve fits the data.

I don't pretend to understand how that one works.

Lastly, the program will also calculate the 90% upper and lower confidence bounds for β, as in "I'm 90% confident that I don't understand a lick of this, more or less."

If you want to see the reliability trend, you can also view graphs of the failure intensity and the MTTF (mean time to failure) over time. Select one of those plot types from the Plot Type selector above the graph.

Plot Type selector

(You might wanna go back to the Cumulative Number of Failures plot when done.)

Next up is the Power Law model for grouped failures, also known as case two of the power law. Like the regular Power Law model, this is used for continuously-operating systems, but differs in that the failure times are not known exactly. Instead, data consists of discrete time intervals, and at the end of each interval is a known number of failures.

Add a new Growth Set (Add, Growth Set menu item, if you've forgotten), and enter the following data:

Click OK when done. Pretend this data represents the hyper-electronic transducent geoplasmic nodal collectimizer component of our disgronificators. If that component fails, we don't know about it until we open up the disgronificator and check. Every 200 hours we check all our disgronificators and see how many have failed. We make some repairs and continue. We want to know if the changes to our repair process are improving reliability, so we enter the number of failures that each inspection detects, and let the software figure it out for us.

Select the "Power Law, Grouped Failures" option from the Method selector drop down list. As before, the program will calculate λ and β, test for growth and goodness of fit, and give upper and lower confidence bounds. However, here the test for growth and goodness of fit are both based on  the chi-squared distribution, which in itself was a bit of a challenge. It was tricky finding an algorithm that would calculate chi-squared critical values for an arbitrary number of degrees of freedom. In the end, I found an old open-source math library for C that I adapted. Nothin' like debugging 20-year-old uncommented code written by Some Guy on the Internet.


The last method is an adaptation of the power law for discrete systems, known as the Crow Discrete model, for its developer, Dr. Larry H. Crow. Now, let me tell you something about this method. I curse the day Dr. Crow was born. The grief he and his discrete power law method has caused me is greater than that of Andromache felt for her beloved Hector. Allow me to explain.

Like the Power Law model, the Crow Discrete method finds λ and β such that E[N(T)] = λT^β. However, this method is suitable for one-shot systems, such as a rocket launch, which are either successful or failures. Testing takes place over successive stages. Each stage consists of several trials, each of which is either successful or not. After each stage, a design change is made that will impact reliability. We want to see if there's an increasing trend in reliability.

No problem! says Dr. Crow. That's easy! The answer will be the values of λ and β that satisfy these equations:

Yeah, have fun with that.

Where k is the number of stages,
t is the cumulative number of trials through stage i,
n_i is the number of trials in stage i, and
m_i is the number of failures in stage i.

"AAARRRGGGH!" says I. Of course, nowhere did Dr. Crow give any advice on how to solve these equations.

To put this problem graphically, if you plot the sum of the absolute values of the two equations for a range of λ and β, and invert the graph, you get this:

λ = 0.02 to 2.0, β = 0.01 to 1.0

Your mission is to algorithmically find the highest peak of that graph. Oh, but it's not even so simple as try a range of λ and β, find the highest peak and iterate towards the solution from there. Oh, no no no. Those peaks have fractal-like complexity. It's most likely you'll just hit the side of a mountain and discount it as having too low a value.

λ = 0.565 to 0.664, β = 0.73 to 0.83

My first attempt at solving this involved sitting in the corner of the room and crying. That didn't work. My second attempt involved sacrificing a goat to Baal. Didn't work either.

After Googling around, I stumbled upon a paper that suggested a genetic algorithm approach. I followed that, implemented it, and hey! it worked! Kind of. Sort of. Well, so long as I didn't blow on it too hard, it sometimes gave me the correct answer for one particular sample data set. But most of the rest of the time it was hilariously bad. And slow.

But it was cool! It taught me evolutionary programming and genetic algorithms and awesome buzzwords. I was so sad when I had to remove it from my code. I left a haiku in the comments:

A work of genius,
This code was simply awesome,
But superseded.


Anyways, I stumbled along with that method for a while, right up until final beta testing when I realized I would be ashamed if we charged customers for that module. So, with the clock running down until final release, I turned back to the internets.

After some searching, I found an obscure DOD document site that sold a paper that, while first researching the problem, I had only heard of in songs and stories. The paper is called "Maximum Likelihood Estimates for the Discrete Application of the AMSAA (Army Material Systems Analysis Activity) Growth Model", by William P. Clay, written in 1987. (AMSAA is the organization Dr. Crow was working for when he invented the Crow Discrete reliability growth model.) I bless the day Mr. Clay was born.

You see, William Clay not only described an algorithm to solve the Crow model (I didn't completely understand it; my eyes glazed over when he started talking about iteratively solving the double-partial derivatives to find the maximum likelihood estimators), but he provided the code for a computer program written in FORTRAN 77! Let me tell you, never have I been so happy to have taken a semester of FORTRAN in college. "Never use it," my left kidney.

So anyways, after converting a FORTRAN 77 program to C#, for which I think I insulted my compiler for not being object-oriented enough, I got it to run and solve the problem.

So you can appreciate all my hard work, add a new Growth Set with this data:

So suppose we were testing nuclear bombs by dropping them on Moscow. In the first 14 attempts, 5 failed to explode. We made a design change and then tried again. In the next 19 attempts (trials 15-33), 3 failed to explode. After another change, we made 15 more trials (34-48) and had 4 failures. After the last change, we had 20 trials (49-68) and 4 failures.

Click OK, and choose "Crow Discrete" from the Growth method selector.

We now see that we are more reliably nuking Russia.

Don't say I never do anything for you.

Answers to the Programming Quiz:

Tom got the answer to number one. The assignment

(-1 ^ (j - 1))

will always be -1 for any integer value of j, because of Visual Basic's order of operations. The carat operator is applied before the minus operator. This added about 15 minutes to my debugging time.

As for the second question, Dave was right on with his assessment:

"it seems to take a set of real numbers and multiplies them against each other, alternately adding/subtracting the products."

That is, in fact, exactly what it does, and had he known the law in question, you probably would have gotten it from here. The law (or maybe formula would be a better term) is the general form of the Addition Law for n Independent, Non-exclusive Events, and is used in the field of probability to determine the chances that at least one of n events occurs, given the probabilities of each individual event.

To explain this, let me first define independent and non-exclusive. Two events are independent if the occurrence of one does not affect the occurrence of the other. Example: a die rolling 6 and a coin landing heads.

Two events are exclusive if they cannot occur at the same time. Example: a coin landing heads, and the same coin landing tails.

So the probability of either of two independent, non-exclusive events occurring is given by the equation:

P(A+B) = P(A) + P(B) - P(A)∙P(B)

Where P(A+B) is the probability that either A or B occurs, P(A) is the probability that A occurs independently, and P(B) is the probability that B occurs independently. (I am using + notation because I can't find the stupid Union symbol in the character map.)

This can be easier to understand with a Venn diagram.

If we want to find the area covered by both ellipses (the probability that at least one of the events occurs), we start by adding the area of each ellipse (P(A)+P(B)). But then we have added the intersection of the two events twice, so we need to subtract it out once (–P(A)∙P(B)).

We can determine the probability that any of three independent, non-exclusive events occurs with the equation:

P(A+B+C) = P(A) + P(B) + P(C) - P(A)∙P(B) - P(A)∙P(C) - P(B)∙P(C) + P(A)∙P(B)∙P(C)

To understand this, let's look at another Venn diagram.

To find the area covered by all three circles, we begin by adding in the area of each circle (the independent probability of each event, P(A) + P(B) + P(C)), but then we have added some regions more than once. To compensate, we can subtract out the intersection of every two circles (–P(A)∙P(B) – P(A)∙P(C) – P(B)∙P(C)). But now we have completely subtracted out the center region, the intersection of all three circles (P(A)∙P(B)∙P(C)).

I think you see where this is going for four events, etc. It gets very complex. In fact, the general form of the Addition Law, for n independent, non-exclusive events is:

(I)

In essence, we add in the independent probability of each event, subtract the products of each combination of two events, add the products of each combination of three events, subtract the products of each combination of four, etc., until we finally add or subtract (depending on if there are an odd or even number of events) the product of every event probability. This is what my VB code was attempting.

(There is an alternate, less complex way of attempting this. Instead of a brute force method, we can instead consider the inverse probabilities. Instead of finding the area covered by the circles/ellipses, we can find the area not covered by them, and subtract that from 1. To wit:

De Morgan's theorem states: P(A+B)' = P(A'∙B') = P(A')∙P(B')

That is, the probability that neither A nor B occurs is the probability that A does not occur and B does not occur. So if we subtract that from 1, we get:

P(A+B) = 1 – P(A')∙P(B') = 1 – (1 – P(A))∙(1 – P(B))

Which, for the general form, expands to:

(II)

)

So what brought this all up?

Probability theory figures heavily in fault tree analysis. A fault tree is, essentially, a collection of events (each with an associated probability of occurrence, specifically a probability of failure) connected via logic gates, such as AND and OR gates (and sometimes VOTE gates).

What? Explanation of the logic gates? Ok, we can do that.

This is an AND gate. It denotes that the output occurs if all of the input events occur.

An OR gate denotes that the output occurs if any of the inputs occurs.

A VOTE gate denotes that the output occurs if at least n of the inputs occurs (in this case 2).

(These symbols may seem familiar to those of you who took computer architecture. The symbols are the same ones used in logic gate diagrams, but turned sideways.)

There are also other miscellaneous gates, such as a NOT gate or XOR gate, but that's a bit beyond the scope of a simple introduction.

So the basic premise of fault tree analysis is we want to know about how often (the probability) a TOP event will occur, based on the occurrence of basic events. The TOP event is usually some hazard that we would like to prevent. A basic event is a failure or occurrence which may cause the TOP event to occur.

So, to create a fault tree, we connect these basic events to the TOP event through intermediate logic gates. These logic gates represent how the failures combine to cause a TOP event occurrence. A small fault tree might look like this:

The way to read this is by looking at the logic symbol of each gate to determine which of its input events must occur to cause the gate to occur. For example, if EVENT5 and EVENT6 occur, this will cause GATE3 to occur, because it's an AND gate and requires the occurrence of all input events. If GATE3 occurs, GATE2 will also occur, because it's an OR gate and the occurrence of any of its inputs will cause its occurrence. Now, if any two of the inputs to GATE1 occurs, say EVENT1 and EVENT3, then GATE1 will occur (because it's a VOTE gate). And if both GATE1 and GATE2 occur, then TOP1 will also occur. So we could say that TOP1 will occur if EVENT1, EVENT3, EVENT5, and EVENT6 all occur at the same time. We call this occurrence of events a minimal cut set. (I might talk more about cut sets in the future, if you wish, but as it is, this post is running too long.) Another minimal cut set might be EVENT2.EVENT3.EVENT4. (Note, EVENT2.EVENT3.EVENT4.EVENT5 is also a cut set, but not minimal. Do you see the difference?)

I will skip over how to produce the full list of minimal cut sets for this tree, and just list them all. They are: 

EVENT1.EVENT2.EVENT5.EVENT6

EVENT1.EVENT2.EVENT4

EVENT1.EVENT3.EVENT5.EVENT6

EVENT1.EVENT3.EVENT4

EVENT2.EVENT3.EVENT5.EVENT6

EVENT2.EVENT3.EVENT4

By this point, you can probably see how to produce them yourself. Now each cut set can be thought of as an event, itself. So, basically, what we have here is a list of non-exclusive, independent events, the occurrence of any of which will cause the TOP event.

Sound familiar?

The obvious way to solve this, then, would be to use... the Addition Law!

And this is exactly what the software (that company I work for) writes does. You tell it the events, how they're logically connected, and some quantitative data, such as a failure rate or MTTF (mean time to failure), and it will calculate:

1) The probability of each event occurring,

2) The probability of each cut set occurring (simple: the product of the probabilities of each event in the cut set)

3) The probability of occurrence of the TOP gate (using the Addition Law).

Note how quickly the addition law gets complex. This small and simple tree produces six cut sets, which would require 6 + 15 + 20 + 15 + 6 + 1 = 63 product terms in the addition law. Just the other day I was working with a fault tree that had 152 cut sets for one intermediate-level gate. This is why our software will use approximation methods and all sorts of other tricks to be able to solve the tree. In fact, one of the approximation methods involves using de Morgan's theorem and equation II. This is called the Esary-Proschan approximation method, but you didn't really want to know that.

If you would like, I could talk a little bit more about fault tree analysis later, and give real examples, instead of odd abstractions. But this post has gone on too long.

And now you know what I do for a living!

So back to the topic at hand, why was I writing a VB program to implement the Addition Law? Basically, I needed to verify the results of a fault tree. By hand.