Using a Monte Carlo Program to Define Process Control Parameters

Using a Monte Carlo Program to Define Process Control Parameters

© Ken Osborn 2013

Abstract: Profiler is a program written in Excel® Visual Basic® and uses a Monte-Carlo simulation to fit process data to a Normal distribution. Output includes raw data and simulation statistics and outliers are identified and flagged. If the data are not skewed and show a good fit except for the outliers, process control limits can be calculated from the Normal curve average and standard deviation.

The Problem
Controlling any production process like weight of material added to a mix, treatment and discharge of a waste product, or analysis of a material requires that the variance be predictable. When a process exceeds the established variance it is out-of-control. Process control parameters include the process average, standard deviation, and control limits. Control limits are frequently set at the average plus or minus three standard deviations.

The difficulty with this approach is that if the statistics used to set the control limits are generated when the process is out-of-control they will overestimate the acceptable range of variance by including outliers that may not be initially suspect. This may be particularly true at process startup when there is no baseline of data.

The Solution
To avoid the problems of outliers driving the process control statistics, the control limit statistics can be derived from a Normal distribution fitted to the process data. The Normal distribution will exclude outliers and more closely reflect the variance in the process output when it is under control.

Profiler, a program written in Excel, fits a Normal curve to a data set using a Monte Carlo simulation. The fit is iterative starting with an initial guess which checks for fit against the 25th and 75th percentiles of the data. Adjustments are made to the mean and standard deviation of the simulation until a pre-set error is achieved based on the 75th percentile. The error is optimized for low error vs time for convergence. Output consists of 100 data points, simulation mean, standard deviation, 25th, 75th, and 99th percentiles, and flags all values in the original data set exceeding one standard deviation of the 99th percentile.

An Example
Data plotted in red in Chart 1 are from a waste discharge stream controlling for mercury. Clearly some of the values are outliers and do not represent the process when it is in control.


Chart 1: Mercury in discharge vs Monte Carlo simulation

The blue line represents the Monte Carlo fit to a Normal distribution. The majority of the mercury discharge data exhibit a close fit to the simulation curve with outliers clearly visible. The statistics for the mercury and simulation results (Table 1) show the divergence from normalcy at the extreme ends of the distribution. Of course, we are only concerned about the high outliers that would drive the statistics if retained. Using the raw data, a 3-sigma upper control limit (3S UCL) calculates as 240 while the simulation UCL is 141.

Since the raw data clearly include results from the process when it was out of control, outliers should be excluded. But which ones? Using the raw data 3S UCL leaves several high values that could misrepresent the process when in control. Conversely, using the simulation UCL may remove more values than justified. Using an alternative determination of 1 standard deviation above the 99%tile (UCL_PTile), the simulation UCL_PTile is 133. If the raw data 3S UCL is used, there are six outliers (240, 250, 290, 300,310, and 400). If the simulation UCL_PTile value is used, there are ten outliers (175, 180, 190, 200, 240, 250, 290, 300, 310, 400). These are the values plotted in red at the extreme right in Chart 1.


Table 1: Comparison of Monte Carlo and Raw Data Statistics

When the 3S UCL is recalculated for the raw data excluding the outliers exceeding the simulation UCL_Ptile, the statistics for the raw data are comparable to the simulation (Table 2). In addition, whether the 3S UCL or the UCL_PTile is used makes little difference. Once the outliers have been identified, the data represent a process that is predictable.



Table 2: Process statistics after removal of outliers

Conclusion

Statistical control of a process requires establishing realistic control limits based on data that represent the process in control. At startup the process may not be in control and include outliers. Fitting the data to a Normal curve from which the control statistics are derived provides a simple mechanism to identify the outliers and establish statistical process control limits. A Monte Carlo simulation can be used to fit the data to a Normal curve.

If you would like a copy of the 'Profiler' program, please write me. You may use the program for your own purposes or post on your blog with credit. You may not use the program for commercial purposes. Write me at kozborn@sbcglobal.net.

Using a Monte-Carlo Simulation to Model Natural Selection

Using a Monte-Carlo Simulation to Model Natural Selection

Ken Osborn © 2012

Darwin postulated that natural selection was the driver of the evolution of species diversity[1].   Elements of this model include 1) a selective agent as a reproductive threshold, 2) genetic variability (plasticity), and 3) a continuing process of genetic variation with each succeeding population.   In the Darwinian model of evolution, natural selection provides a reproductive threshold.  Individuals within a population whose genetics allow them to leave more offspring when challenged by the threshold determine the genetic direction of the population and the process leads to a “more fit” population, where fitness is defined in terms of the population’s ability to adapt to the threshold.  If the selective agent threshold is too high, the population may become extinct. 

