Friday, March 6, 2015

Lines in Pleasant Places

Lines in Pleasant Places

By Bobby Neal Winters


Lord, you alone are my portion and my cup;
   you make my lot secure.
The boundary lines have fallen for me in pleasant places;
   surely I have a delightful inheritance.
I will praise the Lord, who counsels me;
   even at night my heart instructs me.
--Psalms 16:5-7, NIV


It is rumored that deep within the Ozarks there is a place where those interested go to learn the art of making moonshine.  It is called the Moonshine Academy.  Every year a new class of students enter the academy, but not everyone decides to come back the following year.  The proprietors of the Academy are concerned about this because their aim is spread the art of moonshine making as far and wide as possible.  Those who leave the Academy tend to become Revenuers. They would like to know whether there might be a way to tell ahead of time which of their new students might be at risk of becoming Revenuers.
When each student comes to MA, a record is made of the length of their beards and the number of missing teeth.  The proprietors wish to make predictions based on those two characteristics, having noted that those who have become successful as moonshiners have certain standards with respect to them.  

Some Math

There are two groups, those who succeed and those who don’t. Call the group that succeeds the Elect and denote them by E. We will denote the rest by R, for Revenuer.   We will begin the discussion by assuming a simpler situation than what we have and supposing there is only one characteristic or, to use the standard nomenclature, random variable to measure.  Call it X.
Random variables are functions that take values on populations.  We will write X|E for X restricted only to members of the Elect and X|R for X restricted only to revenuers.  We hope that we won’t be piling on too much notation if we let mean(X|E) denote the mean of X for members of the Elect and let mean(X|R) be the mean of X for members of the Revenuers.
Suppose that we find mean(X|R)
Overlapping Normals.png
In the figure, we note mean(X|R)=3 and mean(X|E)=5.
In our dreams, we want a magic number C (for cutoff). Say a student named Maynard walks in and his random variable is X(Maynard).  We would like to be able to say that if X(Maynard)> C, then Maynard would be one of the Elect and if X(Maynard)
As Yogi Berra said, “Predicting is hard, especially about the future.” In making making a prediction about Maynard’s future, using only the value of X(Maynard), a procedure would be to simply let this magic number C be the midpoint between mean(X|E) and mean(X|R), respectively, and predicting using the method described above, with the realization that you are going to be wrong some of the time.
This is the same as saying predicting Maynard to be in the group E if X(Maynard) is closer to the mean of X on E than it is to the mean of X on R and to be in R otherwise.

Some More Math

Recall that the folks at the Moonshine Academy had two random variables they kept on each of their students: beard length and number of missing teething.  Mathematicians like to abstract, so we can think of each of the students at MA as a point in space and the values of the random variables coordinates for that point.  We can then plot them out on an X-Y axis.  The results might look like those below:
MA1.plot.png
As you see here, the data separates into two groups.  There is one in the lower left hand corner that appears to be centered at about (1,1) and another in the upper right that appears to be centered at about (2,2).
Notice that these two groups are more clearly separated on the 2-dimensional plane than either would be if we simply ignored one of the random variables.   Ignoring the vertical component, for example, would cause collisions between quite a few members of the two groups and there would be considerably more overlap between the 1-dimensional projections of the two groups.
If only there were some way to rotate the picture so that the line connecting the centers of these two groups so that it aligned with the horizontal.

Even More Math

There is, in fact, a way to perform such a rotation. One could apply a rotation matrix.  That’s right, matrices have uses beyond just things to be row-reduced to solve linear systems of equation.  Indeed, there real purpose is that of geometric transformation.  If you multiply a line by a matrix (by which I mean each point on that line), the result will be either a line or a single point.  If the determinant of that matrix is not zero, then the result will always be a line.  What is more, if you multiply a parallelogram of area A by a matrix of determinant D, you will get a parallelogram of area DxA.  And that’s what determinants are all about. Pretty cool if you ask me.
In this particular case, we probably could get away with multiplying by a rotation matrix that we pick by inspection.  However, there is a standard procedure, which we will describe in excruciating detail later, wherein we can obtain the matrix.  At this point, let’s pretend we’ve gone gone through it and obtained a matrix M.  We then take the a point with coordinates (x1,x2) and put it in the row matrix Y=[x1, x2].  We can then obtain the coordinates of the transformed point as YxM, using matrix multiplication.   Applying this to every point gives the picture below:
MA1.rotated.png
Note that these two blobs are fairly well separated. One could simply forget about the vertical component  and project the points onto the horizontal axis in order to go through the procedure described at the first of this article.
I like to think of the groups in the above picture as spongy clouds.  They have centers that you can obtain simply by taking the averages of all the of coordinates of each.  The centers of the groups above are exactly the same as the rotated versions of the  centers of the groups in the previous picture.

