Statistics - Study Guide

Chapter 4 Statistics In this course, we will be examining large data sets coming from different sources. There are two basic types of data: • Deterministic Data: Data coming from a specific function- each point has a unique image. If we know the function, we can predict exactly the data terms that we see. • Probabilistic Data: This data cannot be predicted exactly; with similar inputs, we might find very different outputs. The data is random in some sense. In this case, we can only hope to model some characteristics of the data, and not the exact data points themselves. We saw this earlier in the N −armed bandit problem, where we saw, with each pull of the arm, we would get different outputs. Thus, we were interested in extracting general information about the processes- the expected values, or mean payoff of each machine. In data analysis, most data falls in between the two classes- that is, data from most sources involves a deterministic part, and a random part (perhaps from measurement errors). We will not do a comprehensive listing of statistical procedures here. Instead, we introduce this material so that we have some statistical language to work with later on. 4.1 Measuring Variability Some data sets naturally have variation- For example, suppose I measure the temperature at the airport in Walla Walla on January 29th of each year. Here is the data for the maximum temperature recorded on January 29th for the years 1 2 CHAPTER 4. STATISTICS shown: 2008 2007 2006 2005 2004 39 30 55 50 57 2003 2002 2001 2000 1999 51 34 45 28 55 1998 1997 1996 1995 1994 46 30 18 61 43 1993 1992 1991 1990 1989 41 59 39 48 55 1988 1987 1986 1985 1984 54 53 36 32 52 Given the historical record, we might try to predict the temperature for next January 29th. What might our answer be? • The maximum is 61, the minimum is 18 and the average is 44.44 degrees, so we might guess that the temperature should be approximately 44.44 degrees. • We could make bins out of the data and look at the frequency of temperatures in each bin. For example, suppose we take 5 equally spaced intervals between 18 and 61. In each of these intervals, we count how many temperatures are present (you should verify a couple of them!): Interval 18 − 26.6 26.6 − 35.2 35.2 − 43.8 43.8 − 52.4 52.4 − 61 Count (out of 25) 2 4 5 6 8 Percentage 8 16 20 24 32 We might now say that 44.4 degrees was probably too cold- It looks like 32 percent of the time, the temperature on January 29th was between 52.4 and 61 degrees, and this “bin” is the most likely. We could continue the process of “binning” if we had a lot of data- In Figure 4.1 we use a computer simulation using random data to see what happens as we add more and more bins (or intervals). On the vertical axis, we are measuring the frequency rather than the percentage. By adding more data and more subintervals, it looks like the probabilities are approaching a continuous function. This function is what we define next: 4.2 Functions that Define Data The basic way of defining non-deterministic data is through the use of a probability density function1 , or p.d.f. Rather than defining a function in terms of inputs and outputs, a p.d.f. defines the probability of certain events occurring. To be more specific, a function f (x) is said to be a probability density function if it satisfies the following conditions: 1 Some texts use different vocabulary for discrete v. continuous data- We will not, as it should be clear from the context 4.2. FUNCTIONS THAT DEFINE DATA 3 Figure 4.1: The result of adding more and more data from a computer simulation. The vertical axis is the frequency count for each bin. We see that the histogram is beginning to approximate a continuous function. 1. f (x) is always non-negative. R∞ 2. −∞ f (x) dx = 1 3. The probability of an event between values x = a and x = b is given by: Z Pr(a ≤ X ≤ b) = b f (x) dx a From the definition, we see that the probability of any specific number is zero. Furthermore, in practice we won’t be dealing with continuous functions (although we might try modeling them); rather, we will be looking at discrete intervals. Therefore, suppose that our function f is non-zero on a finite interval (another way to say this is that f has bounded support). Then we can break the interval (or intervals) up into subintervals, and consider the probability of data occurring in each of these. In this case, the p.d.f. is discrete, and our definition changes somewhat: A discrete p.d.f. will be a finite set of numbers, {P1 , P2 , . . . , Pk }, so that: 1. Pi is non-negative, for all i. Pk 2. i=1 Pk = 1 3. The probability of a data value occurring in subinterval i (or bin i) is Pi . Note that in this case, each of our “events” (data in subinterval i) are disjoint, as the probability of landing on an endpoint is zero. In the N −armed bandit 4 CHAPTER 4. STATISTICS problem, we had some experience with these- in that case, Pi = π(i), which was the probability of choosing machine i. In Matlab, we visualize a PDF via the hist command, which produces a histogram of the input. Let’s take a look at some template probability distributions: • Example 1: The Uniform Distribution – The Continuous Version:  f (x) = 1 b−a , 0, if a ≤ x ≤ b otherwise – The Discrete Version (using N bins over the same interval):   b−a b−a 1 Pr a + (i − 1) ≤x≤a+i = = Pi , i = 1, 2, . . . , N. N N N – In Matlab, to obtain an m × n array of numbers from a uniform distribution over [0, 1], we type x = rand(m, n) or, for a single number: x = rand. Here is an example you should try. In this case, we use 5000 data points and 20 intervals (try changing them around): x=rand(5000,1); hist(x,20); • Example 2: The Normal (or Gaussian) Distribution – The Continuous Version: The Normal distribution with mean µ and variance σ 2 (to be defined shortly) is defined as:   1 −(x − µ)2 f (x) = √ exp σ2 2π σ This is the common “bell-shaped curve”; the constant in the front is needed to make the integral evaluate to 1. Note that the normal distribution with zero mean and variance 1 simplifies to: 2 1 f (x) = N (0, 1) = √ e−x 2π – In Matlab, we can obtain values from a normal distribution with zero mean and unit variance by x = randn(m, n) or x = randn. Try the same commands as before: x=randn(5000,1); hist(x,20); To change the mean, just add a constant. To change the standard deviation, just multiply: 4.3. MEASURING THE CENTRAL TENDENCY OF DATA 5 x=randn(5000,1); y=5*x+8; hist(y,20); mean(y) std(y) • Example 3: The p.d.f. commonly used to model the human voice is called the double Laplacian distribution:  Ke−|x| , −1 ≤ x ≤ 1 f (x) = 0, otherwise In the exercises, you’ll be asked to determine the value of K, and we’ll also see how the shape of the Laplacian compares to a normal distribution. In the meantime, you can visualize this in Matlab- The file laughter is built into Matlab: load laughter sound(y,Fs); %Plays the file over the speaker, if available hist(y,50) 4.2.1 The cumulative distribution function (cdf ) The cumulative distribution function (a.k.a. distribution function) is defined via the pdf: Z x f (t) dt F (X) = Pr(−∞ < X < x) = −∞ We note that: • By the Fundamental Theorem of Calculus, part I, F (x) is the antiderivative of f (x). • F (x) is strictly increasing, going from a minimum of zero to a maximum of 1. • We saw the discrete version of this when we were working with the N −armed bandit; we used the Matlab function cumsum to create the cumulative distribution. • To minimize confusion between terms, we’ll refer to F (x) as the cdf. 4.3 Measuring the Central Tendency of Data The most basic way to characterize a data set is through one number- the average, which can be defined a number of ways. Let us assume we have some data: {x1 , x2 , . . . , xm }: 6 CHAPTER 4. STATISTICS • The Sample Mean is: m x= 1 X xk m k=1 We’re all familiar with this definition- it’s just the average of the data points. For those who have had a course in statistics, you will recall that this is different than the theoretical population mean, which requires knowledge of the underlying p.d.f., which we seldom have. In that case, the definition (a.k.a. the Expected value) is given by: X µ = E(X) = xf (x) all x You might remember this formula by seeing that this is really a weighted average, with those events most probable getting a higher weight than events less probable. Also, if your data is being drawn independently from a fixed p.d.f., then the sample mean will converge to the population mean, as the number of samples gets very large. Suppose we have m vectors in IRn . we can similarly define the mean, just replace the scalar xk with the k th vector: m x= 1 X (k) x m k=1 In Matlab, the mean is a built-in function. The command is mean, and the output depends on whether you input a vector or a matrix of data: – For vectors, mean(x) outputs a scalar. – For an m × n matrix X, ∗ ∗ ∗ ∗ mean(X), the default, is the same as mean(X,1) mean(X,1); returns a row of n elements mean(X,2); returns a column of m elements The Matrix Grand Mean is the mean of all the data in the matrix (so the output is one number). Exercise: Show, algebraically, that the result of mean(mean(X)) will be the same as taking the average of all the matrix elements. We will look more closely at these commands below, but first we finish a few definitions: • The Sample Median is a number so that exactly half the data is above that number, and half the data is below that number. Although the median 4.3. MEASURING THE CENTRAL TENDENCY OF DATA 7 does not have to be unique, we follow the definitions below if we are given a finite sample: If there are an odd number of data points, the median is the middle point. If there is an even number of data points, then there are two numbers in the middle- the median is the average of these. Although not used terribly often, Matlab will perform the median as well as the mean: median(X); • The Mode rarely if ever comes into use as a measure of central tendency, but it is the value taken the most number of times. In the case of ties, the data is said to be multi-modal. 4.3.1 Some Matlab: There are times when we want to define an array as a block of numbers, then build a new array out of it. For example a matrix A is given and two different ways of building it in Matlab follow:     1.3 1.3 1.3 1.3 1.3  0.1   0.1 0.1 0.1 0.1    b= A=  2.1   2.1 2.1 2.1 2.1  −1.0 −1.0 −1.0 −1.0 −1.0 A = b ∗ ones(1, 4); A = repmat(b, 1, 4) Example: Build an array A by repeating the sub-array given in B three times down and two times across (so that the array A is 6 × 6).  B= 1 0 0 −1 2 1  SOLUTION (Check in Matlab): B=[1 0 2; 0 -1 1]; A=repmat(B,3,2) 4.3.2 Centering and Double Centering Data Let matrix A be n × m, which may be considered n points in Rm or m points in IRn . If we wish to look at A both ways, a double-centering may be appropriate. The result of the double-centering will be that (in Matlab), we determine  so that mean(Â, 1) = 0, mean(Â, 2) = 0 The algorithm is (in Matlab): 8 CHAPTER 4. STATISTICS %Let A be n times m [n,m]=size(A); rowmean=mean(A); A1=A-repmat(rowm,n,1); colmean=mean(A1,2); Ahat=A1-repmat(colmean,1,m); or, equivalently: %Let A be n times m [n,m]=size(A); colmean=mean(A,2); A1=A-repmat(colmean,1,m); rowmean=mean(A1,1); Ahat=A1-repmat(rowmean,n,1); Proof: For the first version (row mean first): Let A1 be the matrix A with the row mean b subtracted:   a11 − b1 a12 − b2 · · · a1m − bm  a21 − b1 a22 − b2 · · · a2m − bm    A1 =   .. ..   . . an1 − b1 an2 − b2 · · · anm − bm with n bj = 1X aij n i=1 Now define c as the column mean of A1 . Mean subtraction of this column results in the Â, written explicitly as:   a11 − b1 − c1 a12 − b2 − c1 · · · a1m − bm − c1  a21 − b1 − c2 a22 − b2 − c2 · · · a2m − bm − c2     =   .. ..   . . an1 − b1 − cn an2 − b2 − cn · · · anm − bm − cn By definition, the column mean of  isP zero. Is the new row mean zero? It is clear that the new row mean is zero iff k ck = 0, which we now show: n X Proof that ck = 0 k=1 We explicitly write down what ck is: m ck = 1 X (akj − bj ) m j=1 and substitute the expression for bj , m n 1 X 1X aij ck = akj − m j=1 n i=1 ! m = m n 1 X 1 XX akj − aij m j=1 mn j=1 i=1 4.4. MEASURING SPREAD- THE VARIANCE 9 Now sum over k: n X k=1 ck = n X k=1 n   m m X n X X 1 1  akj − aij  = m j=1 mn j=1 i=1 m m n n XX 1 XX akj − aij = 0 m mn j=1 i=1 j=1 k=1 It may be clear that these two methods produce the same result (e.g., row subtract first, then column subtract or vice-versa). If we examine the (i, j)th entry of Â, Âij = aij − bj − ci = aij − n m m X n X 1X 1 X akj − aik + ars n m r=1 s=1 k=1 k=1 Therefore, to double center a matrix of data, each element has subtracted from it its corresponding row mean and column mean, and we add back the average of all the elements. As a final note, this technique is only suitable if it is reasonable that the m × n matrix may be data in either IRn or IRm . For example, you probably would not double center a data matrix that is 5000 × 2 unless there is a specific reason to do so. Exercise: Experimentally verify the results of this section by performing (in Matlab) three ways of double-centering the data on a random 6 × 8 matrix A (it should not already have mean zero in either columns or rows- you should check that first). • Subtract the row mean of A, then compute and subtract the column mean. • Subtract the column mean of A, then compute and subtract the row mean. • Compute the row mean and column mean and overall mean of the matrix A. Subtract row means, then column means, then add in the overall mean. 4.4 Measuring Spread- The Variance The number that is used to describe the spread of the data about its mean is the variance: Let {x1 , . . . , xm } be as defined above. Then the Sample Variance is: m s2 = 1 X (xk − x)2 m−1 k=1 where x is the mean of the data. If we think of the data as having zero mean, and placing each data point in a vector of length m, then this formula becomes: 1 kxk2 m−1 Why does m − 1 appear in the fraction? s2 = 10 CHAPTER 4. STATISTICS • People who have had a course in statistics know that this makes s2 an “unbiased estimator” of the population variance. • In these notes, the actual variance will not be important- rather, we will be interested in comparing variances across data sets in which case the constant used for normalizing will not matter. While we’re defining the terms, we should distinguish between the sample variance and the actual population variance (as we did in the mean). The population variance is defined as the expected value of (x − µ)2 , E((x − µ)2 ) = X (x − µ)2 f (x) all x However, we seldom know the population variance in practice, so we will only use the sample variance. The Standard Deviation is the square root of the variance. Let’s take some template data to look at what the variance (and standard deviation) measure: Consider the data: 2 1 1 2 − , − , 0, , n n n n If n is large, our data is tightly packed together about the mean, 0. If n is small, the data are spread out. The sample variance is:   1 4+1+0+1+4 5 2 s = = 2 2 4 n 2n and this is in agreement with our heuristic: If n is large, our data is tightly packed about the mean, and the standard deviation is small. If n is small, our data is loosely distributed about the mean, and the standard deviation is large. Another way to look at the standard deviation is in linear algebra terms: If the data is put into a vector of length m (call it x), then the standard deviation is: s= √ 4.4.1 kxk m−1 Covariance and Correlation Coefficients If we have two data sets, sometimes we would like to compare them to see how they relate to each other. Definition: Let X = {x1 , . . . , xm } , Y = {y1 , . . . , ym } be two data sets with means µx , µy respectively. Then the covariance of the data sets is given by: m Cov(X, Y ) = s2xy = 1 X (xk − µx )(yk − µy ) m k=1 4.4. MEASURING SPREAD- THE VARIANCE 11 There are exercises at the end of the chapter that will reinforce the notation and give you some methods for manipulating the covariance. In the meantime, it is easy to remember this formula if you think of the following: If X and Y have mean zero, and we think of X and Y as vectors x and y, then the covariance is just the dot product between the vectors, divided by m: Cov(x, y) = 1 T x y m We can then interpret what it means for X, Y to have a covariance of zero: x is orthogonal to y. Continuing with this analogy, if we normalized by the size of x and the size of y, we’d get the cosine of the angle between them. This is the definition of the correlation coefficient, and gives the relationship between the covariance and correlation coefficient: Definition: The correlation coefficient between x and y is given by: ρxy Pm 2 σxy (xk − µx )(yk − µy ) = pPm k=1 = Pm 2 2 σx σy k=1 (xk − µx ) · k=1 (yk − µy ) Again, thinking of X, Y as having zero mean and placing the data in vectors x, y, then this formula becomes: ρxy = xT y = cos(θ) kxk · kyk 1 This works out so nicely because we have a m in both the numerator and denominator, so they cancel each other out. We also see immediately that ρxy can only take on the real numbers between −1 and 1. Some interesting values of ρxy : If ρxy is: 1 0 -1 Then the data is: Perfectly correlated (θ = 0) Uncorrelated (θ = π2 ) Perfectly (negatively) correlated (θ = π) One last comment before we leave this section: The Covariance and Correlation Coefficient only look for linear relationships between data sets! For example, we know that sin(x) and cos(x) (as functions, or as data points sampled at equally spaced intervals) will be uncorrelated, but, because sin2 (x)+ cos2 (x) = 1, we see that sin2 (x) and cos2 (x) are perfectly correlated. This difference is the difference between the words “correlated” and “statistically independent”. Statistical independence (not defined here) and correlations are not the same thing! We will look at this difference closely in a later section. 12 CHAPTER 4. STATISTICS 4.5 The Covariance Matrix If we have p data points in IRn , we can think of the data as a p × n matrix. Let X denote the mean-subtracted data matrix (as we defined previously). A natural question to ask is then how the ith and j th dimensions (columns) covaryso we’ll compute the covariance between the i, j columns to define: p 2 σij 1X = X(k, i) · X(k, j) p k=1 Computing this for all i, j will result in an n×n symmetric matrix, C, for which: 2 Cij = σij In the exercises, we have you show that we can conclude that C can be computed using the definition below: Definition: Let X denote a matrix of data, so that, if X is p × n, then we have p data points in IRn . Furthermore, we assume that the data in X has been mean subtracted. Then the covariance matrix associated with X is given by: C= 1 T X X p In Matlab, it is easy to compute the covariance matrix. For your convenience, we repeat the mean-subtraction routine here: %X is a pxn matrix of data: [p,n]=size(X); m = mean(X); Xm = X-repmat(m,p,1); C=(1/p)*X’*X; Matlab also has a built-in covariance function. It will automatically do the mean-subtraction (which is a lot of extra work if you’ve already done it!). C=cov(X); If you forget which sizes Matlab uses, you might want to just compute the covariance yourself. It assumes, as we did, that the matrix is p × n, and returns an n × n covariance- HOWEVER, it will divide by p − 1 rather than by p. This is not a big issue for the applications we will be considering- you may use either method, but be aware that there are differences in the actual algorithms. 4.6 Exercises 1. By hand, compute the mean and variance of the following set of data: 1, 2, 9, 6, 3, 4, 3, 8, 4, 2 4.6. EXERCISES 13 2. Obtain a sampling of 1000 points using the uniform distribution: and 1000 points using the normal distribution: x=rand(1000,1); y=randn(1000,1); Compare the distributions using Matlab’s hist command: hist([x y],100) and print the results. You’ll note that the histograms have not been scaled so that the areas sum to 1, but we do get an indication of the nature of the data. 3. Compute the value of K in the double Laplacian function so that f is a p.d.f. 4. Next, load a sample of human voice: load laughter If you type whos, you’ll see that you have a vector y with the sound data. The computers in the lab do have sound cards, but they don’t work very well with Matlab, so we won’t listen to the sample. Before continuing, you might be curious about what the data in y looks like, so feel free to plot it. We want to look at the distribution of the data in the vector y, and compare it to the normal distribution. The mean of y is already approximately zero, but to get a good comparison, we’ll take a normal distribution with the same variance: clear load laughter whos sound(y,Fs); %This only works if there’s a good sound card s=std(y); x=s*randn(size(y)); hist([x y],100); %Blue is "normal", Red is Voice Print the result. Note that the normal distribution is much flatter than the distribution of the voice signal. 5. Compute the covariance between the following data sets: x −1.0 y −1.3 −0.7 −0.7 −0.4 −0.1 −0.1 0.5 0.2 0.5 0.8 1.1 1.7 2.3 (4.1) 6. Let x be a vector of data with mean µ, and let a, b be scalars. What is the mean of ax? What is the mean of x + b? What is the mean of ax + b? (NOTE: Formally, the addition of a vector and a scalar is not defined. Here, we are utilizing Matlab notation: The result of a vector plus a scalar is addition done componentwise. This is only done with scalars- for example, a matrix added to a vector is still not defined, while it is valid to add a matrix and a scalar). 14 CHAPTER 4. STATISTICS 7. Let x be a vector of data with variance σ 2 , and let a, b be scalars. What is the variance of ax? What is the variance of x + b? What is the variance of ax + b? 8. Show that, for data in vectors x, y and a real scalar a, Cov(ax, y) = aCov(x, y) Cov(x, by) = bCov(x, y) 9. Show that, for data in x and a vector consisting only of the scalar a, Cov(x, a) = 0 10. Show that, for a and b fixed scalars, and data in vectors x, y, Cov(x + a, y + b) = Cov(x, y) 11. If the data sets X and Y are the same, what is the covariance? What is the correlation coefficient? What if Y = mX? What if Y = mX + b? 12. Let X be a p×n matrix of data, where we n columns of p data points (you may assume each column has zero mean). Show that the (i, j)th entry of 1 T X X is the covariance between the ith and j th columns of X. HINT: It p might be convenient to write X in terms of its columns, X = [x1 , x2 , . . . , xn ] Also show that p1 X T X is a symmetric matrix. 13. This exercise shows us that our geometric insight might not extend to high dimensional space. We examine how points are distributed in high dimensional hypercubes and unit balls. Before we begin, let us agree that a hypercube of dimension n has the edges: (±1, ±1, ±1, . . . , ±1)T so, for example, a 2-d hypercube (a square) has edges: (1, 1)T , (−1, 1)T , (1, −1)T , (−1, −1)T (a) Show that the distance (standard Euclidean) from the origin to a √ corner of a hypercube of dimension d is d. What does this imply about the shape of the “cube”, as d → ∞? (b) The volume of a d−dimensional hypersphere of radius a can be written as: Sd ad Vd = d where Sd is the d−dimensional surface area of the unit sphere. 4.7. LINEAR REGRESSION 15 First, compute the volume between hyperspheres of radius a and radius a − . Next, show that the ratio of this volume to the full volume is given by:   d 1− 1− a What happens as d → ∞? If we have 100,000 data points “uniformly distributed” in a hypersphere of dimension 10,000, where are “most” of the points? 4.7 Linear Regression In this section, we examine the simplest case of fitting data to a function. We are given two sets of dataX = {x1 , x2 , . . . , xn } Y = {y1 , y2 , . . . , yn } We wish to find the best linear relationship between X and Y . But what is “best”? It depends on how you look at the data, as described in the next three exercises. 1. Exercise: Let y be a function of x. Then we are trying to find m and b so that y = mx + b best describes the data. If the data were perfectly linear, then this would mean that: y1 = mx1 + b y2 = mx2 + b .. . yn = mxn + b However, most of the time the data is not actually, exactly linear, so that the values of y don’t match the line: mx + b. Thus we have an error: E1 = n X |yk − (mxk + b)| k=1 (a) Show graphically what this error would represent for one of the data points. (b) E1 is a function of two variables, m and b. What is the difficulty in determining the minimum2 error using this error function? 2 We could solve this problem by using Linear Programming, but that is outside the scope of this text. 16 CHAPTER 4. STATISTICS (c) Define the squared error as: Ese = n X (yk − (mxk + b))2 k=1 Why is it appropriate to use this error instead of the other error? (d) Ese is a function of m and b, so the minimum value occurs where ∂Ese =0 ∂m ∂Ese =0 ∂b Show that this leads to the system of equations: (the summation index is 1 to n) X X X m x2k + b xk = xk yk X X m xk + bn = yk (e) Exercise: Write a Matlab routine that will take a 2 × n matrix of data, and output the values of m and b found above. The first line of code should be: function [m,b]=Line1(X) and save as Line1.m. 2. Exercise: If we treat x as a function of y (i.e., x = f (y)), how does that change the equations above? Draw a picture of what the error would represent in this case. Write a new Matlab routine [m,b]=Line2(X) to reflect these changes. 3. Exercise: The last case is where we treat x and y independently, so that we don’t assume that one is a function of the other. (a) Show that, if ax + by + c = 0 is the equation of the line, then the distance from (x1 , y1 ) to the line is |ax1 + by1 + c| √ a2 + b2 which is the size of the orthogonal projection of the point to the line. This is actually problem 53, section 11.3 of Stewart’s Calculus text, if you’d like more information. (HINT: The vector [a, b]T is orthogonal to the line ax + by + c = 0. Take an arbitrary point P on the line, and project an appropriate vector to [a, b]T .) Conclude that the error function is: E= n X (axk + byk + c)2 k=1 a2 + b2 4.7. LINEAR REGRESSION 17 (b) Draw a picture of the error in this case, and compare it graphically to the error in the previous 2 exercises. (c) The optimimum value of E occurs where ∂E ∂c = 0. Show that if we mean subtract X and Y , then we can take c = 0. This leaves only two variables. (d) Now our error function is: E= n X (axk + byk )2 k=1 a2 + b2 Show that we can transform this function (with appropriate assumptions) to: n X (xk + µyk )2 E= 1 + µ2 k=1 (for some µ), and conclude that E is a function of one variable. (e) Now the minimum occurs where dE dµ = 0. Compute this quantity to get: µ2 A + µB + C = 0 P P 2 P 2 xk yk , xk , yk . This is a where A, B, C are expressions in quadratic expression in µ, which we can solve. Why are there (possibly) 2 real solutions? (f) Write a Matlab routine [a,b,c]=Line3(X) that will input a 2 × n matrix, and output the right values of a, b, and c. 4. Exercise: Try the 3 different approaches on the following data set, which represents heights (in inches) and weight (in lbs.) of 10 teenage boys. (Available in HgtWgt.mat) X Y 69 138 65 127 71 178 73 185 68 141 63 122 70 158 67 135 69 145 70 162 Plot the data with the 3 lines. What do the 3 approaches predict for the weight of someone that is 72 inches tall? 5. Exercise: Do the same as the last exercise, but now add the data point (62, 250). Compare the new lines with the old. Did things change much? 6. Matlab Note: Consider using the command subplot to plot multiple graphs on the same figure. For example, try the following sequence of commands: x=linspace(-8,8); y1=sin(x); y2=sin(2*x); 18 CHAPTER 4. STATISTICS y3=sin(x.*x); y4=sin(exp(-x)); subplot(2,2,1); plot(x,y1); subplot(2,2,2); plot(x,y2); subplot(2,2,3); plot(x,y3); subplot(2,2,4); plot(x,y4); 4.8 The Median-Median Line: The median of data is sometimes preferable to the mean, especially if there exists a few data points that are far different than “most” data. 1. Definition: The median of a data set is the value so that exactly half of the data is above that point, and half is below. If you have an odd number of points, the median is the “middle” point. If you have an even number, the median is the average of the two “middle” points. Matlab uses the median command. 2. Exercise: Compute (by hand, then check with Matlab) the medians of the following data: • 1, 3, 5, 7, 9 • 1, 2, 4, 9, 8, 1 The motivation for the median-median line is to have a procedure for line fitting that is not as sensitive to “outliers” as the 3 methods in the previous section. Median-Median Line Algorithm • Separate the data into 3 equal groups (or as equal as possible). Use the x−axis to sort the data. • Compute the median of each group (first, middle, last). • Compute the equation of the line through the first and last median points. • Find the vertical distance between the middle median point and the line. • Slide the line 1/3 of the distance to the middle median point. 4.8. THE MEDIAN-MEDIAN LINE: 19 3. Exercise: Understand the commands in the program mmline. There are several new matlab functions used. Below are some hints as to how we’ll divide the data into three groups. • If we divide a number by three, we have three possible remainders: 0, 1, 2. What is the most natural way of seperating data in these three cases (i.e., if we had 27, 28 or 29 data points)? • Look at the Matlab command rem. Notice that: rem(27,3)=0 rem(28,3)=1 rem(29,3)=2 • Look at the Matlab command sort. The full command looks like: [s,index]=sort(x) The output vector s will be x sorted. The vector index will contain which order the original indices were in. That is, x(index)=s • We can therefore sort x first, then sort y according to the index for x. 4. Exercise: Try this algorithm on the last data set, then add the new data point. Did your lines change as much? 5. Exercise: Consider the following data set [?] which relates the index of exposure to radioactive contamination from Hanford to the number of cancer deaths per 100, 000 residents. We would like to get a relationship between these data. Use the four techniques above, and compare your answers. Compute the actual errors for the first three types of fits and compare the numbers. County/City Index Deaths Umatilla 2.5 147 Morrow 2.6 130 Gilliam 3.4 130 Sherman 1.3 114 Wasco 1.6 138 Hood River 3.8 162 Portland 11.6 208 Columbia 6.4 178 Clatsop 8.3 210