For example, if color variation within a population of moths allows some moths to escape predator detection then predation pressure would act as a selective agent that would result in a change in distribution and abundance of color expression within the moth population.  Of course, the difficulty of finding camouflaged moths could act on the predator population as a selective agent driving changes that improved their abilities to detect prey.  Another example would be the effect of climatic variation on a population of Pika.  If genetic variation within the Pika population was insufficient to provide some individuals who could survive and reproduce in the presence of shifting summer temperatures, the Pika population would become extinct. 

Mathematically these elements can be modeled.  I have written a program in Excel using a Monte-Carlo process to generate a population of values that can evolve into another population of values in a “natural selection” process.  In the model, two numbers define a starting population: 1) the population average (mean) which represents Darwinian fitness, and 2) population relative standard deviation which represents genetic plasticity.  Once a threshold is set, the program generates a new population of values using the mean and relative standard deviation (RSD).  The Monte-Carlo equations generate a Normal distribution of values with mean = 0 and standard deviation = 1.  These values are converted using the population mean (fitness) and RSD (genetic plasticity).  Population values that exceed the threshold become the basis for the next population.  The process can be repeated with or without a new threshold. 

There are several observations I have made from this model.  Some seem reasonable and others are counter-intuitive.  All are potential candidates for testing, though some might be trivial and others cost-prohibitive.  I leave the testing for others to pursue. 

Figure 1 is an example of the program output.  The distribution of values will be referred to as a population.  Fitness is the average value of the population.  Spread is the standard deviation.  Plasticity is the relative standard deviation (RSD).   Adaptation is the percentage of a distribution that exceeds a threshold.  Options allow for specifying an initial and challenge threshold, where the challenge threshold is set at a level just below extinction.  The extinction threshold is one that exceeds the highest value in a population.  Survivors refer to population values that exceed a threshold.

The blue curve is for the initial population, Generation 1.  The vertical axis represents values of the population and the horizontal axis is the count of values (i.e., there are 100 individuals in the population).   The values have been sorted to provide a less complex visualization of the population distributions.  Individual values in Generation1 vary from a low of 159 to a high of 257 with a fitness mean of 203 and a plasticity RSD of 10.  Values of Generation 1 that exceed the initial threshold of 230 are used to create Generation 2.  A challenge threshold (250) was used to create Generation 4 from Generation 3 (not shown).  The table above the plot summarizes the population statistics for each succeeding generation. 

Figure 1: Monte-Carlo Evolution Program Output Example  (Run 1)

The most obvious observation is that there is variability around the mean fitness value.  Generation 1 ranges from a minimum of 159 to a maximum of 257.  Note that the relative range between maximum and minimum values decreases with each succeeding generation.  I will discuss this more.  That there is variability within a population is a trivial observation given that these numbers are generated using a random Monte-Carlo process, but it is worth noting. 

The second observation is that for a given set of initial fitness, genetic plasticity, and selective threshold the final outcome is only generally predictable.  Each run of the program will generate a slightly different outcome.  Figure 2 shows the second run with the same starting inputs for mean, plasticity, and threshold values as the first run shown in Figure 1.  The changes are small with a change in mean fitness for Generation 1 from 203 to 204, 241 to 240 for Generation 2, and 255 to 254 for Generation 4.

Figure 2: Run #2 for Initial Settings of Mean =200, RSD = 10, and Threshold = 230. 

Again, because this is a Monte-Carlo process this is a trivial observation.  The driver of genetic plasticity is genetic mutation, a constraint driven but random processes. A biological population with a specified fitness level and genetic plasticity would not respond exactly the same way to each instance of a natural selection challenge even if it were the exactly the same repeated challenge.  If this conclusion is generally valid, the implication is that smaller population size would lead to lower adaptability. 

The third observation is that for a given fitness and plasticity, there is a maximum threshold.  Setting the threshold above the maximum causes the program to crash.  I refer to this as the extinction threshold.  Because it is a Monte-Carlo process, there is variability in the extinction threshold.   A setting that causes a population crash in one run may not result in a crash in the next run.  In fact, the challenge threshold for Generation 4 of 250 was initially set high and was expected to result in a program crash (extinction).  It did on the third run, as seen in Figure 3. Generation 3 had no survivors to contribute to Generation 4. 

Figure 3: Setting the threshold above extinction

What happens if the genetic plasticity is increased?  If the initial fitness is fixed, changing the plasticity would be comparable to having two biological populations with the same overall average fitness but with different ranges from maximum to minimum fitness.  Figure 4 shows what happens when the plasticity for the first run is doubled from 10% to 30%.  Figure 4 shows that the fitness for Generations 2 and 4 has increased relative to run 1 even though the thresholds have not changed.  Generation 4, for example, went from a mean fitness of 256 to 285.  The increase for Generation 1 from 199 to 205 is within the run-to-run variability and not significant though the increase of the population maximum for an individual value from 243 to 334 is.  Note that it would not be expected for the fitness of Generation 1 to change with an increase (or decrease) in plasticity as it was not challenged with a threshold. 

