Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have some data from a study that I am trying to fit with a rational function of arbitrary form:

$y(x) = \frac{\sum_{i=n}^m A_ix^{i/q}}{\sum_{j=o}^p B_jx^{j/r}}$

where $i$, $j$, $k$, $l$, $m$, $n$, $o$, $p$, and $q$ are all integers. The goal is to produce a function of this form with as few terms as possible that still provides a reasonable match to the data. For instance, I have some data that looks like this:

enter image description here

where the black line is a fit to the data, and the fit is a function of the following form:

$y(x) = \frac{-0.31+1.2 \sqrt{x}-1.2 x }{1.0+1.3 \sqrt{x}-4.2 x}$

My current method of generating these functions is to produce a rational function of a high order, and then to successively eliminate terms if their coefficients are small compared to other terms in the expression. My algorithm is pretty ad-hoc, and the cutoffs, maximum order to use, and fractional powers to include are all parameters that one can pass as inputs to the function. I was curious if anyone had a better way to produce a compact rational fit to an arbitrary set of data that perhaps didn't rely on such an ad-hoc approach.

Here is the code I currently use to generate these fits:

compactRationalFit[data_, x_, coeffthresh_: 0.1, minorder_: 0, 
   maxorder_: 4, powerfrac_: 1] := 
  Module[{lcoeffthresh, bestfit, cld, cln, csc, i, dd, nn, norderlist,
     dorderlist, ndelpos, ddelpos, result, totdel, norigdelems, 
    norignelems},
   dorderlist = norderlist = Range[minorder, maxorder, 1/powerfrac];
   norignelems = Length[dorderlist];
   norigdelems = Length[dorderlist];
   totdel = 1;
   cld = ConstantArray[1, Length[dorderlist]];
   cln = ConstantArray[1, Length[norderlist]];
   While[totdel != 
      0 && ! (Length[dorderlist] == 1 && Length[norderlist] == 1),
lminorder = Min[norderlist~Join~dorderlist];
norderlist -= lminorder;
dorderlist -= lminorder;
lncoeffthresh = coeffthresh*Length[norderlist]/norignelems;
ldcoeffthresh = coeffthresh*Length[dorderlist]/norigdelems;
lncoeffthresh = coeffthresh;
ldcoeffthresh = coeffthresh;
cld[[Flatten[Position[cld, x_ /; x < 0]]]] = 1;

bestfit = 
 Normal[NonlinearModelFit[SetPrecision[data, 20],(*{*)
   Total[Table[
       ToExpression[
        "nn" <> 
         ToString[PaddedForm[i, 2, NumberPadding -> "0"]]], {i, 
        Length[norderlist]}] Table[x^i, {i, norderlist}]]/
    Total[Table[
       ToExpression[
        "dd" <> 
         ToString[PaddedForm[i, 2, NumberPadding -> "0"]]], {i, 
        Length[dorderlist]}] Table[x^i, {i, dorderlist}]](*,
   Apply[And,Table[ToExpression["dd"<>ToString[PaddedForm[i,2,
   NumberPadding->"0"]]]>0,{i,Length[dorderlist]}]]}*), 
   Table[{ToExpression[
       "dd" <> ToString[PaddedForm[i, 2, NumberPadding -> "0"]]], 
      cld[[i]]}, {i, Length[dorderlist]}]~Join~
    Table[{ToExpression[
       "nn" <> ToString[PaddedForm[i, 2, NumberPadding -> "0"]]], 
      cln[[i]]}, {i, Length[norderlist]}], x, 
   AccuracyGoal -> 3, PrecisionGoal -> 3, MaxIterations -> 1000, 
   ConfidenceLevel -> .95, 
   VarianceEstimatorFunction -> (Mean[#^2] &)]];
    cld = 
     Table[Coefficient[Expand[Denominator[bestfit]], x, i], {i, 
       dorderlist}]; 
    cln = Table[
      Coefficient[Expand[Numerator[bestfit]], x, i], {i, 
       norderlist}];
    cln = Delete[cln, Position[cln, x_ /; x == 0]]; 
    cld = Delete[cld, Position[cld, x_ /; x == 0]];
    ndelpos = 
     If[Length[cln] == 1, {}, 
      Position[cln, x_ /; Abs[x] < lncoeffthresh Max[Abs[cln]]]];
    ddelpos = 
     If[Length[cld] == 1, {}, 
      Position[cld, x_ /; Abs[x] < ldcoeffthresh Max[Abs[cld]]]];
    totdel = Length[Flatten[ndelpos]] + Length[Flatten[ddelpos]];
    norderlist = Delete[norderlist, ndelpos]; 
    dorderlist = Delete[dorderlist, ddelpos];
    cln = Delete[cln, ndelpos]; cld = Delete[cld, ddelpos];
    csc = First[cld];
    result = 
     SetPrecision[Total[cln/csc*Table[x^i, {i, norderlist}]]/
      Total[cld/csc*Table[x^i, {i, dorderlist}]], 2];
    ];
   Return[result];
   ];

Here, I have made the simplifying assumption that the numerator and the denominator have the same limits/fractional order, i.e. $n=o$, $p=m$, and $q=r$. A call to this function (using the default parameters) looks like this:

compactRationalFit[myData, x]

with the restriction that the number of data points + 1 should not exceed the total number of terms in the numerator and denominator combined.

Some of the other ideas I have been thinking of involve trying all possible combinations simultaneously, and selecting the one with the least amount of terms that satisfies an error threshold, but this involves finding fits to potentially thousands of possible combinations of power series. My current method is relatively fast, as it only requires at most $q^2(m-n)(m-n-1)/2$ iterations, but often times the fits it finds are not perfect, as rational functions are prone to zeroes in the denominator that can produce sharp discontinuities in-between data points:

enter image description here

If anyone could come up with a more robust solution that results in better, more compact fitting formulae, I would be grateful.

share|improve this question
6  
Have a look here. See in particular the section "A mixed discrete−continuous optimization". It does this sort of thing with polynomials. Should not be too hard to extend to rational functions. (There is a more up to date version of this work that will be freely available, but it has not yet appeared.) –  Daniel Lichtblau Jun 8 '12 at 18:31
1  
NonlinearModelFit may not be the best bet as a rational function can be made linear in its arguments, with some modifications. Specifically, look at how FunctionApproximations`RationalInterpolation works. It specifically uses LinearSolve. –  rcollyer Jun 8 '12 at 19:19
    
I gave RationalInterpolation a try, in general it seems to produce worse fits with more terms than my algorithm. I think this is because my algorithm eliminates terms and allows the power series to start in either the denominator or numerator from a power greater than 0. For instance, the data I fitted above has only 3 terms in the denominator/numerator, whereas I had to use 4 terms in the RationalInterpolation approach. –  Guillochon Jun 8 '12 at 20:47
1  
Is it possible to upload some sample data online somewhere, so we can experiment with fitting equations to it? –  Matt W-D Jul 27 '12 at 17:56
    
@MattW-D The sample data can be anything, the particular dataset I was fitting to is not more problematic than any other dataset. Probably the easiest way to construct a test dataset is to create a rational function with random coefficients, and then add some Gaussian noise: Total[(2.*RandomReal[]-1.)x^i, {i,1,5}]/Total[2.*RandomReal[]-1.)x^i,{i,1,5}] + 0.1*RandomVariate[NormalDistribution[]]. The resulting fit need not have the same power series with the same coefficients, just as long as it closely matches the dataset without having any sharp discontinuities. –  Guillochon Jul 30 '12 at 19:14

1 Answer 1

Try using a complex independent variable and fitting your function to the magnitude of the rational polynomial. This avoids the sharp discontinuities between data points due to poles in the denominator. With integer exponents the coefficients will still be real; I don't know what happens with fractional exponents.

myData = Function[x, {x, (-0.31 + 1.2*Sqrt[x] - 1.2*x)/(1. + 1.3*Sqrt[x] - 
    4.2*x)}] /@ Flatten[{Range[0.5, 1.5, 0.03], 2., 2.5}]; 
ratpoly[n_, m_] := Sum[a[k]*$s^(n - k), {k, 0, n}]/Sum[b[l]*$s^(m - l), 
  {l, 0, m}] /. b[m] -> 1; 
fitfnc = ratpoly[5, 5] /. $s -> I*w; 
vars = Cases[fitfnc, (a | b)[_], Infinity]; 
fit = fitfnc /. FindFit[myData, fitfnc, vars, w, NormFunction -> 
  (Abs[Norm[#1]] & ), Method -> "NMinimize"]; 
Show[{ListPlot[myData, PlotRange -> All, PlotStyle -> 
 Directive[PointSize[0.015], Opacity[0.5], ColorData[1, 3]]], 
Plot[Abs[fit], {w, 0.5, 2.5}, PlotRange -> All]}]
share|improve this answer
1  
Welcome to our site! Could you possibly add an explanation of how this approach addresses the principal question of using as few terms as possible? –  whuber Feb 9 '13 at 15:33
2  
You could improve your answer by adding some Mathematica code that illustrates the approach you are recommending. –  m_goldberg Feb 9 '13 at 17:16
    
Finally got around to testing this, it works pretty well. There are two concerns: Firstly, it complains with this error: Less::nord: "Invalid comparison with 0.480873 +0.` I", and secondly it claims to hit the MaxIterations in NMinimize. Increasing the MaxIterations actually returns a poorer solution! It would be great if this were a bit more robust, but it's definitely better than what I proposed (and in far fewer lines of code to boot!) –  Guillochon Apr 25 '13 at 22:05

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.