The Process from Thirty-thousand Feet

I am now going to give you the big picture by telling a few lies.  Later in the article, if you are still awake, I will provide as many details as I can.
Suppose that a new student comes to seek admission to the Moonshine Academy.  The folks at the Admissions Office will measure his beard and count the number of teeth his as missing. (Sometimes they will count the number of teeth he has and subtract that from 32 because it is easier, but I digress.)   They will then measure the distance between that point and the centers of the two groups.  The student will be associated with the group to whose center he is closer.  That is if he is closer to the center of the Elect group he will be predicted to be Elect; if he is closer to the center of the Revenuer group he will be predicted to be a Revenuer.
The word predicted is very important in the paragraph above.  Predictions can be wrong.  Indeed, even when you use a procedure on data where you know the answers, some of the predictions will be wrong.  For example, for the data above we get the following table:



Predicted
Actual
E
R
E
799*
1
R
2
198*

This says that my model predicted 801 (799 + 2) be in the Elect; of these 799 were but two actually flunked out and became Revenuers.  It predicted 199 (1+198) to become Revenuers.  The model was right in 997 (799*+198*) of its predictions for its founding data. (That is why those numbers are marked with a *.)   These are incredibly good results and this will certainly not be the case with all data sets.

What Happened to the Rotation?

You might have noticed that I’d made a big deal about rotating the matrices to better see how nicely separated the groups were along the horizontal axis, but when I described the process the groups weren’t actually rotated.  This is because you don’t have to rotate for the test. The rotation does not change the distances.  What would be the closest center in the rotated picture would be closest in the unrotated picture.
However, once the initial test is run and you want to make predictions about candidates for next year it is more convenient to have a formula.  For the data above the formula would be: f(bl,mt)=0.982*bl+0.188*mt. Here, bl=beard length and mt=missing teeth.  It is easier to determine whether f(bl,mt)>2.330 than it would be to calculate the distance to each of the centers of the groups.  
And, as a student who would want to look good with respect to these standards, it is easier to tell that you’d be best advised to spend your time growing your beard as opposed to pulling your teeth.

How Do We Get That Formula?

To get to that formula, we must enter into the jungle of matrices where many enter and few emerge unchanged.
The first matrix I want to introduce you to is the covariant matrix. There is far more to know about it than what I am going to tell you.  This is how you get a covariant matrix.  Let X be a matrix each of whose columns are values of a random variable.  So column one will contain the beard length of each of our students and column two will contain the number of missing teeth of the same student. If you have more random variables to study, you will have more columns.  For us, with the beard lengths and the missing teeth of 1000 students, we have X as a 1000x2 matrix.
Now, calculate the mean of each column and subtract that number from each of the columns.  Call this new matrix Y. So each element of Y will be the difference between the value of a particular random value for a particular student and the mean of that random variable for all the students.  Now the covariant matrix of X, cov(X), is transpose(Y)*Y.  For us it will be a 2x2 matrix.  Note that cov(X) is symmetric; it is important.
Now, forget about this particular matrix because we are not going to use it, but remember how we got it.
Let XE be the data for the Elect and let XR be the data for the Revenuers. Let pE be the proportion of the data that is from the Elect and pR be the proportion from the Revenuers.  Let Sw=pE*cov(XE)+pR*cov(XR). Note that Sw is symmetric; it is important.
Now remember that first matrix X, the one I told you to forget about.  We are going to fix it.  Instead of being the value of the random variable minus the mean whole column, let it mean of the column for the data group minus the mean of the whole column.  That is to say if the student from row i is in the Elect, his entry in the Beard Length column will be the average beard length for the Elect minus the average beard length for everybody.  Then Sb=cov(X).  Note that Sb is a symmetric matrix; this is important.
Now let M=inverse(Sw)*Sb.  Inverse(Sw) is symmetric so M is symmetric.  This is important. Write it on your hand with a sharpie.
We are about to digress.

