The solution request is only for the question entitled "Groupwork (population example) on page 7) in the pdf attached. I set a short deadline period because it is just one question.
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 arrive 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 error. The error 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, or a ‘full’ derivative? This gives us a nonlinear algebraic 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. 36 -1.0000 18 Feb. 39 -0.8796 21 March 48 -0.5184 29 April 60 -0.0368 39 May 71 0.4047 49 June 79 0.7258 58 July 84 0.9264 64 Aug. 83 0.8863 63 Sept. 75 0.5652 54 Oct. 63 0.0836 42 Nov. 52 -0.3579 34 Dec. 41 -0.7993 24 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, ... 4 0.9264, 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 for 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 carrying 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 error. • 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 be for each? Would an upper bound be useful to specify? Year 1860 1870 1880 1890 1900 1910 1920 1930 Population 31.4 35.6 50.2 63.0 76.2 92.2 106.0 123.2 Year 1940 1950 1960 1970 1980 1990 2000 2010 Population 132.2 151.3 179.3 203.2 226.5 248.7 281.4 308.7 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, r 4 % p(2): carrying 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 lowerbound=[0, 0, 0]; % lower bound for parameters in initGuess vector 10 options=optimset('TolFun',1e�8,'Display','off'); % setting tolerance for solving, and turning output off 11 12 [p, error] = lsqcurvefit(@logistic,initGuess,decades,population, lowerbound,[],options); % getting r, K, p0 coefficients that minimize error 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 error for this model is calculated in the same way we’ve used above: E = nX i=1 (population(i)� modelPopulation(i))2 Find the error of the model used in the MATLAB code. We see this error is a lot larger than errors with previous models (like the scaled temperature model). Why does this error seem so much larger than previous errors? Can we judge a model based purely on its error? What considerations would we need to take into account when judging the ‘goodness’ of a model based on its error? (2) Use the model to predict the US population in 2020 (160 years after 1860).