clear
reset(0)
echo on
%This example demonstrates the use of matrix
%solution techniques for solving linear regression
%problems.  As a simple example, suppose we have
%measured the concentration of a reservoir as a
%function of time in order to determine the mass
%transfer coefficient of a membrane.  For the
%exact geometry, see the notes.  We have the
%observed concentrations and times:
t=[0,1,2,4,8];
%and
c=[2.1,1.55,1.38,1.13,1.016];
pause
%We can plot this up:
figure(1)
plot(t,c,'o')
set(gca,'FontSize',18)
xlabel('time')
ylabel('concentration')
hold on
%Note that the concentration approaches an equilibrium
%value at large times.  We set things so that the
%equilibrium value is 1.0.
pause
%In order to use linear regression, our equation
%must be linear in the modelling parameters (not
%necessarily in time, however).  We accomplish this
%by transforming the problem:
%   (c-ceq)=(c0 - ceq)exp(-hAt/V)
%becomes:
%   ln(c-ceq)=ln(c0-ceq) - (hA/V)t
%where ln(c0-ceq) and -(hA/V) are the new modelling
%parameters.  We can solve for this directly:
pause
%Letting x2=-(hA/V) and x1=ln(c0-ceq).  We get:
n=length(t);
ceq=1.0;
y=log(c-ceq);
x2=(y*t'-sum(y)*sum(t)/n)/(t*t'-sum(t)^2/n)
x1=(sum(y)-x2*sum(t))/n
pause
%We can plot up this result:
plot(t,ceq+exp(x1+x2*t))
%Which more or less goes through the data points.
hold off
pause
%There is more than one way to solve the system
%of equations.  The regression problem we are
%solving can be turned into a set of linear equations.
%We define the -modelling functions- as a matrix a:
a=[ones(n,1),t']
%and we define the right hand side (the measured
%data points that we are fitting the model to) by
%the column vector b:
b=y'
%The solution vector of modeling parameters x obtained
%for linear regression can be found from the solution
%to the problem (a'*a)*x = a'*b.  These are known as the
%normal equations for linear regression.
pause
%If we solve this problem using the matlab equation
%solver, we get:
x=a\b
%which is the same as the direct solution approach.  Note
%that the matlab operator "\" solves this problem in the
%least squares sense.
pause
%We can also do this explicitly by using the normal equations:
k=inv(a'*a)*a';
x=k*b
%which again yields the same result.  As we shall see in the
%next class, this approach has advantages, because it shows
%that the modelling parameters are just a linear combination
%of the observations - and thus it is simple to figure out
%the matrix of covariance (error) of these parameters from
%the error in the data!
pause
%As a final note, what do you do if you don't know ceq?
%In this case, you can no longer linearize the model, and
%thus you can't use linear regression.  You still can solve
%it, though: all you have to do is define an objective
%function as the square of the residual for the non-linear
%model, and then use a non-linear minimization routine such 
%as fminsearch to find the optimum parameters.  We'll look
%at such a problem in the last lecture.
echo off