clear
echo on
format compact
%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')

%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:
figure(1)
hold on
plot(t,ceq+exp(x1+x2*t),'r')
legend('data','model')
hold off
%Which more or less goes through the data points.
pause
%You will note that the fitted curve isn't very smooth!
%This is because we evaluated it only at a the same points
%were we measured the data points.  We can do much better
%than this, however!  Let's define a vector tm:
tm=[min(t):.1:max(t)]';
%and plot up the smooth curve:
pause
figure(1)
hold on
plot(tm,ceq+exp(x1+x2*tm),'k')
legend('data','model','smooth model')
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'*a)\(a'*b)
%which is the same as the direct solution approach.
%The advantage of this approach is that we can easily
%extend it to any number of modelling parameters.  Doing
%it this way, the "\" operator would use PLU decomposition
%(Gaussian elimination) for this square problem.
pause
%We can also directly solve it via the \ operator, which uses QR
%for this overdetermined problem:
x=a\b
%which again is the same result.
pause
%And, finally, we can define a matrix "k" from the normal
%equations and use it to solve the problem:
k=inv(a'*a)*a'
x=k*b
%This is useful, because the matrix k shows how each of
%the fitting parameters depends on the elements of the
%data set (they aren't weighted uniformly!).  You also only
%need to compute "k" once for many possible vectors of b.
%This approach will also be very useful when we do error 
%calculations!
echo off