Nonlinear Regression

Problem 15.11.

The following model is used to represent the effect of solar radiation on the photosynthesis rate of aquatic plants:

where P = the photosynthesis rate (mg m^-3d^-1), P_m = the maximum photosynthesis rate (mg m^-3d^-1), I = solar radiation (µE m^-2s^-1), and I_sat = optimal solar radiation (µE m^-2s^-1). Use nonlinear regression to evaluate P_m and I_sat based on the following data:

Solution:

First, a M-file function must be created to compute the sum of squares: 

function f = fSSR(a,Im,Pm)

Pp = a(1)*Im/a(2).*exp(-Im/a(2)+1);

% a(1) = Pm, a(2) = Isat

f = sum((Pm-Pp).^2);

In command mode, the data can be entered as:

>> I = [50 80 130 200 250 350 450 550 700];

>> P = [99 177 202 248 229 219 173 142 72];

Employ initial guesses of 200 for the coefficients. The minimization of the function is then implemented by:

>> a = fminsearch(@fSSR, [200, 200], [], I, P)

a =

  238.7124  221.8239

The best-fit model is therefore

The fit along with the data can be displayed graphically as:

>> Ip = linspace(min(I),max(I));

>> Pp = a(1)*Ip/a(2).*exp(-Ip/a(2)+1);

>> plot(I,P,'o',Ip,Pp)

>> xlabel('I'),ylabel('P')

>> title('Effect of solar radiation on photosynthesis rate')

General Linear Least-Squares Regression

Problem 15.10.

Three disease-carrying organisms decay exponentially in a seawater according to the following model:

p(t) = Ae^-1.5t + Be^-0.3t + Ce^-0.05t 

Use general linear least-squares to estimate the initial concentration of each organism (A, B and C) given the following measurements:

Solution:

Firstly, note that the given model is actually a multiple linear function, where e^-1.5t, e^-0.3t and e^-0.05t are x₁, x₂ and x₃ respectively. And there is no intercept. Let's solve the problem using MATLAB!

The data can be entered as column vectors and the [Z] matrix can be set up as below:

>> t = [0.5 1 2 3 4 5 6 7 9]';

>> p = [6 4.4 3.2 2.7 2 1.9 1.7 1.4 1.1]';

>> x1 = exp(-1.5*t); x2 = exp(-0.3*t); x3 = exp(-0.05*t);

>> Z = [x1 x2 x3]

Z =

    0.4724    0.8607    0.9753

    0.2231    0.7408    0.9512

    0.0498    0.5488    0.9048

    0.0111    0.4066    0.8607

    0.0025    0.3012    0.8187

    0.0006    0.2231    0.7788

    0.0001    0.1653    0.7408

    0.0000    0.1225    0.7047

    0.0000    0.0672    0.6376

Then, the coefficients can be generated using Eq.(15.10) in textbook:

>> a = (Z'*Z)\(Z'*p)

a =

    4.1375 = A

    2.8959 = B

    1.5349 = C

Thus, the least-squares fit is p(t) = 4.1375e^-1.5t + 2.8959e^-0.3t + 1.5349e^-0.05t .

Best-fit line and the data can be plotted as

>> pp = Z*a;

>> plot(t,p,'bo',t,pp,'k-')

>> xlabel('t'),ylabel('p'),title('Decay Exponential Model')

General Linear Least-Squares Regression

Here are the M-files to implement quadratic (second-order polynomial) regression, multiple linear regression and general linear least-squares regression.

Quadratic Regression M-file

function [a,r2,syx] = polyregr(order,x,y)

% polyregr: least-squares regression curve fitting (quadratic)

%   [a, r2] = polyregr(x,y): Least squares fit of curve line to data by solving the normal equations

% input:

%   x = independent variable

%   y = dependent variable

% output:

%   a = vector of coefficients->> a(1):intercept,a(2):slope x,a(3):slope x2

%   r2 = coefficient of determination

%   syx = standard error of estimate

n = length(x);

if length(y)~=n, error('x and y must be same length'); end

% error checking to have complete quantity of data point

x = x(:); y = y(:); % convert to column vectors

sx = sum(x); sy = sum(y);

sx2 = sum(x.^2); sxy = sum(x.*y);

sx3 = sum(x.^3); sx4 = sum(x.^4); sx2y = sum(x.*x.*y);

matrix = [n sx sx2; sx sx2 sx3; sx2 sx3 sx4];

right = [sy sxy sx2y]';

a = matrix \ right;

% to calculate error statistics

Sr = sum((y-a(1)-(a(2)*x)-(a(3)*x.^2)).^2);

