********************************************************************** Homework Assignment A$ CSE 150, UCSD Fall 2002 Due: Saturday, November 23, 11:59 PM ********************************************************************** >>> SUMMARY: This assignment is meant to give you some experience with the behavior of the EM algorithm specialized to modeling density functions using mixtures of circular Gaussians. This is just the beginning! I expect to have a better version of this assignment later this weekend. THe second part of this assignment will be to use a perceptron to classify the points from one_gauss.dat and annulus.dat. annulus.dat forms a ring in 2-d (like a donut). one_gauss.dat forms the donut hole. This is usually a problem that a perceptron can't solve, but once we fit this data with gaussians, it should be trivial. (with some errors due to noise). If this sounds mysterious at this point, that's because I haven't explained much of this yet! But I think you have enough to at least start on the EM algorithm. >>> RULES: You may discuss your program with other students or the TA or the professor. However, you may not copy anybody else's program. In writing your code, you are on your own. However, I would prefer if you work in teams of two. >>> ASSIGNMENT: I will provide you with data generated from two circular gaussians. I will also provide you with data similar to the "annular" data given in Figure 2.14 of Bishop. (this is the donut). I will also give you data generated by one gaussian (the donut hole). Your job: Write code to implement the EM algorithm and experiment with fitting mixtures of circular Gaussians to this data. In order to observe the algorithm as it runs, I would like you to save the center and standard deviations of your Gaussians at each iteration of the algorithm. Plot your results using matlab's standard plotting tools as the algorithm progresses, on top of the data it is fitting. >>> WHAT TO TURN IN details of submission, what to name your functions, etc., will become available "soon." My first guess: Please tar your well-documented code files, including plots of the movements of Gaussians as they fit the two distributions, using two different initial parameter settings. Include a README file that: 1. Documents your approach, and any deviations from the algorithms given. 2. Explains the file organization of your code. 3. Explains how to run your code on the data. 4. Explains how to view the plots you made. 5. Documents any known bugs in your code. Here is the EM algorithm in pseudocode: 0. start with m gaussians. Initialize their parameters randomly, as follows: pick random points to be the centers (means) pick the standard deviations to be the distance to the nearest other center. pick the prior on each component to be 1/m. 1. E-step: for each data point x and each mixture component j, compute P(j|x) = p(x|j)P(j)/p(x); [note: p(x) = SUM_j p(x|j)P(j)] what is p(x|j)? that is the result of plugging x into the formula for a gaussian with the parameters of gaussian j. I.e., p(x|j) = 1/[(2PI*sigma^2)^d/2] exp(-||x-mu||^2/(2 sigma^2)) (where d is the dimensionality of the data - in our case, 2, and we are assuming all of the gaussian components are circular, and mu and sigma are the mu and sigma for component j.) Since we are working in 2d (d=2), this can be written more simply as: p(x|j) = 1/(2PI*sigma^2) * exp(-||x-mu||^2/(2 sigma^2)) ||x - mu||^2 is just the euclidean distance squared, i.e., for two-d data, (x_1 - mu_1)(x_1 - mu_1) + (x_2 - mu_2)(x_2 - mu_2) 2. recompute the mean, standard deviation, and prior on each component: mu_j [NEW] = SUM_n P(j|x_n)[OLD]*x_n ----------------------- SUM_n P(j|x_n)[OLD] where SUM_n means, sum over all of the data points, x_n is the nth data point. (sigma_j)^2 [NEW] = 1/d SUM_n P(j|x_n)[OLD]*||x_n - mu_j[NEW]||^2 ------------------------------------- SUM_n P(j|x_n)[OLD] P(j)[NEW] = 1/N * SUM_n P(j|x_n) [OLD] where N is the number of data points. Go to step 1. (repeat this, say, 10 times). That's it! COMMENTS: Note that if you run this with ONE component on the donut hole data, in one iteration, you should get the standard ML estimate for the parameters. This is a good way to check your code! You should get a probability of ONE for the component, and the algorithm should give you the sample mean and sample std deviation (divided by 2) for this two-d data. Another way to check your code is to compute the log likelihood of the data after every iteration. This should ALWAYS INCREASE IF YOUR CODE IS WORKING. In fact, we would really like it if your code did this!!! ________________________________________________ PART II 1. Run EM with k=7 for the donut, k=1 for the hole, separately. This is called "class conditional density estimation," and is reasonable given that I am about to give you a supervised training problem. Concatenate the parameters for the resulting 8 gaussians into one set of parameters. If you keep the data in an array for each parameter, this should be trivial. You could do this interactively if you like. Now, I will give you the same data from the donut and the hole, in one file, with a class label, 1 or -1. E.g., x y class [these strings won't be in the file] -0.419368 -0.393973 -1 1.019422 -0.914783 1 ... 2. read in the data, and compute P(x|j) for each (x,y) data point and each of the 8 components. Therefore, you will have 8 numbers for each data point. Note that since you computed these densities separately for each class in step 1 above, these should sum to 2 rather than 1. That doesn't matter. train a perceptron to classify the 8 numbers into the class 1 or -1, depending upon the class label. How is this done? The perceptron, as I will explain in class today, is a linear threshold unit. In this case, it will have 8 inputs and one output. Its parameters are 8 weights plus a threshold, which may be initialized to 0. For each point, you compute: net_input = (SUM (for j=1 to 8) of weight_j*input_j) - threshold. class = SIGN(net_input) if class =/= training label, then: foreach j weight_j = (training_label - class)*input_j; threshold = (training_label - class)*(-1); else /* do nothing */ This is repeated, say, 100 times, OR, until convergence. Convergence occurs when no updating is done. Unfortunately, you should not converge on this data, as it is noisy. That is, some points will actually overlap between the two classes. Report the number of incorrect classifications on every epoch (every pass through all of the data points). This number should drop quickly, and then bounce around at some low number because of the noise. The End.