Instructions: For each of the following exercise, create an M-file to store the MATLAB commands. Copy and paste the M-file into a text document. Include in the text document the pictures produced by MATLAB. Resize and crop the pictures so that they do not take up too much space. If the question requires written answers, include them in the text file in the appropriate location. Make sure you clearly label and separate all the Exercises. Preview the document before printing and remove unnecessary page breaks and blank spaces.
1. Consider a set of data points (1, y1), (2, y2), ...(100, y100) where the y's are random numbers with a normal distribution with mean 0 and variance 1 (note that this implies that the probability of the y- values being less than zero is the same as the probability of being greater than zero).
(a) Based on the distribution of the points, what do you expect is the best straight-line model that can be fit to the data?
(b) To confirm numerically your answer to part (a), generate the data as follows:
x = [1:1:100]'; %note the prime
y = randn(size(x)); % a column vector of 100 standard normal values.
Plot the values as dots with the command plot(x,y,'.'). Follow Example 1 to find the best linear fit using MATLAB. Answer the following questions:
(i) What values do you obtain for the slope and y-intercept? Give the values in scientific notation with five digits of accuracy (use format short e).
Plot the linear fit together with the points (use q = x;)
(ii) Do the values of the slope and y-intercept confirm your answer to part (a)? Explain.
2. Download the file co2.dat from the web page. This file contains two columns of numbers: the year from 1959 to 2012 and the average annual carbon dioxide concentration in the atmosphere (in parts per million), as measured at the Mauna Loa observatory in Hawaii. The data come from the Earth Systems Research Laboratory (Global Monitoring Division) of the National Oceanic and Atmospheric Administration (which is the government agency charged with developing weather and climate forecast models; if you are interested in additional information and data set you can visit https://www.esrl.noaa.gov/gmd/ccgg/trends ). The estimated standard error in the observations is 0.12 ppm.
Save the file in the MATLAB working directory. Then load it with
dat = load('co2.dat');
Note: if you do not save the file in your working directory, you need to include the whole path to load it into MATLAB. For instance, if you saved your file in the C:\tmp directory you can load it by entering load('C:\temp\co2.dat').
This creates the 49 × 2 array dat whose first column is the year and whose second column is the CO2 concentration.
You may find it convenient to create two separate data vectors, as follows:
year = dat(:,1);
conc= dat(:,2);
Plot the data with
plot(year,conc,'o')
(a) Follow Example 1 to compute the best-fit line for these data. What are the coefficients c1 and c2? Give the values with five digits of accuracy. Plot the line in black ( use q = year;), together with the original data. Use 'linewidth',2 and axis tight in your plot.
(b) Find the best-fit quadratic polynomial to the CO2 data (make sure you change appropriately the matrix X). What are the coefficients c1, c2 and c3? Give the values with five digits of accuracy. Plot the fitted curve in red (use 'linewidth' ,2) together with the data points and the linear fit from part (a). Add a legend to the graph.
3. Among the important inputs in weather-forecasting models are data sets consisting of temperature values at various parts of the atmosphere. These values are either measured directly with the use of weather balloons or inferred from remote soundings taken by weather satellites. A typical set of RAOB (weather balloon) data is given next. The temperature T in Kelvins may be considered as a function of p, the atmospheric pressure measured in decibars. Pressure in the range from 1 to 3 decibars correspond to the top of the atmosphere, and those in the range from 9 to 10 decibars correspond to the lower part of the atmosphere.
|
p
|
1
|
2
|
3
|
4
|
5
|
6
|
7
|
8
|
9
|
10
|
|
T
|
222
|
227
|
223
|
233
|
244
|
253
|
260
|
266
|
270
|
266
|
Enter the pressure values as a column vector p (use the colon : operator), and enter the temperature values as a column vector T.
The goal of this problem is to find the cubic polynomial y = c1 + c2x + c3x2 + c4x3 that gives the best least square fit to the data.
(a) Follow Example 1 to find the coefficients of the best cubic fit (make sure you appropriately modify the matrix X).
Plot the data together with the cubic fit. Make sure you define your vector q so that the graph of the cubic is nice and smooth.
(b) Let X be the matrix you used in part (a). Here is an easier way to evaluate and plot the cubic polynomial:
c = X\T
c = c([4:-1:1]); q = 1:0.1:10;
z = polyval(c,q); figure plot(q,z,p,T,'o');
How do the values of c compare to the ones you found in part (a)?
How does the plot compare to the one you found in part (a)?
Note that the system Xc = T is overdetermined. In this case, the "\" command finds the least-squares solution to the system. In fact, for overdetermined systems, the "\" command is equivalent to solving the Normal Equations using the Cholesky decomposition, just like you did in the previous problems. Thus the vector c = X\T is the same one you found in part (a) by solving the normal equations.
The command polyval(c,q) evaluates the polynomial with coefficients given by the vector c in decreasing powers of q. Because the powers are decreasing we had to rearrange the entries in c using the command c = c([4:-1:1]);. Type help polyval to learn more.
Attachment:- Least Squares.pdf