St = sum((y-(sy/n)).^2);

r2 = (St-Sr)/St;

syx = sqrt(Sr/(length(x)-length(a)));

% create plot of data and best fit line

xp = min(x):0.1:max(x);

yp = zeros(1,length(xp));

    for k = 1:length(yp) 

        for l = 0:order

            yp(k) = yp(k) + a(l+1)*(xp(k))^l;

        end

    end

plot(x,y,'o',xp,yp),xlabel('T (K)'),ylabel('cp (J mol-1 K-1')

grid on

Multiple Linear Regression M-file

function [a,Sr,r2,syx] = multilinear(x1,x2,y)

% multilinear: multiple linear regression curve fitting

% input:

%   x1 = independent variable 1

%   x2 = indepencent variable 2

%   y  = dependent variable (function of x1 and x2)

% output:

%   a   = vector of coefficient->> a(1)= intercept,a(2)= slope x1,a(3)= slope x2

%   Sr  = sum of square of residuals

%   r2  = coefficient of determination

%   syx = standard error of estimate

n = length(y);

if length(x1)~=n, error('x1 and y must be same length'); end

if length(x2)~=n, error('x2 and y must be same length'); end

x1 = x1(:); x2 = x2(:); y = y(:); % convert to column vectors

sx1 = sum(x1);sx2 = sum(x2);sy = sum(y);sx1x2 = sum(x1.*x2);

sx1y = sum(x1.*y); sx2y = sum(x2.*y);

sx1sq = sum(x1.*x1); sx2sq = sum(x2.*x2);

%calculate for coefficients [a]

A = [n sx1 sx2;sx1 sx1sq sx1x2;sx2 sx1x2 sx2sq];

b = [sy;sx1y;sx2y];

a = A\b;

%calculate for Sr

yp = a(1)+a(2)*x1+a(3)*x2;

Sr = sum((y-yp).^2);

%calculate for syx

syx = sqrt(Sr/(n-length(a)));

%calculate for r2

r2 = 1-(Sr/sum((y-mean(y)).^2));

end

General Linear Regression M-file

function [a, r2, syx] = genlinreg(x1,x2,y)

% genlinreg: general linear least-squares regression model 

%   [a, r2, syx] = genlinreg(x,y): Least squares fit of line to data by solving the normal equations

% input:

%   x1 = independent variable 1

%   x2 = independent variable 2 = zeros(size(x1)) for polynomial regression

%   y  = dependent variable

% output:

%   a   = vector of coefficients->> a(1):intercept,a(2):slope for x1/x,a(2):slope of x2/x^2

%   r2  = coefficient of determination

%   syx = standard error of estimate

n = length(y);

if length(x1)~=n, error('x1 and y must be same length'); end

if length(x2)~=n, error('x2 and y must be same length'); end

x1 = x1(:); x2 = x2(:); y = y(:); % convert to column vectors

if x2 == zeros(size(x1))

    Z = [ones(size(x1)) x1 x1.^2]; 

else 

    Z = [ones(size(y)) x1 x2];

end 

a = (Z'*Z)\(Z'*y); %solve for the unknown coefficients by utilizing backslash operator 

% calculate for coefficient of determination

Sr = sum((y-Z*a).^2);

St = sum((y-mean(y)).^2);

r2 = 1-Sr/St;

% calculate for standard error of estimate

syx = sqrt(Sr/(n-length(a)));

end

General Linear Least-Squares

Problem 15.6.

Use multiple linear regression to derive a predictive equation for dissolved oxygen concentration as a function of temperature and chloride based on the data from Table P15.5. Use the equation to estimate the concentration of dissolved oxygen for a chloride concentration of 15 g/L at T = 12 °C. Note that the true value is 9.09 mg/L. Compute the percent relative error for your prediction. Explain possible causes for the discrepancy. 

Table P15.5 Dissolved oxygen concentration in water as a function of temperature (°C) and chloride concentration (g/L)

Solution:

The multiple linear regression to evaluate is 

*y is dissolved oxygen concentration (mg/L).

The [Z] and y matrices can be set up using MATLAB commands:

>> %enter the data T, c & y to be fit

>> T = [0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30]';

>> c = [0 0 0 0 0 0 0 10 10 10 10 10 10 10 20 20 20 20 20 20 20]';

>> y = [14.6 12.8 11.3 10.1 9.09 8.26 7.56 12.9 11.3 10.1 9.03 8.17 7.46 6.85 11.4 10.3 8.96 8.08 7.35 6.73 6.20]';

>> %create the [Z] matrix

>> Z = [ones(size(T)) T c];

>> %solve for the coefficients of the least-squares fit

