Fitting Multiple-Response Data

Bruce Miller
brucem@wolfram.com
Last edit: November, 1999 

Mathematica has a built-in function Fit and a Standard AddOn Package  Statistics`NonlinearFit` for  finding the least-squares fit of an expression with adjustable parameters to data.  They fit only one expression at a time, however, and it is sometimes convenient to fit more than one expression to a set of data, with the parameters being adjusted to minimize the combined error in all the expressions being fit.  This package includes a fairly general function  for performing such fits.  

The example below shows the basic technique.  Two functions, y1[x1,x2] and y2[x1,x2], are to be fit to {x1,x2,{y1,y2}} data using the same parameters for the two functions.  This is most easily done by building an expression for the total error in both functions in terms of the parameters,  then finding the parameter values which minimize that total error.   This is a general technique  used for fitting irregular functions or data.  (One example from the Technical Support FAQ site is  in http://support.wolfram.com/Math/Statistics/NonlinearOptimization.html.)

Acknowledgement

This package and accompanying notebook are significantly expanded versions of a notebook published in The Mathematica Journal, Volume 7, Issue 3, as part of the Technical Support column.

Create sample data

Make some fake data (with a little noise).     I need the line
         dat2 = Flatten[dat,1];
because the Table command makes a matrix of points, giving an extra level of braces ({}) around rows of points.  The rest of my example needs the data in {...,point,point,..} form.  The  Flatten eliminates the division of the matrix into rows of points and give just one long list of {x1,x2,{y1,y2}} points.

dat = Table[{x1,x2, 
              {1  + 2.2*x1 + 3.3*x2^2 + Random[Real,{-0.01,0.01}],
               Sqrt[1+1.1*x1+2.2*x2] + Random[Real,{-0.01,0.01}]}
            },
            {x1, 0, 4},{x2,0,3}];

dat2 = Flatten[dat,1];

The first points in the data list are

Take[dat2, 2]

The basic technique

Define the two functions describing the expressions to fit.  The function arguments are the free parameters to be fitted and the independent variables.

y1[a_,b_,c_,x1_,x2_]:= a + b*x1 + c*x2^2;
y2[d_,f_,x1_,x2_]:= Sqrt[1+d*x1+f*x2];

Define an error function which returns chi-squared at this point given these parameter values.  Run this function over all the data points, add the individual point chi-squared errors together, and find the parameters which minimize the total.

erfunc[a_,b_,c_,datum_]:=
             (y1[a,b,c,datum[[1]],datum[[2]]]-datum[[3,1]])^2 +
             (y2[a,b,datum[[1]],datum[[2]]]-datum[[3,2]])^2

chisq = Map[erfunc[a,b,c,#]&,dat2];

The Apply[Plus[..  in  FindMinimum below turns the list of squared errors   chisq  into a sum, that is, 
         {error^2,error^2,error^2,..} => error^2+error^2+error^2+..
If the FindMinimum complains about Complex values, put the first argument inside an Abs  or Chop,
 
             FindMinimum[ Abs[Apply[Plus,chisq]],..
 
This will eliminate tiny imaginary components which sometimes arise in finite-precision arithmetic.

soln = FindMinimum[Apply[Plus,chisq],{a,1},{b,1},{c,1}]

The output from FindMinimum is a list containing the final value of the function being minimized and a sublist of replacement rules for the fitted parameters.  We can use the result to create final fitted functions for calculations or plotting.

y1[x1_,x2_] = y1[a,b,c,x1,x2] /. soln[[2]]

y2[x1_,x2_] = y2[a,b,x1,x2] /. soln[[2]]

Plot3D[y1[x1,x2],{x1,0,4},{x2,0,3},Lighting->False,
                     ColorFunction->Hue];

The mrfit package

To make the above procedure easier to use, the algorithm can be placed into a function so it can be called with different data and expressions.  The  job performed by the errfunc function above is in the separate function onept (one-point).  The mrfit (Multiple Response Fit) function returns a list of rules for the fitted parameters.  The file MRFit.m has the code and additional comments.

<< MRFit.m

?mrfit

This function takes two options, in addition to the options it can pass on to FindMinimum.

Options[mrfit]

?mrfitErrorFunction

?mrfitReport

For convenience, I put the list of expressions to fit into a symbol.

functofit={a+b*x1+c*x2^2, Sqrt[1+a*x1+b*x2]};

Call the function.   Note that the parameters need initial values.

solu=mrfit[dat2,functofit,
         {x1,x2},{{a,1},{b,1},{c,1}}]

Symbolic functions can be created easily from these fitted values.  The  mrfitReport option can be used to return the expression directly.

{f[x1_,x2_],g[x1_,x2_]}=functofit /. solu

f[2,1.5]

mrfit[dat2, 
      functofit, {x1, 
        x2}, {{a, 1}, {b, 1}, {c, 
          1}}, 
   mrfitReport -> {mrfitFittedModels}]

This routine can also be used to fit one expression at a time to {x1,..,y} or {x1,..,{y}}  data.  Here I also show a different format for the parameters (useful if the models cannot be symbolically differentiated) and two options which are passed to FindMinimum.  When two starting values are given,  the Method option is needed - either Method->Gradient for Mathematica 4   or  Method->FindMinimum  for Mathematica 3.

dat4=Table[{x,BesselJ[0,x]},{x,-2,2}];
mrfit[dat4,{a*x^2+b*x+c},
         {x},{{a,{0,2}},{b,{0,2}},{c,{-2,2}}},
         Method->Gradient,AccuracyGoal->3]

Suggestions and caveats

The mrfit function wants the models, independent variables, and parameter specifications in lists, even if there is only one.  This means that if there is only one parameter to fit, it must be doubly nested (i.e., {{a,1}}). 

If you want a quick estimate of the goodness-of-fit, use the   mrfitReport option to print the final value of the sum of squared errors (the "residual" error) along with the parameter values.  Remember to extract the rules from the combined output before using them.

mrfre=mrfit[dat4,{a*x^2+b*x+c},
         {x},{{a,1},{b,0},{c,1}},
         mrfitReport ->
           {mrfitFittedParameters,mrfitResidualError}]

{afit,bfit,cfit} = {a,b,c} /. (mrfitFittedParameters /. mrfre )

To use these functions with Complex data and models,  call mrfitc or call mrfit with the option mrfitErrorFunction->oneptc.  This will necessitate specifying two starting values for each parameter, as in the first dat4 example above.  

If you have a different way to calculate the error (for small-N Poisson statistics, for example), just write your own version of the onept function and use the mrfitErrorFunction to pass its name to  mrfit. 

I deleted the General::spell1 messages from the code examples above. 

Nonlinear fitting codes seek minima in the total chi-squared error function, not distinguishing between local minima and global minima.  Careful selection of parameter starting values and model functions may be necessary for satisfactory results.

Additional tests and demos for MRFit.m package

Create sample data

dat2f = Flatten[
            Table[
            {x1,x2, x3,
               { 1 + 2.2*x1 + 3.3*x2^2 + 1.5*Sin[1.5 x3] +
                                        Random[Real,{-0.01,0.01}],
                Sqrt[1+1.1*x1+2.2*x2] + Random[Real,{-0.01,0.01}]}
             },
             {x1,0,4},{x2,0,3},{x3,0,4}],
          2];

dat3f = Flatten[
            Table[{x1,x2, x3,
              { 1 + 2.2*x1 + 3.3*x2^2 + 1.5*Sin[1.5 x3] +
                                        Random[Real,{-0.01,0.01}],
               Sqrt[1+1.1*x1+2.2*x2] + Random[Real,{-0.01,0.01}],
               2.2*x1*Exp[-3.3*x3]}
             },
             {x1,0,4},{x2,0,3},{x3,0,4}],
         2];

The first points in the data list dat2f are

Take[dat2f, 2]

The 32nd point in the dat3f list:

Part[dat3f,32]

Two- and three-model fits

For convenience, I put the lists of expressions to fit into symbols.

functofit2 = {1 + a*x1 + b*x2^2 + c*Sin[d x3] ,
                Sqrt[1+e*x1+a*x2] };
 
functofit3 = {1 + 2.2*x1 + b*x2^2 + c*Sin[d x3] ,
               Sqrt[1+e*x1+a*x2] ,
               a*x1*Exp[-b*x3] };

Call the function.   Note that the parameters need initial values.

solu2=mrfit[dat2f,functofit2,
         {x1,x2,x3},{{a,1},{b,1},{c,1},{d,1},{e,1}}]

solu3=mrfit[dat3f,functofit3,
         {x1,x2,x3},{{a,1},{b,1},{c,1},{d,1},{e,1}}]

Symbolic functions can be created easily from these fitted values.  The  mrfitReport option can be used to return the expression directly.

\(mrfit[dat3f, 
      functofit3, 
\ \ \ \ \ \ \ \ \ {x1, x2, 
        x3}, {{a, 1}, {b, 1}, {c, 1}, {d, 1}, {e, 
          1}}, 
\ \ \ mrfitReport -> {mrfitFittedModels, mrfitResidualError}]\)

Redoing the example with simple {x, y} data, this time comparing original and fitted curves and demonstrating options.  In the Plot, the solid red curve is the original Bessel function and the dashed blue-green curve is the parabola fit to five points.  The points are drawn with a Prolog option so the curves will overlay them.

beslfit=mrfit[dat4,{a*x^2+b*x+c},
         {x},{{a,{0,2}},{b,{0,2}},{c,{-2,2}}},
         Method->Gradient,AccuracyGoal->3,
         mrfitReport->True]

Plot[{BesselJ[0,x],mrfitFittedModels /. beslfit},{x,-2.2,2.2},
      PlotStyle->{{Hue[0]},{Hue[0.6],Dashing[{0.02}]}},
      Prolog->{PointSize[0.03],Map[Point,dat4]} ];