Eigen What?

A discussion of eigenvalues and eigenvectors would take us too far afield.  One or more of us would be dead by the time I finished.  Suffice it to say that the theory is beautiful in the same way the Arctic is: The beauty is real, but seeing it face to face comes at a cost.
It is enough for us to know the following. If M is a symmetric matrix, and it is, then there is a matrix P and a diagonal matrix D such that M=transpose(P)*D*P.  Furthermore, the matrix P will consist of the eigenvectors of M, whatever the hell they are, and the diagonal entries of D will be the eigenvalues of M, whatever the hell they are.  Take the eigenvector corresponding to the eigenvalue of largest absolute value and multiply it by the coordinates of the student, and you have your formula.  For us that is [bl , mt] * transpose([0.982, 0.188]).  To get the critical point C, one applies this formula to the midpoints of each of the groups, and takes the average of the two.


References

“Linear Discriminant Analysis-A Brief Tutorial,”  S. Balakrishnama and A. Ganapathiraju, Institute for Signal and Information Processing

Discriminant Analysis, William. R. Kleka, A Sage University Paper, Series: Quantitative Applications in the Social Sciences, 1980

Friday, January 23, 2015

M&M ANOVA

Khan Academy does a nice job of explaining ANOVA at this link.  This is in fact where I learned it.  Below I have a nice application of ANOVA using M&Ms that I would like to share.

There are numerous tables below which can be glided over without any loss of understanding.  Indeed, if you just read the prose in between the tables you will be far better off.

In the Fall of 2014, I assigned a series of activities to my Elementary Statistics involving M&Ms.  These activities begin here. There were six groups of students involved and each group took a sample of ten bags of M&Ms.  These samples are listed below.

Group 1
Color/Bag 1 2 3 4 5 6 7 8 9 10
Red 2 5 2 4 3 7 3 1 7 5
Orange 7 3 6 5 6 7 5 7 6 9
Yellow 2 3 0 1 2 0 2 4 3 4
Green 3 0 6 6 3 2 1 3 4 2
Blue 3 5 4 2 4 7 4 5 4 5
Brown 1 1 1 0 1 3 1 1 1 3

Group 2
Color/Bag 1 2 3 4 5 6 7 8 9 10
Red 2 2 3 2 2 3 1 6 3 3
Orange 3 4 0 4 4 4 3 2 2 3
Yellow 2 5 5 1 1 1 0 1 2 4
Green 3 1 5 4 4 3 4 2 4 2
Blue 3 2 2 4 3 3 6 1 4 4
Brown 3 2 1 1 2 2 2 4 1 1
Group 3
Color/bag 1 2 3 4 5 6 7 8 9 10
Red 1 1 1 0 2 3 0 4 0 0
Orange 2 0 1 2 5 6 5 3 1 0
Yellow 1 1 2 0 3 5 0 7 1 3
Green 1 3 1 2 2 1 5 4 1 1
Blue 2 2 4 2 4 5 5 2 2 1
Brown 1 0 0 1 2 0 3 1 2 1
Group 4
Color/Bag 1 2 3 4 5 6 7 8 9 10
Red 4 5 1 3 1 0 2 2 3 2
Orange 1 1 7 3 3 6 3 5 4 6
Yellow 1 4 6 0 3 2 3 0 3 4
Green 2 3 1 8 4 4 5 6 1 3
Blue 5 1 1 5 4 2 3 3 3 2
Brown 2 3 4 1 3 3 2 1 2 0
Group 5
Color/Bag 1 2 3 4 5 6 7 8 9 10
Red 2 3 1 5 5 2 0 1 5 5
Orange 2 0 1 3 2 3 3 1 5 3
Yellow 6 2 5 2 3 4 7 7 1 5
Green 2 5 5 7 1 2 4 6 5 1
Blue 1 5 6 2 4 5 4 1 2 3
Brown 3 6 0 1 3 2 3 2 1 1
Group 6
Color/Bag 1 2 3 4 5 6 7 8 9 10
Red 1 2 2 4 3 3 1 2 1 7
Orange 3 3 1 1 3 1 3 2 1 2
Yellow 3 5 3 3 2 2 3 6 4 4
Green 3 1 7 3 3 6 2 2 4 2
Blue 4 4 1 4 3 3 3 4 4 3
Brown 2 1 4 1 4 2 4 2 5 1
This is a nice collection of real data and my thought was to make the most of it.  As a sample size of ten is small, my thought was to pool the data, but before this can be legitimately done, it must be justified.  One might argue that since all of the samples were taken from M&M's that might be justification enough, but I had lingering doubts.  What if proportion of M&M color is not consistent from batch to batch? What if M&Ms are put out in a variety of Fun Sizes?  What if my students had just royally goofed?  In order to be careful, I decided that after having taught elementary statistics for twenty years it was time to learn ANOVA.

