A Self-Organizing Model for Logical Programming EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 3, No. 2, 2010, 163-173 ISSN 1307-5543 – www.ejpam.com Email address: norman.j.morin@frb.gov http://www.ejpam.com 163 © 2010 EJPAM All rights reserved. A Self-Organizing Model for Logic Regression Jerry Farlow Department of Mathematics, University of Maine, Orono, Maine __________________________________________________________________________ Abstract: Logic regression, as developed by Ruczinski, Kooperberg, and LeBlanc [1], is a multivariable regression methodology that constructs logical relationships among Boolean predictor variables that best predicts a Boolean dependent variable. More specifically, it finds a regression model of the form 0 1 1 2 2[ | ] m mg E Y b b L b L b L= + + + where both the coefficients 0 1, ,..., mb b b and the logical expressions , 1,...,jL j m= are determined, whereby a logical expression one means logical relationships among the predictor variables, such as 1 2" ,X X are true but not 5 "X , or 3 5 7" , ,X X X are true but not 1X or 2 "X In paper [1] the authors investigate the use a simulated annealing algorithm. In this paper. the Group Method of Data Handling (GMDH) is used. 2000 Mathematics Subject Classifications: 62 Key Words and Phrases: Logic regression, GMDH, algorithm, self-organizing methods. __________________________________________________________________________ 1. Introduction In this paper we study the possibility of using Ivakhenko’s Group Method of Data Handling (GMDH) algorithm for finding the best logical expressions jL among Boolean predictor variables 1 2, , , mX X X which best predicts a dependent Boolean variable Y . By logical expressions jL , we mean Boolean expressions such as 1 2" ,X X are true but not 5 "X , or 3 5 7" , ,X X X are true but not 1X or 2 "X or " if 1X and 2X are true then 5X is true. Before showing how this algorithm is performed, since the GMDH algorithm is not well- known, we summarize Ivakhnenko’s basic algorithm for ordinary regression. J. Farlow / Eur. J. Pure Appl. Math 164 2. Basic Group Method of Data Handling (GMDH) Algorithm 2.1 Basic GMDH Method If one were to anthropomorphize, one might say the GMDH algorithm builds a mathematical model similar to the way biological organisms are created through evolution. That is, starting with a few basic primeval forms (i.e. equations); one grows a new generation of more complex off-springs (equations) and then allows for a survival-of-the- fittest principle to determine which new off-springs survive and which do not. The idea is that new generations of off-springs (equations) are better suited to model the real world than earlier ones. Continuing this process for more generations, one finds a collection of models that hopefully describes the problem at hand. The process is stopped once the model begins to “over-fit” the real world, thus stopping when the model reaches some level of optimal complexity. In 1966, Ukrainian cyberneticist, A.G. Ivakhnenko, discouraged by the fact that many mathematical models require the modeler to know things about the real world that are difficult or impossible to know, produced a heuristic self-organizing model, called the Group Method of Data Handling algorithm. For more information see Farlow, [1] or view one of the many websites which discuss applications of the technique. Ivakhnenko’s website at http://www.gmdh.net/ will give the reader a complete history of the method. The basic GMDH algorithm can be broken into a few distinct steps. Step 1 (constructing new variables 1 2 ( ,2), ,..., C mz z z ) The basic GMDH algorithm begins with regression-type data points of the form 1 2 1 2, , ,..., , , ,...,i i i imy x x x i n= , where the n observations are subdivided into two groups, the first nt observations are called the training observations and the remaining nc n nt= − observations are called the checking observations. See the data in Figure 1. Normally, about half the observations are chosen to be in each group. Y X Training set Y 1x 2x  mx 1y 11x 12x  1mx 2y 21x 22x  2mx      nty ,1ntx ,2ntx  ,nt mx Checking set 1nty + 1,1ntx + 1,2ntx +  1,nt mx +      ny ,1nx ,2nx  nmx Figure 1: Regression-Type Data for the GMDH Algorithm http://www.gmdh.net/� J. Farlow / Eur. J. Pure Appl. Math 165 Now, for each ( , 2) ( 1) / 2C m m m= − pair of distinct variables ,i jx x one finds the least- squares regression polynomial for y of the form 2 2 i j i j i jy A Bx Cx Dx Ex Fx x= + + + + + (1) (i.e. find , , , , ,A B C D E F ) from the observations in the training set These ( 1) / 2m m − regression surfaces are illustrated in Figure 2. Figure 2: Computed Quadratic Regression Surfaces Now evaluate each of the ( 1) / 2m m − regression polynomials at all n data points and store these values (new generation of variables) in the respective columns of a new array, say Z. The evaluation of the first regression polynomial and its values in the first column of Z is illustrated in Figure 3. Figure 3: Evaluating the Quadratic Regression Polynomials J. Farlow / Eur. J. Pure Appl. Math 166 The object is to keep only the best of these new variables and this is where the observations in the checking set come into play. Step 2 (screening out the least effective variables) This step replaces the original variables (columns of X) by those columns of Z that best predict y, based on the checking set observations. That is, for each column j of Z we compute the root mean square (or some measure of association) jr given by Eq. (2). 2 2 1 2 1 ( ) , 1, 2,..., ( , 2) n i ij i nt j n i i nt y z r j C m y = + = + − = = ∑ ∑ (2) and then select those columns of Z that satisfy jr R< , where R is some prescribed number. The number of columns of Z that replace columns of X may be larger or smaller than the number of columns of X, although often one keeps the number of columns of X constant at m. Note that the test of goodness of fit jr was summed over the observation in the checking set. Step 3 (test for optimality) From Step 2 we find the smallest of the 'jr s and call it RMIN. Then, each time one completes a generation or iteration, one plots the value of RMIN on a graph as shown in Figure 4. Figure 4: Determining the Optimal Polynomial J. Farlow / Eur. J. Pure Appl. Math 167 Experiments have shown that RMIN decreases for a few generations (the author’s experience is maybe 3-5 iterations) and then begins to increase, the reason being the model gets better and better but eventually starts to over-fit the data. Hence, the rule is to stop the algorithm when the RMIN curve reaches its minimum, and then select the column with the minimum jr value of the final array Z as the best predictor. When the GMDH algorithm stops, the columns of Z (in particular the column of Z that has the smallest jr value) contains the computed values of a high-order polynomial of the form in Eq (3) 1 1 1 1 1 1 1 m m m m m m i i ij j j ijk i j k i i j i j k y a b x c x x d x x x = = = = = = = + + + +∑ ∑∑ ∑∑∑  (3) known as the Ivakhnenko polynomial. At each iteration the degree of the Ivakhnenko doubles, and for a p-th order regression polynomial the number of terms in the polynomial will be ( 1)( 2) ( ) / !m m m p m+ + + . For example, if one starts with 10m = input variables 1 2 10, ,...,x x x and the algorithm is continued for 8 generations, the Ivakhnenko polynomial would be a polynomial in 1 2, ,..., mx x x of degree 82 256= . A sample term might involve the variables 2 4 9 11 9 3 1 3 4 6 7 10x x x x x x . Step 4 (Applying the results of the GMDH Algorithm) One doesn’t actually compute the coefficients in the Ivakhenko polynomial, but saves the regression coefficients A,B,C,D,E,F at each generation. Hence, to evaluate the Ivakhnenko polynomial one simply carries out repeated compositions of these quadratic expressions. Figure 5 illustrates this process. Figure 5: Evaluation of the Ivankhenko Polynomial J. Farlow / Eur. J. Pure Appl. Math 168 3. GMDH Algorithm Applied to Logic Regression Step 1: (Divide Observations into Training and Checking Sets) Starting with n observations of m Boolean predictor variables 1 2, ,..., mX X X and a dependent Boolean variable Y, we subdivide the observations into nt training set observations and nc n nt= − checking set observations. For each of the ( 1) 2 2 m m m  − =    distinct pairs { }, : 1,.. 1, 1,i jX X i m j i m= − = + of predictor variables, we find the logical function that best predicts the dependent variable Y from among the 16 binary functions in Table 1. Table 1: Sixteen Boolean Functions 0(0000) FFFF never true ----- 1(0001) FFFT not ( 1X or 2X ) 1 2X X∧ 2(0010) FFTF 2X but not 1X 1 2X X∧ 3(0011) FFTT not 1X 1X 4(0100) FTFF 1X but not 2X 1 2X X∧ 5(0101) FTFT not 2X 2X 6(0110) FTTF 1X or 2X but not both 1 2X X∨ 7(0111) FTTT not ( 1X and 2X ) 1 2X X∨ 8(1000) TFFF 1X and 2X 1 2X X∧ 9(1001) TFFT 1X is 2X 1 2X X≡ 10(1010) TFTF 2X 2X 11(1011) TFTT If 1X then 2X 1 2X X⇒ 12(1100) TTFF 1X 1X 13(1101) TTFT if 2X then 1X 1 2X X⇐ 14(1110) TTTF 1X or 2X 1 2X X∨ 16(1111) TTTT always true ----- We illustrate the process with the data set in Table 2, which has 7n = total observations, 4nt = training observations, 7 4 3nc n nt= − = − = , and 4m = variables J. Farlow / Eur. J. Pure Appl. Math 169 Table 2: Sample Data Training Set ---------- Checking Set Y 1X 2X 3X 4X 1 0 0 1 0 0 1 0 0 1 0 0 1 1 1 1 1 1 0 1 1 0 1 0 1 0 1 0 0 1 0 0 0 1 1 Note that the dependent variable 1X did not predict Y in the first- and second observations (1s and 0s do not match), but do predict Y in the 3rd and 4th observations (1s and 0s match). Note, too that the logical relation 1 2X X∧ correctly predicts Y in the 1st, 2nd, and 3rd observations, but not the 4th, hence the recorded values 1, 1, 1, and 0 in the respective column, and a 3 in the bottom row illustrating the number of correct predictions. If one carries out this analysis for all 16 logical expressions of 1 2,X X , one arrives at the results in Table 3. Table 3: Correct and Incorrect Predictions of Y from 1 2,X X . 0 (0000) false 1 (0001) 1 2X X∧ 2 (0010) 1 2X X∧ 3 (0011) 1X 4 (0100) 1 2X X∧ 5 (0101) 2X 6 (0110) 1 2X X∨ 7 (0111) 1 2X X∨ 0 1 0 1 0 1 0 1 1 1 1 1 0 0 0 0 1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 2 3 1 2 1 2 0 1 8(1000) 1 2X X∧ 9(1001) 1 2X X≡ 10(1010) 2X 11(1011) 1 2X X⇒ 12(1100) 1X 13(1101) 1 2X X⇐ 14(1110) 1 2X X∨ 15(1111) tautology 0 1 0 1 0 1 0 1 1 1 1 1 0 0 0 0 1 1 0 0 1 1 0 0 1 1 1 1 1 1 1 1 3 4 2 3 2 3 1 2 Note the best logical predictor is 9(1001), which represents the logical function 1 2X X≡ . Table 3 is then computed for each of the ( 1) / 2 6m m − = pairs of dependent variables, recording only the last row of totals which represent the number of correct predictions of a given logical expression. Doing this we arrive at Table 4, where the second row in the table lists the best logical predictor of Y from among the 16 logical functions of the given pair of J. Farlow / Eur. J. Pure Appl. Math 170 variables, and the bottom row 3 gives the fraction (and percentage) of times the given logical function accurately predicted Y . Table 4: Results for the first iteration. Variable pairs 1 2,X X 1 3,X X 1 4,X X 2 3,X X 2 4,X X 3 4,X X Function 1 2X X≡ 1 3X X∧ 1 4X X∨ 2 3X X⇐ 2 4X X⇒ 3 4X X∨ # correct 4/4(100%) 3/4(75%) 4/4(100%) 2/4(50%) 3/4(75%) 2/4 (50%) The next step is to select the best m functions from Table 4 (4 in this example) that predicts Y . Step 2: (Replace the Original Data with the Best Logical Estimates) We now evaluate the best m logical functions in Table 4, which represent the best logical predictors of Y based on the observations in the training set. In this example they are 1 2X X≡ , 1 3X X∧ , 1 4X X∨ , 2 4X X⇒ which yields 4, 3, 4, and 3 correct predictions. Evaluating these functions at the n observations, we arrive at the 4 4n m× = × array in Table 5, which we call XNEW and are the evaluated best predictors of Y . Note that in the column under 2 4X X⇒ for the observations in the checking set, the values are 1, 1, and 1, which are the logical evaluations of the corresponding observations in the original data matrix X in Table 2. We now replace the data matrix X of independent variables with the newly computed matrix XNEW. Table 5: New Computed Data XNEW Replaces X . Training Set -------------- Checking Set 1 2X X≡ 1 3X X∧ 1 4X X∨ 2 4X X⇒ 1 0 1 0 0 0 0 1 0 1 1 1 1 1 0 1 1 0 0 1 0 0 1 1 0 0 0 1 Step 3: (When to Stop: Goodness of Fit) Steps 1 and 2 define the algorithm. To determine when the process is stopped, we compute the number of correct predictions RMAX of Y at the end of Step 2 for those observations in the checking set. This is an easy thing to do since best predictions reside in the first column of XNEW (and also the first column of X ). Experiments have shown that RMAX increases for a few generations and then begins to decrease as the model begins to over-fit the data. Hence, the rule is to stop the algorithm when the RMAX curve reaches its maximum and select the first column of XNEW as the best predictor. When the algorithm J. Farlow / Eur. J. Pure Appl. Math 171 stops, the columns of XNEW contain the computed values of the best m logical predictors of Y . This defines the algorithm. 4. Computations We tested the algorithm’s ability to find best logical relations of several independent variables from among 50m = independent variables with 250n = observations ( 150nt = , 100nc = ). All independent variables were binary 0 and 1 variables, each having probability 0.5. The dependent variable Y was generated as a logical expression of the dependent variables for some observations, and random 0s and 1s for other observations. In this way we could determine how well the GMDH algorithm works. To carry out the generation of the dependent variable Y , we first select a number 0 1q≤ ≤ , then generate a uniform random number r between 0 and 1, and then the observations , 1,...,iy i n= computed by the rule specified logical function when random 0 or 1 when i r q y r q ≥ =  < Note that when 0q = the dependent observations iy are the values of specified logical function of the dependent variables, and when 1q = the computed values of Y are random 0s and 1s and have no relation to the independent variables. The goal was to determine how well the algorithm could pick out the chosen logical relation when 0 1q< < . Experiment 1: The value 1q = was chosen with logical function 15 35X X⇒ . It was not surprising that at the first iteration the maximum logical function was 15 35X X⇒ and that it predicted the dependent variable 150 out of 150 times in the training set and 100 out of 100 times in the checking set. The second through fifth place logical relations were (2) 15 42X X⇒ , (3) 15 49not both and X X , (4) 15 16not both and X X , and (5) 16 35X X⇒ which is not surprising since the chosen logical relation 15 35 15 35X X X X⇒ ≡ ∪ so the variables 15 35 and X X would naturally come into play. The following Table 6 shows the results for the same logical function for different values of q after 1 iteration. Table 6: Correct predictions of Y after 1 iteration for different values of q . Value of q 1 0.9 0.8 0.7 0.6 Best logical functions 15 35X X⇒ 15 35X X⇒ 15 35X X⇒ 15 35X X⇒ 15 35X X⇒ Fraction of correct predictions in the training set 150/150 (100%) 143/150 (95%) 134/150 (89%) 132/150 (88%) 124/150 (83%) J. Farlow / Eur. J. Pure Appl. Math 172 Fraction of correct predictions in the checking set 100/100 (100%) 97/100 (97%) 89/100 (89%) 84/100 (84%) 79/100 (79%) Experiment 2: Using the same data and values of q as in Experiment 1, but with the new logical function ( ) ( )5 10 15 20y X X X X= ∪ ∩ ∪ , we arrived at the results in Table 7 after two iterations. In these cases the algorithm reached its maximum predictions after only two iterations. Note that when 0.5q = the maximum logical function was close to the entered logical function but not exact. This is not surprising since only 50% of the dependent variables were computed from ( ) ( )5 10 15 20iy X X X X= ∪ ∩ ∪ . Table 7: Results after 2 Iterations Value of q 1 0.8 Best logical function ( ) ( )5 10 15 20X X X X∪ ∩ ∪ ( ) ( )5 10 15 20X X X X∪ ∩ ∪ Fraction of correct predictions in the training set 150/150 (100%) 140/150 (93%) Fraction of correct predictions in the checking set 100/100 (100%) 87/100 (87%) Value of q 0.6 0.5 Best logical function ( ) ( )5 10 15 20X X X X∪ ∩ ∪ ( ) ( )10 20 10 15X X X X∪ ∩ ∪ Fraction of correct predictions in the training set 120/150 (80%) 110/150 (73%) Fraction of correct predictions in the checking set 74/100 (74%) 71/100 (71%) 5. Conclusions From the experiments performed we conclude that the GMDH algorithm is an effective tool in finding logical relations among predictor variables. J. Farlow / Eur. J. Pure Appl. Math 173 References [1] I. Ruczinski, C. Kooperberg, and M. LeBlanc. Logic Regression. Journal of Computational and Graphical Statistics 12(3), 475-511, 2003. [2] S. J. Farlow. Self-organizing Methods in Modeling. GMDH Type Algorithms, S. J. Farlow, editor. Marcel Dekker , 1984.