How to guess the correct fitting function to some data?

Question

Imagine we are plotting some date points $\left(x_i,f(x_i)\right)$ that we obtained experimentally, and that we want to know what $f(x)$ is. The way to do this is to use some software and try to fit the data to some guessed function. For example, if the behavior of the data points looks like exponential decay we then choose an exponential decaying function ..etc.

My question is: sometimes the data points are perfectly fitted to the exponential decaying trial function only on a certain region, but then the rest of the points show deviation away from the trial function.

1. How to proceed in this case? is there a catalog (something like mathematical tables) for functions and their plots that one can use as a guide?
2. Is there a systematic way to get the best fit instead of that trial and error method?

Answer 1

For example, if the behavior of the data points looks like exponential decay we then choose an exponential decaying function

No, if your theory predicts an exponential decay, then you use an exponential decay function. Or if your theory predicts a linear relationship, you use a linear fit. You really shouldn't have to guess which kind of function to use, because the "proper" way to analyze data is to test its consistency with some particular model, and the model tells you what kind of curve to expect.

This is something I feel isn't emphasized enough (or at all) in lab classes and such: if you don't have a model, the value of your data analysis is significantly diminished. In other words, just noticing that your data fit e.g. an exponential curve doesn't mean much by itself. It might hint at what sort of theory you should be looking at, but that's just taking shots in the dark (so to speak), not a proper analysis.

Answer 2