I first wanted to get an good confidence interval for the average number of M&Ms per bag. I calculated that using the data for each group, finding the bag by bag total.  I put that into the following table:

Bag/Group 1 2 3 4 5 6
1 18 16 8 15 16 16
2 17 16 7 17 21 16
3 19 16 9 20 18 18
4 18 16 7 20 20 16
5 19 16 18 18 18 18
6 26 16 20 17 18 17
7 16 16 18 18 21 16
8 21 16 21 17 18 18
9 25 16 7 16 19 19
10 28 17 6 17 18 19
I then calculated the mean for each of the groups individually and the grand mean of the total pooled data.  Using this, I calculated the sum of the squares for differences within each of the groups and the sum of the squares for differences between the groups.  Those calculations are in the table below:

DATA
Bag/Group 1 2 3 4 5 6
1 18 16 8 15 16 16
2 17 16 7 17 21 16
3 19 16 9 20 18 18
4 18 16 7 20 20 16
5 19 16 18 18 18 18
6 26 16 20 17 18 17
7 16 16 18 18 21 16
8 21 16 21 17 18 18
9 25 16 7 16 19 19
10 28 17 6 17 18 19
GrandMean Means
17.07 20.7 16.1 12.1 17.5 18.7 17.3

SSW 
1 2 3 4 5 6
7.29 0.01 16.81 6.25 7.29 1.69
13.69 0.01 26.01 0.25 5.29 1.69
2.89 0.01 9.61 6.25 0.49 0.49
7.29 0.01 26.01 6.25 1.69 1.69
2.89 0.01 34.81 0.25 0.49 0.49
28.09 0.01 62.41 0.25 0.49 0.09
22.09 0.01 34.81 0.25 5.29 1.69
0.09 0.01 79.21 0.25 0.49 0.49
18.49 0.01 26.01 2.25 0.09 2.89
53.29 0.81 37.21 0.25 0.49 2.89
Sums
156.10 0.90 352.90 22.50 22.10 14.10

SSB
1 2 3 4 5 6
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
13.20 0.93 24.67 0.19 2.67 0.05
Sums
132.01 9.34 246.68 1.88 26.68 0.54
The SSW sums to 568.6 and the SSB sums to 417.1.  The numerator has m-1=5 degrees of freedom as we are comparing m=6 groups. The denominator has m*(n-1)=6*(10-1)=54 degrees of freedom as each of those groups took a sample of size n=10.  This gives an F test-statistics of F=7.92.  The critical number for those degrees of freedom with a significance level of  alpha=0.10 is 1.957.  As 7.92 is greater than 1.957, we must conclude that these samples are not all drawn from the same population.

This came as something of surprise to me.  As an educator of over 30 years experience, I immediately suspected student error.  Looking at the SSW and SSB table above, I noted that the numbers from group 3 were considerably larger than the rest.  I was curious as whether and how they had erred.  Discerning this was easy because I had had the students document their process.  In looking at the documentation from group 3, I found the follow photograph:




