2

I was trying to fit this dataset:

#Mydataset damped sine wave data
#X ---- Y
45.80   320.0
91.60   -254.0
137.4   198.0
183.2   -156.0
229.0   126.0
274.8   -100.0
320.6   80.0
366.4   -64.0 
412.2   52.0
458.0   -40.0
503.8   34.0
549.6   -26.0
595.4   22.0
641.2   -18.0

which, as you can see by the plot below, has the classical trend of a damped sine wave: enter image description here

So i first set the macro for the fit

f(x) = exp(-a*x)*sin(b*x)

then i made the proper fit

fit f(x) 'data.txt' via a,b
iter      chisq       delta/lim  lambda   a             b            
   0 2.7377200000e+05   0.00e+00  1.10e-19    1.000000e+00   1.000000e+00

Current data point
=========================
#               = 1 out of 14
x               = -5.12818e+20   
z               = 320            

Current set of parameters
=========================
a               = -5.12818e+20   
b               = -1.44204e+20   

         Function evaluation yields NaN ("not a number")

getting a NaN as result. So I looked around on STackOverflow and I remembered I've already have had in the past problems by fitting exponential due to their fast growth/decay which requires you to set initial parameters in order not to get this error (as I've asked here). So I tried by setting as starting parameters a and b as the ones expected, a = 9000, b=146000, but the result was more frustrating than the one before:

fit f(x) 'data.txt' via a,b
iter      chisq       delta/lim  lambda   a             b            
   0 2.7377200000e+05   0.00e+00  0.00e+00    9.000000e+03   1.460000e+05
         Singular matrix in Givens()

I've thought: "these are too much large numbers, let's try with smaller ones". So i entered the values for a and b and started fitting again

a = 0.01
b = 2
fit f(x) 'data.txt' via a,b
iter      chisq       delta/lim  lambda   a             b            
   0 2.7429059500e+05   0.00e+00  1.71e+01    1.000000e-02   2.000000e+00
   1 2.7346318324e+05  -3.03e+02  1.71e+00    1.813940e-02  -9.254913e-02
   * 1.0680927157e+137  1.00e+05  1.71e+01   -2.493611e-01   5.321099e+00
   2 2.7344431789e+05  -6.90e+00  1.71e+00    1.542835e-02   4.310193e+00
   * 6.1148639318e+81   1.00e+05  1.71e+01   -1.481123e-01  -1.024914e+01
   3 2.7337226343e+05  -2.64e+01  1.71e+00    1.349852e-02  -9.008087e+00
   * 6.4751980241e+136  1.00e+05  1.71e+01   -2.458835e-01  -4.089511e+00
   4 2.7334273482e+05  -1.08e+01  1.71e+00    1.075319e-02  -4.346296e+00
   * 1.8228530731e+121  1.00e+05  1.71e+01   -2.180542e-01  -1.407646e+00
   * 2.7379223634e+05   1.64e+02  1.71e+02    8.277720e-03  -1.440256e+00
   * 2.7379193486e+05   1.64e+02  1.71e+03    1.072342e-02  -3.706519e+00
   5 2.7326800742e+05  -2.73e+01  1.71e+02    1.075288e-02  -4.338196e+00
   * 2.7344116255e+05   6.33e+01  1.71e+03    1.069793e-02  -3.915375e+00
   * 2.7327905718e+05   4.04e+00  1.71e+04    1.075232e-02  -4.332930e+00
   6 2.7326776014e+05  -9.05e-02  1.71e+03    1.075288e-02  -4.338144e+00
iter      chisq       delta/lim  lambda   a             b            

After 6 iterations the fit converged.
final sum of squares of residuals : 273268
rel. change during last iteration : -9.0493e-07

degrees of freedom    (FIT_NDF)                        : 12
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 150.905
variance of residuals (reduced chisquare) = WSSR/ndf   : 22772.3

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = 0.0107529        +/- 3.114        (2.896e+04%)
b               = -4.33814         +/- 3.678        (84.78%)

correlation matrix of the fit parameters:
                a      b      
a               1.000 
b               0.274  1.000 

I saw it produced some result, so I thought it was all ok, but my happiness lasted seconds, just until I plotted the output:enter image description here Wow. A really good one.

And I'm still here wondering what's wrong and how to get a proper fit of a damped sine wave dataset with gnuplot. Hope someone knows the answer :)

opisthofulax
  • 517
  • 8
  • 25

1 Answers1

5

The function you are fitting the data to is not a good match for the data. The envelope of the data is a decaying function, so you want a positive damping parameter a. But then your fitting function cannot be bigger than 1 for positive x, unlike your data. Also, by using a sine function in your fit you assume something about the phase behavior -- the fitted function will always be zero at x=0. However, your data looks like it should have a large, negative amplitude.

So let's choose a better fitting function, and give gnuplot a hand by choosing some reasonable initial guesses for the parameters:

f(x)=c*exp(-a*x)*cos(b*x)
a=1./500
b=2*pi/100.
c=-400.
fit f(x) 'data.txt' via a,b,c
plot f(x), "data.txt" w p

gives

enter image description here

user8153
  • 4,049
  • 1
  • 9
  • 18
  • Smart solution! Thanks a lot. I'm only wondering now: how did you choose the guessed values to set the initial parameters `a`, `b` and `c`? What was your method? – opisthofulax Jun 12 '17 at 17:41
  • 2
    1/a is the "decay length" of your exponential, so looking at your data a value of a few hundred seems about right. The period of oscillation of the cosine curve seems to be about 100, and b should be 2*pi divided by that number. At x=0 one would expect a large negative amplitude, and -400 looks about right. – user8153 Jun 12 '17 at 17:52