%
%   FXY_FIT.M   Simple example of multiple regression  --  that is, fitting a curve
%                  to a function of two variables, say f(x,y)
%
%   For this example, we assume that the function f varies versus both x and y.  A 
%   small sample of data are given (see below) that has 9 sets of triplets, fi, xi, 
%   and yi, where the xi and yi values are the independent variables and fi is the 
%   measured value of f(x,y) for the given xi, yi pair.
%
%   In this exercise we assume a mathematical model that has a simple linear 
%   dependence on both x and y
%         f(x,y) = a1 + a2*x + a3*y
%   and the goal is to find the coeffs a1, a2, and a3 that gives the best match of 
%   the model to the actual experimental data.  The normal equations for this 
%   situation are set up and solved and the coeff of determination is computed for 
%   the given data set.
%
%   Program written by J. R. White, UMass-Lowell (Dec. 2013)
%

      clear all,   close all,   nfig = 0;
%
%   base data
      x = [0 1.0 1.1 2.1 2.3 2.9 3.1 4.3 3.7];  N = length(x);
      y = [0 1.0 2.2 1.3 1.8 1.4 2.1 1.2 2.6];  
      f = [15.1 17.9 12.7 25.6 20.5 35.1 29.7 45.4 40.2];
%
%   generate the P and q matrices for the normal eqns.:  Pa = q
      P = [  N       sum(x)      sum(y)   ;
           sum(x)   sum(x.*x)   sum(x.*y) ;
           sum(y)   sum(x.*y)   sum(y.*y)];
      q = [sum(f)   sum(x.*f)   sum(y.*f)]';
%
%   solve the normal eqns and define anonymous function fxy
      a = P\q;     fxy = @(x,y) a(1) + a(2)*x + a(3)*y;
%
%   do error analysis for model A
      fm = sum(f)/N;           % mean value of f
      St = sum((f-fm).^2);     % sum of squares of deviation (original data)
      fe = fxy(x,y);           % evaluate math model at expt points
      Sr = sum((f - fe).^2);   % sum of squares of residuals (model - org data)
      rsq = (St - Sr)/St;      % coeff of determination (approx unity for good model)
%
%   plot data (use more points for smooth curve)
      xx = linspace(0,5,51);       yy = linspace(0,3,31);
      [XX,YY] = meshgrid(xx,yy);   FF = fxy(XX,YY);
%   default view
      nfig = nfig+1;   figure(nfig)
      plot3(x,y,f,'ro','LineWidth',2),grid, hold on
      surf(XX,YY,FF),colormap(white)
      title(['FXY\_FIT:  Model Prediction Versus Experimental Data (r^2 = ',...
             num2str(rsq,'%5.3f'),')'])
      xlabel('x values'),ylabel('y values'),zlabel('f(x,y)')
      legend('Experimental Data', 'Math Model','Location','Best')
%   side view
      nfig = nfig+1;   figure(nfig)
      plot3(x,y,f,'ro','LineWidth',2),grid, hold on
      surf(XX,YY,FF),colormap(white),view(-9,8)
      title(['FXY\_FIT:  Model Prediction Versus Experimental Data (r^2 = ',...
             num2str(rsq,'%5.3f'),')'])
      xlabel('x values'),ylabel('y values'),zlabel('f(x,y)')
      legend('Experimental Data', 'Math Model','Location','Best')
%
%   Note: you can rotate either of these views to see how the measured data are
%         distributed on both sides of the plane defined by the math model
%
%   end of problem 