Figure 4: Increasing genetic plasticity

At first it would appear that an increase in plasticity may not be of value as it not only leads to an increase in the maximum values in the population (243 to 334 for Generation 1), but also decreases in the minimum values (from 153 to 104 for Generation1).  Figure 5 used the input conditions for Run 3 with an increase in the Challenge Threshold from 250 to 280.

  

Figure 5: Increasing the Challenge Threshold after an increase in plasticity but a fixed initial fitness

The Challenge Threshold is set to affect only Generation 4 and the results are consistent with that condition.  The mean fitness of Generation 2 shows a modest increase, but is within the run-to-run reproducibility.  Generation 4, however, shows a much larger increase from 285 to 300 and both the max and min values have increased.  The significance for a biological population is that even if the population has some individuals that are well below the mean fitness, the population as a whole may be more responsive to environmental challenges. 

In this model of evolution, is genetic plasticity or genetic fitness a better predictor of adaptability?  I have rewritten the model to allow for multiple runs with a single set of input values for fitness, plasticity, and threshold.  The percentage of values exceeding a threshold (PECT) is a measure of adaptability. 

The model was run twenty times for each set of plasticity and fitness.  Column 1 is the Run number, column 2 is population 1 with fitness of 220 and a plasticity of 10, and column 4 is population 2 with fitness of 190 and a plasticity of 40.  The data are presented in Table 2 (note that there is no Table 1).  In this example, population 2 with the higher plasticity has the higher PECT (57.8 vs 11).  

This is somewhat counter-intuitive.  Consider if the players in a football team at school A were very evenly matched and all could bench press 200 plus or minus 10 pounds while players at school B could press 190 plus or minus 40 pounds.  You might think that the team with the well matched players and the higher average bench press capabilities would be more capable of achieving a 250 pound bench press challenge.  But in fact only 11% of team A could do that while 58% of team B could meet that challenge.

This mathematical model of evolution is a simple one and does not capture the complexity of biological evolution.  For example, predator-prey interactions, group selection, and eco-system feedback loops are not part of this model.  Given its simplicity, however, it has provided a surprising number of testable hypotheses:

1. Increased genetic plasticity increases adaptability.
2. There are limits to adaptability for a given plasticity and genetic fitness: a selection threshold that is too high will result in extinction. 
3. In an environment where the threshold is changing, adaptation to a challenge threshold is greater for the population that has been subjected to the higher non-extinction threshold below the challenge threshold.
4. As the selection threshold increases, genetic plasticity of a population decreases.
5. Genetic plasticity can be a better predictor of adaptability than population fitness.

If you would like a copy of the 'Evolution' program, please write me.  You may use the program for your own purposes or post on your blog with credit.  You may not use the program for commercial purposes. Write me at kozborn@sbcglobal.net.  

---

[1] Charles Darwin, "On the Origin of Species by Means of Natural Selection, or the Preservation of Favoured Races in the Struggle for Life," 1859.

Predator vs Prey - A Mathematical Model - Part 2 of a 3 Part Series

In the first installment of my predator-prey model I presented graphs of population growth under three scenarios that included deer and forage but no predators.  Scenario 1 was for a stable environment (i.e, constant carrying capacity): the deer herd rapidly grew from a starting population until it reached the carrying capacity then leveled out.  Scenario 2 was a moderately variable environment and the deer initially grew but the population fluctuated above and below the carrying capacity.  Scenario 3 was for a highly variable environment.  In this last case, the deer herd grew, fluctuated in number, then crashed and died out.  For a more complete review see http://misterkenblog.blogspot.com/2011/12/predator-vs-prey-will-wolves-dominate.html.

What happens when predators are introduced?  Will the deer herd die out sooner or will it stabilize because wolves keep the deer herd in check with the environment?  Can predictions even be made?  Let's see.

Scenario 1: Start with a deer herd well below the carrying capacity, a small number of wolves, and a constant environment.

The starting conditions are an initial deer population (N) of 5000, a growth rate (R) of 0.5(50% increase in deer herd per generation), carrying capacity (K) of 20000, no variance in the carrying capacity (KV=0), 5 wolves, a reproductive capacity of 0.1 (10%), and a predator efficiency (E) of 0.53 (53% - VERY good hunters), and a requirement (S) of 24 deer/wolf/year to sustain wolf pack growth.

Both the deer herd and wolf pack show rapid initial growth, followed by oscillations of population size, and finally a steady state based on the carrying capacity for the deer herd.  The final steady state is well above the starting conditions for both deer and wolves and the final state for the deer is about half of the carrying capacity.

Predation in Stable Environmet

