23
$\begingroup$

I have a series of data points $(x_i,y_i)$ which I expect to (approximately) follow a function $y(x)$ that asymptotes to a line at large $x$. Essentially, $f(x) \equiv y(x) - (ax + b)$ approaches zero as $x \to \infty$, and the same can probably be said of all the derivatives $f'(x)$, $f''(x)$, etc. But I don't know what the functional form for $f(x)$ is, if it even has one that can be described in terms of elementary functions.

My goal is to get the best possible estimate of the asymptotic slope $a$. The obvious crude method is to pick out the last few data points and do a linear regression, but of course this will be inaccurate if $f(x)$ does not become "flat enough" within the range of $x$ for which I have data. The obvious less-crude method is to assume that $f(x) \approx \exp(-x)$ (or some other particular functional form) and fit to that using all the data, but the simple functions I've tried like $\exp(-x)$ or $\dfrac1{x}$ don't quite match the data at lower $x$ where $f(x)$ is large. Is there a known algorithm for determining the asymptotic slope that would do better, or that could provide a value for the slope along with a confidence interval, given my lack of knowledge of exactly how the data approach the asymptote?


This sort of task tends to come up frequently in my work with various data sets, so I'm mostly interested in general solutions, but by request I'm linking to the particular data set that prompted this question. As described in comments, the Wynn $\epsilon$ algorithm gives a value that, as far as I can tell, is somewhat off. Here is a plot:

Asymptotically linear data

(It does look like there's a slight downward curve at high x values, but the theoretical model for this data predicts that it should be asymptotically linear.)

$\endgroup$
  • $\begingroup$ This may well be too elementary - or too vague - for this site, but I figured private beta is the time to try such things out. $\endgroup$ – David Z Dec 6 '11 at 8:53
  • $\begingroup$ No, I think this is a great question. Not everything must be advanced and fancy. Good solutions to simple problems are important. $\endgroup$ – Colin K Dec 6 '11 at 16:42
  • $\begingroup$ @Dan: was replacing $\exp$ really justified? $\endgroup$ – J. M. Dec 6 '11 at 19:14
  • $\begingroup$ Having exps tends to make things harder for me to read, but I'll admit that it was small enough that I shouldn't have done it. $\endgroup$ – Dan Dec 6 '11 at 19:20
  • $\begingroup$ I really don't care either way, I just figured I might as well approve the edits because, well, why not. You get a couple reputation off it, whatever that's worth. $\endgroup$ – David Z Dec 6 '11 at 19:33
13
$\begingroup$

It's a rather rough algorithm, but I'd use the following procedure for a crude estimate: if, as you say, the purported $f(x)$ that represents your $(x_i,y_i)$ is already almost linear as $x$ increases, what I'd do is to take differences $\dfrac{y_{i+1}-y_i}{x_{i+1}-x_i}$, and then use an extrapolation algorithm like the Shanks transformation to estimate the limit of the differences. The result is hopefully a good estimate of this asymptotic slope.


What follows is a Mathematica demonstration. The Wynn $\epsilon$ algorithm is a convenient implementation of the Shanks transformation, and it is built in as the (hidden) function SequenceLimit[]. We try out the procedure on the function

$$\frac4{x^2+3}+2 x+e^{-4 x}+3$$

xdata = RandomReal[{20, 40}, 25];
ydata = Table[(3 + 13*E^(4*x) + 6*E^(4*x)*x + x^2 + 3*E^(4*x)*x^2 + 
      2*E^(4*x)*x^3)/(E^(4*x)*(3 + x^2)), {x, xdata}];

SequenceLimit[Differences[ydata]/Differences[xdata],
              Method -> {"WynnEpsilon", Degree -> 2}]
1.999998

I might as well show off how simple the algorithm is:

wynnEpsilon[seq_?VectorQ] := 
 Module[{n = Length[seq], ep, res, v, w}, res = {};
  Do[ep[k] = seq[[k]];
   w = 0;
   Do[v = w; w = ep[j];
    ep[j] = 
     v + (If[Abs[ep[j + 1] - w] > 10^-(Precision[w]), ep[j + 1] - w, 
         10^-(Precision[w])])^-1;, {j, k - 1, 1, -1}];
   res = {res, ep[If[OddQ[k], 1, 2]]};, {k, n}];
  Flatten[res]]

Last[wynnEpsilon[Differences[ydata]/Differences[xdata]]]
1.99966

This implementation is adapted from Weniger's paper.

$\endgroup$
  • $\begingroup$ Just curious, but why did you the original form of the function, instead of combining all of the terms? $\endgroup$ – rcollyer Dec 6 '11 at 13:37
  • $\begingroup$ It was only for demonstration purposes. :) I included the $\LaTeX$-ed expression so you guys know what the expected answer was supposed to be. What I wanted to demonstrate was that you can do this sort of analysis on some complicated-looking function... $\endgroup$ – J. M. Dec 6 '11 at 13:40
  • $\begingroup$ How close to flat do the points have to be for the algorithm to be effective? $\endgroup$ – rcollyer Dec 6 '11 at 14:01
  • 2
    $\begingroup$ Okay, last question (I swear), can you generate an error bound on the estimate? $\endgroup$ – rcollyer Dec 6 '11 at 18:20
  • 1
    $\begingroup$ That's a bit trickier. I've seen some suggested methods in a few papers, but I confess to not having done experiments with them. (Maybe I should, one of these days.) The book by Brezinski and Redivo-Zaglia has a few pointers you might want to look into. $\endgroup$ – J. M. Dec 6 '11 at 18:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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