How to find the right initial values for non-linear fits of a Sinc²(x) function in GnuPlot?

Question

I have experimental data from optical measurements. In a particular regime, but I will skip the details, the data received must obey a cardinal sine squared law, which is indeed what I observe (approximately):

I want to perform a cardinal sine fit on my data in order to recover the parameters of interest (and their statistical error). I give the data here (there are more than 2K lines, so I'm not going to put it here unless I'm asked to). My code is the following:

sinc(x) = sin(x)/x
f(x) = a * sinc(b * x + c)**2 + d
fit f(x) 'data.txt' using 1:3 via a, b, c, d
plot 'data.txt' using 1:3, f(x) title 'Fit'

Without an initial value, GnuPlot of course has trouble fitting anything, it makes me a line. So I have to choose parameters, as it is the central peak that is prominent, I have based several initial values around it, but no matter what I choose, a will always be negative in the end, in the parameters deduced by GnuPlot, and I don't understand why. How do I modify this code and choose the right parameters to make it work?

Answer 1

I guess your fit does not work, because you have to first shift your x-value (x-c) and then scale it b*(x-c), i.e.

f(x) = a*sinc(b*(x-c))**2 + d

With this you can easily find reasonable starting values:

• a maximum intensity, ~3500
• b the basic function (sin(x)/x)**2 has a FWHM of ~2.8, and your data has a FWHM of about 800. Hence with b*800=2.8 gives approx. b=0.035.
• c x-shift, i.e. max position, ~1100. Since your x-values are integer values from 1 to 2000, with 1100 you will get an argument which is exactly 0. However, then sin(x)/x is NaN and the fitting algorithm doesn't like that: Error: Undefined value during function evaluation. So, simply choose c=1100.5 as starting value.
• d y-shift, first guess would be 0, but the fitting algorithm does not like 0 as starting value, so then try 1.0

You have to decide yourself if this is a good fit. The secondary maxima are not reproduced very well in their position. This might be because the experimental global maximum seems to be lower than expexted value from the function.

Script:

### fit sinc(x)**2 function
reset session

FILE = "SO76186402.dat"

set datafile separator ";"

set samples 1000
sinc(x) = sin(x)/x
f(x) = a*sinc(b*(x-c))**2 + d

# starting values estimates
a = 3500
b = 0.035
c = 1100.5
d = 1.0

set fit nolog
fit f(x) FILE u 1:3 via a,b,c,d

plot FILE u 1:3 w l lw 1.5 lc "blue" ti "data", \
     f(x) w l lc "red" ti "fit"
### end of script

Result:

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = 3477.59          +/- 4.765        (0.137%)
b               = 0.00622796       +/- 9.026e-06    (0.1449%)
c               = 1060.32          +/- 0.2737       (0.02582%)
d               = 17.2389          +/- 2.054        (11.92%)

Addition:

Just in case you need a better fit for the secondary maxima, you could fit the logarithm of the function and of your data.

Script:

### fit sinc(x)**2 function
reset session

FILE = "SO76186402.dat"

set datafile separator ";"

set samples 1000
sinc(x) = sin(x)/x
f(x) = a*sinc(b*(x-c))**2 + d

# estimates
a = 3700
b = 0.0067
c = 1055.3
d = 1.0

set multiplot layout 4,1

    fit f(x) FILE u 1:3 via a,b,c,d
    plot FILE u 1:3 w l lw 1.5 lc "blue" ti "data", \
         f(x) w l lc "red" ti "fit"

    set logscale y
    replot
    
    fit log(f(x)) FILE u 1:(log($3)) via a,b,c,d
    unset logscale y
    replot

    set logscale y
    replot
unset multiplot
### end of script

Result:

Fit of f(x):

• 1st graph linear scale
• 2nd graph logarithmic scale

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = 3477.59          +/- 4.765        (0.137%)
b               = 0.00622794       +/- 9.025e-06    (0.1449%)
c               = 1060.32          +/- 0.2737       (0.02582%)
d               = 17.2365          +/- 2.054        (11.92%)

Fit of log(f(x)):

• 3rd graph linear scale and
• 4th graph logarithmic scale

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = 4070.87          +/- 36.04        (0.8853%)
b               = 0.00667818       +/- 5.032e-06    (0.07535%)
c               = 1051.48          +/- 0.4286       (0.04076%)
d               = 9.99952          +/- 0.2293       (2.293%)