>> a = (Z'*Z)\(Z'*y)

a =

  13.522142857142857

  -0.201238095238095

  -0.104928571428572

Thus, the best fit multiple regression model is

We can evaluate the prediction at T = 12°C and c = 15 g/L and evaluate the percent relative error as

>> yp = a(1)+a(2)*12+a(3)*15

yp =

   9.533357142857142

>> ea = abs((9.09-yp)/9.09)*100

ea =

   4.877416313059866

Thus, the error is considerable. This can be seen even better by generating predictions for all the data and then generating a plot of the predictions versus the measured values. A one-to-one line is included to show how the predictions diverge from a perfect fit.

>> yp = a(1)+a(2).*T+a(3).*c;

>> ymin = min(y);ymax = max(y);

>> dy = (ymax - ymin)/100;

>> ymod = [ymin:dy:ymax];

>> plot(y,yp,'ko',ymod,ymod,'k-')

>> axis square,title('Plot of predicted vs measured values of dissolved oxygen concentration')

>> legend('model prediction','1:1 line','Location','Northwest')

>> xlabel('y measured'),ylabel('y predicted')

The cause for the discrepancy is because the dependence of oxygen concentration on the unknowns is significantly nonlinear. It should be noted that this is particularly the case for the dependency on temperature.

Gauss Elimination

Here are the M-files to implement naive Gauss Elimination and Gauss Elimination with partial pivoting.

Naive Gauss Elimination M-file

function x = GaussNaive(A,b)

% GaussNaive: naive Gauss elimination

%   x = GaussNaive(A,b): Gauss elimination without pivoting.

% input:

%   A = coefficient matrix

%   b = right hand side vector

% output:

%   x = solution vector

[m,n] = size(A);

% m = no of rows, n = no of columns, size of matrix A is m x n

if m~=n, error('Matrix A must be square'); end

% if m not equal to n, error message is displayed as 'Matrix A must be square'

nb = n+1;

% no of columns for augmented matrix = no of columns of the matrix A + 1 column for the vector B

Aug = [A b];

% Coefficient matrix A and right-hand-side vector b are combined in an augmented matrix

% forward elimination

for k = 1:n-1

    % k represents column 1 to (n-1)th column

  for i = k+1:n

      % i represents row 2 to nth row

    factor = Aug(i,k)/Aug(k,k);

    % to compute the factor in order to obtain the upper triangular matrix

    Aug(i,k:nb) = Aug(i,k:nb)-factor*Aug(k,k:nb);

    % to compute the new equations until transform the matrix into upper triangular form

    % the new ith row kth column elements is equal to the old ith row kth column elements 

    % minus the multiplication of factor and kth row kth column elements 

    % since using augmented matrix, need to use nb (n+1) because of additional 1 column vector b

  end

end

% back substitution

x = zeros(n,1);

% solution vector x is a matrix of zero coefficients of n x 1 dimension

% n =  no of equations

x(n) = Aug(n,nb)/Aug(n,n);

% to calculate the last unknown x in nth equation

% by dividing the constant b in last equation with the coefficient of x(n)

% formula is same as Eq (9.12) in textbook

for i = n-1:-1:1

    % i represents no of equations (n) minus 1 with increment of -1 until 1

  x(i) = (Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);

  % to calculate the remaining unknown x's in (n-1)th equation and so on until the first equation

  % formula is same as Eq (9.13) in textbook

end

Gauss Elimination with Partial Pivoting M-file

function x = GaussPivot(A,b)

% GaussPivot: Gauss elimination pivoting

%   x = GaussPivot(A,b): Gauss elimination with pivoting.

% input:

%   A = coefficient matrix

%   b = right hand side vector

% output:

%   x = solution vector

[m,n]=size(A);

% m = no of rows, n = no of columns,size of matrix A is m x n 

if m~=n, error('Matrix A must be square'); end

% if m not equal to n, error message will be displayed as 'Matrix A must be square'

nb=n+1;

% no of columns of augmented matrix = no of columns of matrix A + 1 column of vector b 

Aug=[A b];

% forward elimination

