One can write simple programs, or `scripts,' in MATLAB , which are simply lists of MATLAB commands written in a text file (called m-file for the `.m' file extension MATLAB uses) that can be called from within the MATLAB command window. To create one of these m-files, click on the File menu in the upper left corner of the command window, then `New' and `M-file.' A text-editor window should pop up, in which you can type a sequence of commands. For example, here is an m-file which plots a stability diagram for the ventilation/ model we studied in class (from Edelstein-Keshet):
alpha=1; m=1; beta=linspace(0.02,2,100); % set range for beta lam1=zeros(1,100); lam2=zeros(1,100); % initialize eigenvalues for j=1:100 % notice `...', for line continuation jacob=[ 1-sqrt(m*alpha*beta(j)), -sqrt(m*alpha/beta(j)); ... beta(j), 0]; % Jacobian for this beta eigenvalues = eig(jacob); % Find eigenvalues lam1(j)=eigenvalues(1); % First eigenvalue lam2(j)=eigenvalues(2); % Second eigenvalue end % end of loop plot(beta,abs(lam1),'m',beta,abs(lam2),'y'); hold on, plot(beta,1,'g--'), hold off % plot a dashed line at 1To save your m-file go to the `File' menu and click `Save As.' You must get on the C:/Temp disk, or use your own floppy, then save the m-file as something, say yourname1.m (I will use jim1 below). Within the MATLAB command window, you can now type
cd c:\temp % this changes directories to the location of your m-file jim1and the sequence of commands you put in the m-file will be executed.
EXERCISE 1: Change the above m-file to plot stability for a range of values from zero to 3, with and . Add labels to the horizontal and vertical axes. At what does the equilibrium become unstable? | ||
One of the reasons to use an m-file is the need for loops which will iterate a
certain activity (like predictions of our disease models). Recall that our
simple disease model was given by
inf=zeros(1,7) % initialize infections to be a vector of zeros z=3; inf(1)=1; % set initial condition for infections for n=1:6, % begin disease loop inf(n+1)=inf(n)+z/100*inf(n)*(51-inf(n)); end % of for loop observed=[1 2 7 15 29 46 51]; % a set of observations days=[0:6]; % make an 'x' variable for plot plot(days,observed,'r*',days,inf,'y') xlabel('Day'), ylabel('Infections') % label axes legend('Observed','Predicted') % put a figure legend inThis script will calculate the predicted course of the infection and compare it to the observations graphically. Also you see how to do two nice things to graphs: label them and provide a graph legend. You can actually move the legend! To see this, use your mouse to click on the legend in the figure and move it somewhere else in the plot.
EXERCISE 2: Write a MATLAB m-file to analyze the data from the (simple) class experiments, using the script given above as a base. Once this is done, try some of your data with more complicated disease dynamics. | ||
EXERCISE 3: Add a line or two to the disease m-file to calculate the value of
the
model. The formula for the of two vectors
of data, and is
You will probably want to use the MATLAB command mean; to check out what this does, type help mean. How do the correlation coefficients of the predictions for our disease models match your visual representations? Another approach is to use the MATLAB command corrcoef, which gives a matrix of correlation coefficients. Try sending two column vectors to corrcoef and see how that compares with your home-baked result. |
||
A basic problem we often face is determining what polynomial
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | |
Seconds, | 0 | 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 |
Height, | 14 | 10.6 | 8.5 | 6.7 | 4.9 | 3.5 | 2.2 | 1. | .5 | .1 |
In MATLAB we can accomplish this regression with the following commands:
height = [14 10.6 8.5 6.7 4.9 3.5 2.2 1. .5 .1]; time = [ 0 10 20 30 40 50 60 70 80 90]; polyfit(time,height,2)The first argument of polyfit is a vector of independent variables at which observations were made, the second is a vector of observations, and the third is an integer which defines the order of the polynomial being fit (in this case, 2, meaning quadratic). These commands return the output
ans = 0.0014 -0.2766 13.6509which tells us that , and . To plot this fitted quadratic against our data, we can use the commands
t=linspace(0,90,91); hfit= 13.6509-.2766*t+.0014*t.^2; plot(time,height,'y*',t,hfit,'r'), xlabel('Time'), ylabel('Height')Which looks pretty good - but remember, fits always look good.
EXERCISE 4: A godless capitalist business manager intends to model the revenue
to be gained by selling units of a product with the function
Data from the sales unit of his company for and include the points:
|
||||||||||||||
EXERCISE 5: One way to determine how well a model matches data is to plot the observed data against the predicted data. If these points fall on a line through the origin with slope 1 then the model has done perfectly. Plot your disease data against your model predictions and fit the resulting graph with a line. What is the slope and intercept, and what does this mean about your model? | ||
Often we want to fit data to an exponential model, for example when we are
looking at rates of cooling, radioactive decay, or population growth curves
for bacteria. In these circumstance we expect the observations, , to be
related to the independent variable, , by
Let's do logarithmic regression on the leaky bucket data above:
logh = log(height); polyfit(time,logh,1)which returns the output
ans = -0.0486 3.1558which tells us that if then and . Let's plot this against our data and the quadratic regression to see how it looks:
logfit=exp(3.1558)*exp(-.0486*t); plot(time,height,'y*',t,hfit,'r',t,logfit,'m') xlabel('Time'), ylabel('Height'), legend('Observations','Quadratic','Exponential')Yuck! This looks horrible! What went wrong?
EXERCISE 6: How can you re-do the exponential fitting so that there is no need to fit for ? Does this improve the previous result? Hint: What is the behavior if ? | ||
In Lay's Linear Algebra and Its Applications (pp. 414-420) is
discussed the linear algebra related to least-squares fitting data to
arbitrary linear models
y = [14; 10.6; 8.5; 6.7; 4.9; 3.5; 2.2; 1.; .5; .1]; X = [ 1 0 0^2; 1 10 10^2; 1 20 20^2; 1 30 30^2; 1 40 40^2; 1 50 50^2; 1 60 60^2; 1 70 70^2; 1 80 80^2; 1 90 90^2 ]; XT=X' % This is the transpose in Matlab A = XT*X c = inv(A)*XT*yThe output vector, c, should contain the coefficients of the quadratic polynomial; compare them to your polyfit results above to check your work.
EXERCISE 7: A healthy child's systolic blood pressure (in millimeters of Hg) and
weight (in pounds) are approximately related by the equation
Use the following experimental data to `determine parameters which give the best `fit:'
Estimate the systolic blood pressure of a healthy child weighing 100 pounds. Plot the data and the curve, indicating where the estimate for the 100 pound child lies. Does this estimate seem consistent with the data? |
||||||||||||||||||||