Curve Fitting
Regression with more general functions
Skills to work on in this section:
• Working with nonlinear models to fit data;
• Writing and calling ‘functions’;
• Using ‘fzero’ function from MATLAB;
• Organizational skills in coding.
Yearly temperature fluctuations are often modeled by trigonometric functions, which
lead to more difficult regression analyses. Consider the case of monthly average
maximum and minimum temperatures in Hartford. The table below lists the raw
temperature and then a scaled temperature. In order to model this data with a
simple trig. model, we’ll subtract the mean and divide by the maximum absolute
value to a
ive at the scaled temperatures that vary like sin and cos (between -1 and
+1).
One model for the scaled data is T (m) = sin(m � a), where T is the scaled temper-
ature, m is a scaled month (m 2 [0, 2⇡]), and a is a shift coefficient. With this, we
have reduced the model to a single parameter, a, we need to find which minimizes
some measure of e
or.
The e
or we want to minimize looks like:
E(a) =
nX
i=0
[Ti � (sin(mi � a))]2 .
What derivative do we need to take to minimize E? Is this a ‘partial’ derivative, o
a ‘full’ derivative?
This gives us a nonlinear alge
aic equation to solve for a. In general, nonlinear equa-
tions are difficult (or impossible!) to solve analytically. However, we can solve them
numerically with some of MATLAB’s built-in functions, like the ‘fzero()’ function.
1
Day 2 lecture: More Curve Fitting 2
Avg. High Scaled High Avg. Low Scaled Low
Jan XXXXXXXXXX
Feb XXXXXXXXXX
March XXXXXXXXXX
April XXXXXXXXXX
May XXXXXXXXXX
June XXXXXXXXXX
July XXXXXXXXXX
Aug XXXXXXXXXX
Sept XXXXXXXXXX
Oct XXXXXXXXXX
Nov XXXXXXXXXX
Dec XXXXXXXXXX
Mean ⇠60.92 – 41.25 –
To go from the model with the scaled temperature values back to the ‘normal’ tem-
peratures (that people are used to working with), how do we transform the scaled
function?
Listing 1. Hartford Weather Data
1 highs = [36, 39, 48, 60, 71, 79, 84, 83, 75, 63, 52, 41];
2 avgHigh = mean(highs);
3
4 % scaledHighs = highs � avgHigh;
5 % scaledHighs = scaledHighs./max(scaledHighs);
6
7 months = 1:12;
8
9 aCoefHigh = fzero(@hartfordHighs,1);
10 modelHighs=max(abs(highs � avgHigh))*sin(2*pi*months/12� aCoefHigh)+
avgHigh;
11
12 plot(months,modelHighs)
13 hold on
14 scatter(months,highs,'o','filled')
Day 2 lecture: More Curve Fitting 3
Listing 2. function; hartfordHighs
1 function value = hartfordHighs(a)
2 scaledHighs = [
3 �1.0000, �0.8796, �0.5184, �0.0368, 0.4047, 0.7258, ...
XXXXXXXXXX, 0.8863, 0.5652, 0.0836, �0.3579, �0.7993 ];
5
6 value = 0;
7 for k=1:12
8 m=2*pi*k/12;
9 value = value + (scaledHighs(k)�sin(m�a))*cos(m�a);
10 end
11 end
GROUP WORK: (Temperature example)
Work through the same example using the average low temperatures. Save your code
as .m files and any written/typed work as a .pdf. Submit code and written parts to
Moodle.
(1) Find the scaled average low temperatures.
(2) Create a function similar to function value = hartfordHighs(a) above fo
the average low temperatures in Hartford.
(3) Create a script in Matlab similar to Listing 1 that will take in the low tem-
peratures, scale them, fit a sin curve to the data, and plot the information.
(This should recover the same curve as shown above.)
Day 2 lecture: More Curve Fitting 4
Consider the US population (in millions) in the table on the next page. We will
determine an equation to model the population.
If p(t) is the population at time t, the logistic model of population growth is given by
the ordinary differential equation:
dp
dt
= r p
⇣
1� p
K
⌘
, p(0) = p0 .
The population p(t) needs to satisfy the equation and the initial condition. Since
the differential equation will give us a solution up to a constant, the initial condition
allows us to find a unique solution to the problem. What does this equation say about
different populations?
What happens to the change in the population dpdt if the population is p(t) = 0?
What happens to dpdt if the population is large, like p(t) � K?
What happens to dpdt if the population is moderately sized, say p(t) = K?
What happens to dpdt if the population is less than above, like p(t) ⇡ K/2?
The equation above has an analytic solution of p(t) =
p0 K
(K � p0)e�rt + p0
. Without
actually doing so, how could we show this satisfies the equation above?
In the model and solution above, the coefficients p0 represents the initial population,
K is the ca
ying capacity, and r is the growth rate of the population.
Day 2 lecture: More Curve Fitting 5
With this example, we have a good guess about what the solution looks like– a logistic
model. However, we have to find the coefficients r,K, p0 for the model that minimize
the e
or.
• create a function to output population at some given time t with some given coef-
ficients r,K, p0;
• Use the MATLAB function ‘lsqcurvefit’ to find the optimal choices for the coeffi-
cients r,K, p0;
• The function ‘lsqcurvefit’ requires an initial guess for the coefficients. How could
you determine an appropriate initial guess?
• The function ‘lsqcurvefit’ also has an optional lower and upper bound that it can
satisfy in determining the coefficients. What would an appropriate lower bound
e for each? Would an upper bound be useful to specify?
Year XXXXXXXXXX XXXXXXXXXX
Population XXXXXXXXXX XXXXXXXXXX.2
Year XXXXXXXXXX XXXXXXXXXX
Population XXXXXXXXXX XXXXXXXXXX XXXXXXXXXX
Listing 3. function; logistic
1 function P = logistic(p,t)
2 % inputs: vector p = ( p(1), p(2), p(3) ), scalar t
3 % p(1): growth rate,
4 % p(2): ca
ying capacity, K
5 % p(3): initial population, p0
6 % t: time (in years AFTER 1860)
7 % output: P, population at time t
8
9 P = p(2).*p(3)./((p(2)�p(3)).*exp(�p(1).*t)+p(3));
10
11 end
Day 2 lecture: More Curve Fitting 6
Listing 4. Population Model
1 decades=0:10:150;
2 population=[
3 31.4, 35.6, 50.2, 63.0, 76.2, 92.2, 106.0, 123.2,...
4 132.2, 151.3, 179.3, 203.2, 226.5, 248.7, 281.4, 308.7 ];
5
6 initGuess = [.01, 1000, 3.93]; % initial guess for parameters, r, K,
p0
7
8 optimset lsqcurvefit;
9 lowe
ound=[0, 0, 0]; % lower bound for parameters in initGuess vecto
10 options=optimset('TolFun',1e�8,'Display','off'); % setting tolerance
for solving, and turning output off
11
12 [p, e
or] = lsqcurvefit(@logistic,initGuess,decades,population,
lowe
ound,[],options); % getting r, K, p0 coefficients that
minimize e
o
13
14 modelPopulation=logistic(p,decades); %getting vector of data from
model using the function we created
15
16 plot(decades,modelPopulation)
17 hold on
18 scatter(decades,population,'o','filled')
Day 2 lecture: More Curve Fitting 7
GROUPWORK: (Population example)
Save your code as .m files and any written/typed work as a .pdf. Submit code and
written parts to Moodle.
(1) The e
or for this model is calculated in the same way we’ve used above:
E =
nX
i=1
(population(i)� modelPopulation(i))2
Find the e
or of the model used in the MATLAB code. We see this e
o
is a lot larger than e
ors with previous models (like the scaled temperature
model). Why does this e
or seem so much larger than previous e
ors? Can
we judge a model based purely on its e
or? What considerations would we
need to take into account when judging the ‘goodness’ of a model based on
its e
or?
(2) Use the model to predict the US population in XXXXXXXXXXyears after 1860).