for k = 1:n-1

    % k represents column 1 to (n-1)th column

  % partial pivoting

  [big,i]=max(abs(Aug(k:n,k)));

  % to determine the coefficient with the largest absolute value in the column below the pivot element

  % big = largest element

  % i = index corresponding to the largest element (no of row)

  ipr=i+k-1;

  % new row (with largest element below the pivot element)= old row (with pivot element)+ k - 1

  if ipr~=k

      % if new row doesn't equal to k (row with largest element isn't row 1)

    Aug([k,ipr],:)=Aug([ipr,k],:);

    % to switch the row with largest pivot element to the first row

  end

  for i = k+1:n

      % i represents row 2 to nth row

    factor=Aug(i,k)/Aug(k,k);

    % to compute the factor in order to obtain the upper triangular matrix

    Aug(i,k:nb)=Aug(i,k:nb)-factor*Aug(k,k:nb);

    % to compute the new equations until transform the matrix into upper triangular form

    % the new ith row kth column elements is equal to the old ith row kth column elements 

    % minus the multiplication of factor and kth row kth column elements 

    % since using augmented matrix, need to use nb (n+1) because of additional 1 column vector b

  end

end

% back substitution

x=zeros(n,1);

% solution vector x is a matrix of zero coefficients of n x 1 dimension

% n =  no of equations

x(n)=Aug(n,nb)/Aug(n,n);

% to calculate the last unknown x in nth equation

% by dividing the constant b in last equation with the coefficient of x(n)

% formula is same as Eq (9.12) in textbook

for i = n-1:-1:1

    % i represents no of equations (n) minus 1 with increment of -1 until 1

  x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);

  % to calculate the remaining unknown x's in (n-1)th equation and so on until the first equation

  % formula is same as Eq (9.13) in textbook

end

Linear Least-Squares Regression

Problem 14.14.

An investigator has reported the data tabulated below for an experiment to determine the growth rate of bacteria k (per d) as a function of oxygen concentration c (mg/L). It is known that such data can be modeled by the following equation: 

 Eq.(14.14(a))

where c_s and k_max are parameters. Use a transformation to linearize this equation. Then use linear regression to estimate c_s and k_max and predict the growth rate at c = 2 mg/L.

Solution:

Equation (14.14(a)), which is in the format of saturation-growth-rate model can be linearized by inverting it to give:

Thus, a plot of 1/k (y) versus 1/c² (x) will yield a straight line with a slope of c_s/k_max and an intercept of 1/k_max. The data can be set up in tabular form and the necessary sums computed as in table below.

The means can be computed as:

The slope and intercept can then be calculated with Eqs (14.15) and (14.16) as in textbook.

The least-squares fit of the transformed data is 1/k = 0.09666 + 0.20201(1/c²).

The model coefficients can then be calculated as:

k_max = 1/0.09666 = 10.34554

c_s = 0.20201 x 10.34554 = 2.0899

The growth rate at c = 2 mg/L can be predicted as:

1/k = 0.09666 + 0.20201(1/4) = 0.1471625

k = 6.7952 per day

*The linregr M-file can also be used to determine the least-squares fit:

>> format long

>> c = [0.5 0.8 1.5 2.5 4];

>> k = [1.1 2.5 5.3 7.6 8.9];

>> [a,r2,syx] = linregr(1./c.^2,1./k)

a =

   0.202005196520525   0.096665700618091

r2 =

   0.999569288346182

syx =

   0.007995513530533

The model coefficients can then be calculated as:

>> kmax = 1/a(2)

kmax =

  10.344930969370603

>> cs = kmax*a(1)

cs =

   2.089729813458969

Linear Least-Squares Regression

Here is an M-file to implement linear regression. 

Linear Regression M-file

function [a,r2,syx] = linregr(x,y)

% linregr: linear regression curve fitting

%   [a, r2] = linregr(x,y): Least squares fit of straight line to data by solving the normal equations

% input:

%   x = independent variable

%   y = dependent variable

% output:

%   a = vector of slope, a(1), and intercept, a(2)

%   r2 = coefficient of determination

%   syx = standard error of estimate

n = length(x); % show number of data

if length(y)~=n, error('x and y must be same length'); end

% error checking to have complete number of data points

x = x(:); y = y(:); % convert to column vectors

sx = sum(x); sy = sum(y);

sx2 = sum(x.*x); sxy = sum(x.*y); sy2 = sum(y.*y);

a(1) = (n*sxy-sx*sy)/(n*sx2-sx^2); % a(1) = slope

a(2) = sy/n-a(1)*sx/n; % a(2) = intercept

r2 = ((n*sxy-sx*sy)/sqrt(n*sx2-sx^2)/sqrt(n*sy2-sy^2))^2;

Sr = sum((y-a(2)-a(1)*x).^2);

St = sum((y-mean(y)).^2);

syx = sqrt(Sr/(n-2)); 

% create plot of data and best fit line

xp = linspace(min(x),max(x),2);

yp = a(1)*xp+a(2);

% best fit line equation

plot(x,y,'o',xp,yp),xlabel('x'),ylabel('y')

% to generate plot of x vs y with circle marker and plot for best fit line

grid on