STA113: Lab 1 (Friday, 27 August, 2004)

Starting Matlab

Log onto the acpub unix workstation
Start Matlab by typing
```
 matlab 
```
Input commands to Matlab from its screen interactively and execute Matlab commands simply by typing them in, followed by Enter
```
 
>> 3+2
ans =
     5
```
You can suppress output to the screen by adding a semicolon after the statement
```
>> 3+2;
```
You can assign values to variable names with the equals sign, for future use. The assignment below is of a scalar quantity. Very shortly you'll be assigning vectors.
```
>> x=5;
>> x^2
ans =
    25
```

Vectors can be created with the following syntax: start:increment:end. The following will show you how to create vectors and plot graphs.

>> x = 0:10:50
>> y = x.^2
   % The .^ operator means element-wise squaring, as opposed to squaring
   % the vector, which would not be possible since the dimsions don't work out.
>> plot(x,y,'.')
   % The '.' argument tells MatLab to plot points instead of, say, lines or circles.
>> hold on
   % The hold on command indicates that the present graph should not be erased
   % when the next one is created.
>> plot(x,y+500,'go')   %g is for green and o is for circles
>> plot(x,y+1000,'r-')  %r is for red and - is for line segments.
>> hold off

Draw bar graph for motorcycle data from Example 1.12 (page 22). (The percent symbol is the "comment" symbol in MatLab. Anything following a % in a line is not read by MatLab. Such lines are given here to to indicate what you're doing.)
```
% input data
freq = [41, 27, 20, 18, 3, 11];
manufacturer = {'Honda', 'Yamaha', 'Kawasaki', 'Harley-Davidson','BMW', 'Other'};
total = sum(freq)
```
Note that the line above, because it didn't end in a semicolon, resulted in both the assignment of a number to the variable called "total" and the reporting of that value to you as output. Had you typed that line with a semicolon at the end, then it would have resulted in assignment only without output, just like the two preceding lines did.
```
% draw bar graph (vertical)
bar(freq)
ylabel('Frequency of Motorcycles');
title('Histogram for motorcycle Data');

% "gca" gets the handle for the current axes
set(gca, 'XTickLabel', manufacturer);
```
If you want to play around with some of the graphical features of MatLab--this is NOT expected for this lab but could come in handy later--try the command set(gca). You'll get a lot of characteristics of the graph that can be set manually. For example, set(gca,'color','yellow') makes the background of your graph yellow--probably not something you want, but you never know.
```
% draw relative frequency bar graph (horizontal)
barh(freq./total)
% (the operator ./ is element-wise division in a vector or matrix.)
set(gca, 'YTickLabel', manufacturer);
xlabel('Relative Frequency of Motorcycle');
```
Find the definition for Pareto diagram from exercise 1.30 on page 27.
```
% draw Pareto diagram
[newfreq, index] = sort(-freq);
% (type help sort to see what this command is doing.)
% (or you might be able to figure it out yourself by removing the
%  semicolon to stop suppression of the output.)
newfreq = -newfreq;
bar(newfreq)
set(gca, 'XTickLabel', manufacturer(index));
ylabel('Frequency of Motorcycles');
title('Histogram for motorcycle Data');
```

Click the following link to download the data for Example 1.11 on page 20 exp1-11.dat and save it as a file named "exp1-11.dat" in your home directory.

data = load('exp1-11.dat');  % load data

Summary Statistics

n = length(data)  % sample size               

mean(data)               
median(data)

range(data)                      
max(data) - min(data)
var(data)
sum((data-mean(data)).^2)/(n-1)
std(data)
sqrt(var(data))

prctile(data, 0:25:100)  % five-number summary

% (You might also like to see what happens if you just type 0:25:100 with no semicolon.)
% (Or you might try prctile(data,25).)

Histogram

Two matlab commands "hist" and "histc". Read the help file to learn the difference between the two commands.

hist(data,7);  % draw histogram with 7 classes

edg =[2 4 6 8 12 20 30];
% remove the semicolon if you wish to stop suppression of the
% output and see what the variable edg is that you've created.
[n,bin] = histc(data, edg);  %
rfreq = n/length(data);
width = [edg(2:7)-edg(1:6),1]';
h = bar(edg, rfreq./width , 'histc'); % draw density historgram with classes defined by edg
set(h, 'facecolor', 'green');  % change the color of the bins

% change axes labels and ticks 
set(gca, 'XTick', edg);
xlabel('Bond Strength')
ylabel('Density');

Boxplot

We will draw boxplot for the haircolor data shown in class.

% read data where 'NaN' means missing data
pain= [ 62 60 71 55 48; 63 57 52 41 43; 42 50 41 37 NaN; 32 39 51 30 35];
pain = pain'  % transpose of the original matrix
boxplot(pain)
 
set(gca, 'XTicklabel', {'Light Blonde', 'Dard Blonde', 'Light Brunette', 'Dark Brunette'})
xlabel('Hair Color');
ylabel('Pain threshold score');

Regression

Old Faithful is a geyser in Yellowstone National Park in Wyoming, USA. As it is a major tourist attraction, being able to predict the timing and length of the next interruption would be useful. The Old Faithful dataset contains data about the date of the observation, the duration of an eruption (in minutes) and the time between eruptions (also in minutes).

Use the oldfaith data to investigate the linear relationship between variable DURATION and TIME. We are interested in whether there is a linear relationship between them. And if so, how well does the best fitting straight line predict the TIME between eruptions, given the DURATION of the current eruption.

[date, duration, time]=textread('oldfaith.dat', '%d%f%d', 'headerlines',1);

% In the text above, the %d%f%d argument indicates the data types of
% the three columns.  The first and last columns matlab is to read
% as a signed integer value and the middle column matlab is to read
% as a floating point value.  The 'headerlines' argument indicates
% that the data are preceded by a row of non-data that is to be
% ignored.  This row, of course, contains column titles.  Type
% help textread for more information.

% The predictor variable  is going to be DURATION
% and the response variable  is TIME.
% First make a scatter plot of DURATION vs. TIME
plot(duration, time, '.');

% Now find the least squares regression line

X = [ones(length(duration),1) duration]; 
% X needs to be a matrix with a leading  column of ones.

[b,bint ,r,rint,stats] = regress(time, X, 0.05);  
% b :  estimates for the regression coefficients
% r : residuals
% stats(1) : R-square

b
stats(1)
SSE = sum(r.^2);
SST = var(time)*(length(time)-1);
1-SSE/SST

% Now add the regression line to the scatter plot.
x = [1.5, 5.5];
y = b(1)+ b(2).*x;
hold on;
plot(x , y, 'r-');
hold off;

who
% The command 'who' reports back a list of all the variables MatLab
% recognizes that you have defined.  (Just a useful tip.)

Graded Assignment

Wildlife biologists can fairly accurately determine the length of an alligator from aerial photographs or from a boat. Determining the weight of an alligator from a distance is much more difficult. Wildlife biologists in Florida captured 25 alligators in order to collect data and to develop a model from which weight can be predicted from length. The data set alligator.txt contains the resulting 25 measurements, the first variable is the alligator's weight (in pounds?) and the second is its length (in inches?).

Make histograms of the alligators' lengths and weights. Also determine the mean and median length and weight of the alligators.

Try the following two linear regression models: (i) x = length and y = weight; (ii) x = log(length) and y = log(weight). For each model, draw a regression plot (data points with fitted regression line) and a scatter plot of residual vs. x, and report the R-square value. Which model do you think is better? Why? Write an equation estimating weight explicitly as a function of length. If you used log(length) and log(weight) to determine your model, be sure your final model is an expression for weight, not for log(weight).

Write a Word document containing your graphs and your concise answers to all the questions above. Be sure your graphs have titles and clear axis labels and that the answers to the questions above are clearly stated. Be sure your name is on your paper along with your lab section. Try to keep it limited to one page, but if you must use two pages, staple them together. This is to be turned in on Friday, 10 September at the beginning of the lab section. You will have two lab periods to work on it: 27 August and 3 September. If you realize you will not have enough time to finish the assignment by then, you need to work on it outside of lab time.

Getting Help from MatLab

Nicer, windowed help system
```
 >> helpwin  
```
You can teach yourself virtually everything about MatLab once you learn how to get help. Don't hesitate to experiment and see what you get. Omitting the semicolon at the end of your lines is a good way to see what's going on every time you enter a command.

Two other ways to find help

>> lookfor key_words
>> help command_name

Print Graphics

Directly print from Graphics window : Graphics window > File > Print
Export graph to files (.eps) and then view them from your working directory: "ghostview mygraph.ps &"

Leaving Matlab

```
 
>> save mywork.mat
>> exit 
```
You can find you have a new file called "mywork.mat" in your current directory. Next time when you start matlab at the same directory, you can restore all your results by typing
```
 >> load mywork.mat
```