The possibilities are endless of course. You have to weigh the goodness of the fit against the prediction power of the hypothesis. You can fit any set of data "perfectly" by using as the "fit function" the data itself. This is of course complete nonsense since you would not gain any insight. The more degrees of freedom you admit, the better the fit can be (since, they're more parameters to adjust). With a common $\chi^2$ test, one oftentimes quotes the $\chi^2/\#dof.$, that is the sum of quadratic deviations divided by the number of degrees of freedom.

Ideally, when taking data, you already have an idea what kind of shape the data will have, since you would know the underlying physics process.

If this is not the case there are a number of ways to slice and dice the data to get more insights. In you example of exponential decay, you could try to have a sum of exponentials with different decay constants and different normalizations $\sum_i a_ie^{-x/a_i}$ and see if with a minimal extension your data fits better. (But again, you're losing prediction power).

Another possibility to dig up trends in data is Fourier analysis. This is useful when you suspect that the data has multiple overlapping periodical processes (modes). If this is the case you will see those nicely in a fourier transform of the data.

Still another common possibility is the convolution of two functions. This is for example the case when you're trying to measure a quantity but have a intrinsic uncertainty stemming from the apparatus. In that case the measurement would be a convolution of the true measurement values and the resolution function of the apparatus (often a good approximation is gaussian).

Having a table of functions would not really be useful, I think. In the end you would want to compare your data to some kind of theory, and theories are the better the more they explain (that is "fit data well") with the fewest possible number of assumptions (i.e. "parameters)

EDIT:

A neat example (including code in python) for Fourier Analysis is given here: http://linuxgazette.net/115/andreasen.html (Example 2). There, sunspot data is analyzed and it is found that there is a dominant frequency in the occurence of sunspots.

I'm no expert for a reference, but any introduction to signal processing or data analysis should do. This book by seems to have good rating on Amazon and seems to touch the main points: http://www.amazon.com/Introduction-Signal-Processing-Sophocles-Orfanidis/dp/0132091720/ref=ntt_at_ep_dpt_1

Answer 3

Physics, physics, physics.

You have to understand what is going on (this is what georg is hinting about in the comments): you shouldn't be "just guessing". If you're lucky there will be enough structure in the data to give you hint.

It is hard to be very specific without know the particulars of your situation.

The kinds of things to think about include

• Might there be a random process that fakes your signal? Probably you need a underling constant (i.e. fit to C + f(x)) if fitting in time. Non-random background may have a different structure. What other backgrounds exist? Can you measure them; or failing that can you model them?
• Is your function subtly wrong? Everyone reaches for a Gaussian first when fitting a peak, but if the process is resonant perhaps a Lorentzian would be better.
• Are there multiple channels through which this process can proceed. Maybe one (or more!) that you have neglected should be included in the fit.
• Are you sure your instrument is properly calibrated? Non-lineariteis especially can introduce artifacts in fitting.
• luksen mentions the convolution of your resolution with the data. This will have the biggest effect if the data changes rapidly over the scale of your instruments resolution. Fitting to convolved functions is harder, but can be done.

Note that when you get to the end you should have "the fit works if I add these fudge factors" but "the fit is clear after we account for all the back-grounds that we have discussed". That is you aren't done until you know why all the fudge factor are needed.

Finally, there is a bit of craft to this and you learn it by doing. Worse, most of the experience you will get is somewhat specific to the particular sub-discipline or measurement technique. That's life for you.

Answer 4

The theoretical decay model should be derived before you fit, and then the probability of the fit is given by the Baysian modification of your prior distribution for the parameters of the model. In this, I agree with the other answers.

But there are decay models which are sufficiently general that they can be used to fit large classes of experimental data, without fitting everything. If your data falls of exponentially for a while, then crosses over to falling like a power law with a slowly decreasing power, finally like a reciprocal logarithm, then like log-log, and then like log-log-log, you aren't going to get a good fit from any of these blind methods. But this is extraordinarily rare when studying physical systems--- such a complex decay only occurs in contrived mathematical situations, or when a system is actively moved from one phase to another according to a complicated plan.

Fitting by complete sets of functions

The generic way in which you fit arbitrary data that you feel should be approximated by a smooth curve is to run a best-fit polynomial. The polynomials are dense in the continuous functions, so you can always approximate anything, but you must use a polynomial of lowest order which fits the structure you believe is there. Best-fit lines are most common.

The standard measure of quality of fit is the sum of the squares of the deviation from the fit. Minimizing this gives the least-squares line, and it is also easy to find least squares polynomials of arbitrary order, so long as the order is less than the number of data points.

But this type of fit is wrong for decaying data. In this case you have options:

• If the measurements are very accurate when the quantity is small, you can take the logarithm, and fit the log-linear plot with a polynomial. This will give you the exponential decay as a first approximation, with higher order terms correcting this. The problem with this approach is that you are going to infinite time, and a polynomial always blows up at infinite time very quickly, so that its exponential probably disappears asymptotically faster than a physical function.
• You can try to use a Laplace transform. The Laplace transform of a pure decaying exponential is a delta function, and if you are using a Laplace transform for numerical data, you are just superposing exponentials to get a fit. You can start with a two-exponential fit of the form $f(t) = A e^{-at} + B e^{-bt}$, where you fit the four-constants by least squares (you can do this by steepest descent). This is also a complete expansion, because it is related to the Fourier transform by analytic continuation.
• You can fit to a falling power, which you should do anyway, because this is a common non-exponential decay. For example, if you have a pendulum decaying by air friction, the profile falls off as $1\over A+ Bt$, and this is always badly approximated by decaying exponentials. This is not a complete expansion, but if you find that the result is asymptotically a power, you can write $f(t) = g(t)\over A+Bt$, and fit g by Pade approximations (polynomials over polynomials).

There is no general absolute method, because each asymptotic limit is different, but once you know even just a little bit, you can extract the leading behavior and fit the rest using a complete expansion.

Answer 5

Finding good formulas to fit data is a simple task:
From Cornell Creative Machines Lab download 'Eureqa' (version II, codenamed 'formulize' for LINUX, W$, Mac) and use it. It will give you all the formulas and graphs you need.
If you need more processing power you can pay for cloud computing.

An example with SunSpots (monthly smoothed), $dS=S_{-1}-S_{-2}$:

Sunspots  delta_S   time  
S     dS        t  

58  4.6 1  
62.6    4.6 2  
70  4.6 3  
55.7    7.4 4  
85  -14.3   5  
83.5    29.3    6  
81.6    -1.5    7  
...  
...  
50.6    6.9 3153  
78  6.7 3154  
88  27.4    3155  
96.7    10  3156  
73  8.7 3157  

set TARGET function to be like this

S = f(dS, dS^2, dS^3, t)

it gives the formula ("Correlation Coefficient 0.99097013")

S = 0.4687512939*sma(dS, 34) + 0.321053344*sma((dS^2), 131) + 0.2694580519*sma((dS^2), 406) + -0.04853427812/sma((dS^2), 153) + sma(dS, 73) + sma(dS, 110) - 0.4618182571
where sma(variable,length) is the simple moving average
We can identify several important periods : 34, 73, 110, 131, 153, 406 months and both $dS$ and $(dS)^2$

-- set TARGET function to be like: compute the Derivative of S in function of S

D(S) = f(S)  

it gives the formula ("Correlation Coefficient" 0.96109605) (D(S)) = 0.02477032993 + 2.767538634*S + 0.8820551019*sma(S, 3) - 3.650621347*sma(S, 2) we see that we only need the last three values of monthly sunspot count to achive a good formula

With a Correlation Coefficient 0.99097013 (so close to 1 ;) I expect the formula is predictive to a certain degree.

note: to obtain the sma(variable,length) in excell I used, as an example:
column C = variable (dS)
column F, first row F1 = 76 (length)
next rows in column F will have sma(dS,76)

F2 =AVERAGE(OFFSET($C2;-CELL("contents";F$1)+1;0;CELL("contents";F$1);1))
and propagate to all rows ...