Scenario 2: Add global warming or some other factor to make the environment variable

The variance for the carrying capacity has been increased from 0 to 1000, or a variance factor of 20% (100x1000/5000).

The oscillations in both deer and wolves has increased and no steady state is achieved although after 500 generations it does not appear that neither deer nor wolves are in danger of extinction.

Add Moderate Amount of Environmental Variability

Scenario 3: Scenario 2 with more predation by increasing the wolf pack from an initial 5 to 50

Increasing the initial number of wolves from 5 to 50 does not change the overall dynamics.  It would appear that starting with an initial wolf pack of fewer than what is sustainable has little long term effect.  I leave it for the reader to try other starting wolf pack sizes once I have posted the interactive spreadsheet.

Increase Predation Pressure

Scenario 4: Keep the starting wolf population at 50 and increase the environmental variability

Increasing the variance in carrying capacity to 100% dramatically shifts the oscillations in both the deer herd size and wolf pack numbers. While neither population crashes, they come perilously close.  I leave it to the reader to try more simulations to see if the wolves, or deer, or both go to extinction under these conditions.

Increase Environmental Variability

Scenario 5: Same as scenario 4 with initial carrying capacity cut in half

In this last scenario, the deer herd crashes and the wolves, lacking a food supply, follow.

Reduce Carrying Capacity

Whether the wolves control the deer population or the deer control the wolf population is still an open question, but clearly the environment controls both.  When the deer herd exceeds the carrying capacity of the environment extremes from one year to the next will ultimately result in a population crash.  In the absence of predation, by this model, the deer herd will subsist only if the environment is very stable.  If the environment is not stable (the normal course of events) predation pressure can help stabilize the deer herd by reducing the herd size and the odds that the deer herd will not exceed the carrying capacity are improved.  Of course if the deer herd crashes so will the predators unless they have a reserve food source.  That in fact is the case, but then it becomes a matter of energetics and whether switching to an alternative food supply for the predator is analogous to a drop in carrying capacity for the prey.

 Next month I will post a link to the statistical model so that the reader can try some scenarios and draw their own conclusions.

Strategies for Sampling in a Potato Warehouse

You rent a warehouse in the Mission District of San Francisco for the purpose of storing potatoes.  The storage space is not big with barely enough capacity to hold 5 tons, or 10,000 pounds of potatoes.  The owner of the warehouse is concerned that you might try and exceed that quantity and has placed a penalty clause in the contract if you do.  He also requires an annual inventory of the quantity of potatoes you are storing.  

You hire a high school student to weigh potatoes so you can determine how much is in a warehouse.  To estimate the weight of potatoes in the warehouse without paying the student to weigh each and every potato you've decided that you will have him make no more than 100 weighings, randomly selected.  Once the potatoes are sampled and weighed you will estimate the total weight by multiplying the average weight per gallon of potatoes times the total volume of potatoes in gallons in the warehouse (you know this because you have a calculator and are good at math).

To do this weighing, would it be better to have the student weigh lots of individual potatoes one at a time or weigh composites of several potatoes at a time in fewer collections of a gallon each?  And should you caution the student to take care in weighing?  Does it matter if you are just a little bit over the 10,000 pounds?  Will the weighings more likely show you are over if the weighings are precise?

The answers may not be immediately obvious and depend on at least three variables: 1) the precision of the student in taking the weights and volumes, 2) the difference between the actual 'true' weight of potatoes in the warehouse and the penalty threshold, and 3) whether the actual weight is above or below the threshold.  To test how these three variables can affect the outcome, I have developed A Monte Carlo Model for Sampling Potatoes or Anything Else.  Try it for yourself.

Enter the measurement precision (also known as percent relative standard deviation  or %RSD) in cell B3, the potato penalty threshold in cell B4,  and the expected true weight of potatoes in cell B5.  The outputs are in cells D7 (percentage hits exceeding the penalty threshold for weights of single potatoes) and E7 (percentage hits for collections of potatoes).

The small chart displays overall statistics of average, standard deviation (a measure of variability), RSD (a measure of precision), and number of Hits (number of times the given sampling method exceeds the penalty threshold).  In the example shown here, the averages for both methods exceed 10,000 pounds and the number of hits are 4800 for the single potato approach and 5500 for the collection of potato approach.  The true value was set at 10,000 and precision at 50% (i.e., sloppy technique) for the default settings. 

The long chart displays the first ten samplings of the student under the two different sampling methods and can be used as a check to track his performance.

Try entering different values for RSD, threshold, and true weight and see what you get.  For example, if the student is sloppy does it make any difference if the threshold is above or below the true value?  If the student is very precise can you get closer to the threshold without a penalty?  And don't be timid about the actual weight of potatoes.  You might try putting in more than 10,000 pounds and see if you can avoid the penalty!