The student had been told to use Fun Size M&Ms.  It was assumed that they would plain and that the bags would not be mixed.  We are well tutored in how one spells ass-u-me. 

I would be remiss at this point, however, if I did not say that I had pushed this further. Elementating group 3 does not fix the problem.  The remaining groups are not sampling the same populations and an examination of the documentation of the other groups does not reveal a similar glaring error in methods.   Of all six groups, only 4 and 6 seem to be sampling the same population.

I will be having my class do a similar experiment this semester--with better instructions from the teacher--and after this I will conduct this study again.

Wednesday, January 21, 2015

M&M Binomial

This uses data gathered from the M&M Activity.

Let's begin this with a thought experiment.  Imagine you have the job of filling bags with M&Ms.  One might imagine that there is a huge bin that has been filled with M&Ms and that you are just parcelling them out into the bags.

There is a mountain of M&Ms and a certain proportion of these are red, orange, yellow, green, blue, and brown.  The number is so large that the act of choosing, say, a red M&M on one trial does not appreciably reduce the probability of getting one on the next trial.  All of this argues that, if one is interested in particular colors, each trial is a Bernoulli Trial with a fixed probability of success.

Each M&M is about the same weight and you are aiming to fill your bag to at least a certain weight but not more than one M&M more than that.  This results in the bags having a small variation in number of M&Ms per bag.  All of this means that we have a, more or less, fixed number of M&Ms per bag.

We combine these two observations, and what we have looks astonishingly like a Binomial Experiment.

Using a sample of 10 bags of Plain M&Ms, I investigated whether they did in fact follow the binomial distribution for red M&Ms.  This sample contained a total of 182 M&Ms of which 28 were red.  This yields a proportion p=0.154 of red M&Ms. The ten bags contained from 17 to 19 each. I therefore chose to set n=19.

Given those, I used my data to create a frequency table.  In the table below, the first column is the possible number of successes for a binomial experiment with 19 trials, that is to say the numbers from 0 to 19 inclusive; keeping with standard notation, that column is labeled x.  The second column is the number of bags that had that many red M&Ms.  As we are setting up do to a chi square test to see if the binomial model fits, we have labeled the frequency column with an O.

x O
0 1
1 2
2 2
3 0
4 4
5 0
6 1
7 0
8 0
9 0
10 0
11 0
12 0
13 0
14 0
15 0
16 0
17 0
18 0
19 0

We ask the question of how well this fits with the expectations of the binomial distribution B(19,0.154).  I do the sample calculation for x=3 below:

Putting this calculation into a spread sheet, I obtain:

n O E
0 1 0.42
1 2 1.45
2 2 2.36
3 0 2.44
4 4 1.77
5 0 0.97
6 1 0.41
7 0 0.14
8 0 0.04
9 0 0.01
10 0 0.00
11 0 0.00
12 0 0.00
13 0 0.00
14 0 0.00
15 0 0.00
16 0 0.00
17 0 0.00
18 0 0.00
19 0 0.00

We can then compare the values of the O and the E columns by using the (O-E)^2/E measure.  The results are in the table below:

n O E (O-E)^3/E
0 1 0.42 0.81
1 2 1.45 0.21
2 2 2.36 0.06
3 0 2.44 2.44
4 4 1.77 2.80
5 0 0.97 0.97
6 1 0.41 0.85
7 0 0.14 0.14
8 0 0.04 0.04
9 0 0.01 0.01
10 0 0.00 0.00
11 0 0.00 0.00
12 0 0.00 0.00
13 0 0.00 0.00
14 0 0.00 0.00
15 0 0.00 0.00
16 0 0.00 0.00
17 0 0.00 0.00
18 0 0.00 0.00
19 0 0.00 0.00
Note that the sum of the (O-E)^2/E column is 8.32.  The critical number for the chi square test for this is 30.14.  This is the table look up with 19 degrees of freedom and a significance level of 0.05.

As 8.32 is not larger than 30.14, we cannot reject the null hypothesis of the chi square test.  The null hypothesis is that the model fits.  We've no proven that it is binomial, but we can say that our data is not inconsistent with the binomial distribution B(19, 0.154).

Do this test with the data you collected from the M&M Activity.