Random Events, Trends, and Hot Temperatures

Random Events, Trends, and Hot Temperatures

(c) Ken Osborn Jan 17, 2015

Today's papers (Jan 17, 2015) say 2014 was the hottest on record since 1880. An engineering friend says the trend may be just a random thing. Could that be? Though the data do show a reasonable closeness to an independently produced set of random values with the same average and variability (standard deviation), I decided to test it by comparing the first half of 1880-2014 to the second half.

If random, dividing the data into two halves and ordering the values from lowest value to highest within each half would produce two superimposed sets that would be indistinguishable. If they are part of a trend of increasing temperatures, the sets would not superimpose and the second set would be shifted to the right on the chart. So which is it? Check out the graphs.

 [Data source: http://data.giss.nasa.gov/]


First, let me state that I do understand that random events can generate what appears to be a trend, as in this example:


Chart 1: Example of a 'trend' from a random generation of values 


The Random Walk chart displays what appears to be a strong trend with a correlation coefficient of 0.96, sufficient to earn an award for superior performance by an inebriated soul staggering along a mostly straight path.  I should note that it took over 50 runs of the program to get this result and doubt I could reproduce this particular set again.  


But if we can get a randomly generated collection of results to look like a trend, how do we know that our trend was not generated randomly?  

One way, of course, would be to rerun the events and see if we get the same trend.  However, since I'm looking at data collected from 1880 to 2014 that is not a realistic approach.  Another way would be to take the data, split it into two halves, and compare the two halves.  If the data are a random distribution around a central value, the first half of the data should match the second half of the data.  But first, let's try that with some randomly generated data.


Chart 2: Plot of 100 random values generated by a Monte Carlo simulation 




In chart 2, 100 random values ordered from low to high are plotted against their rank (probability). Half of the values exceed and half are below the average of zero and are symmetrically distributed around the center.  In other words, these data exhibit a nice Normal type distribution.  Note that the average and standard deviation for these data are the same as the temperature data to follow.  




So now if those data are split into two halves and each half is ordered independently from lowest value to highest value and plotted against its rank, what do we get?  


Chart 3:  Comparison of two halves of a randomly generated set of values after ranking each half




The values in blue represent the first half of 100 random values and the ones in red the second half.  Each half was then ranked from lowest to highest value and plotted against the ranking (1 to 50).  The superimposition is not exact, but these are randomly generated results so one should not expect an exact agreement between individual values.  But maybe a real trend would also show something like this.  Let's try it.  


Chart 4: A trend (Y= 10X+5) of 100 values partitioned into two ordered halves 




The trended data set of 100 values was generated from the formula Y = 10X+ 5.  The set was divided into two halves and each half was ordered from lowest value to highest value then plotted against its rank from 1 to 50.  Unlike in chart 3 with conformable data sets, these two sets show no overlap at all.  So how will this work with real data?  


 

Chart 5: Comparing the ranked temperature anomalies from 1880 to 2014 with a random data set  (Source: http://data.giss.nasa.gov/)

The values in green in chart 5 represent the temperature anomalies from 1880 to 2014 ranked from lowest to highest value.  Each value represents the deviation from the average for the 20th century.  The values in blue were randomly generated using the mean and standard deviation from the temperature anomaly set.  They do look as close as two separate runs of a random number generator.  But remember, the real test is to see if the first half of the data (1880 to 1946) matches the second half (1947 to 2014).  Any guesses?  

Chart 6: Comparison of two halves of the 1880-2014 temperature data anomalies


In chart 6, the values in red are for the years 1880 to 1946 and the values in green for 1947 to 2014.  Each set is ordered from its lowest to highest value and plotted against corresponding year.  They do not match and are clearly two separate distributions.  I leave the conclusion to you as to whether these data have been generated by random events.    

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.

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!

The Drunkard's Walk

Random events can be misleading.  Flip a coin.  It has a 50% chance of landing heads and 50% tails.  Now suppose it has landed heads 10 times in a row.  What are the chances it will land heads up on the next toss?  If it is truly an unbiased coin, the odds for a heads up will be 50%.  Past events have no effect on future events.  But what are the odds that a coin will land heads up ten times in a row? 

The chart below shows a plot generated from a random number generation program, Drunkard's Walk.  Try running the program Drunkard's Walk yourself and see what kind of trends you can generate from a random process. Enter any number in cell B8 (Count).  Press return.  Each plotted point is the sum of the previous plotted point plus the random assignment of a one or a minus one.  Are there trends?  Does it always look random?