Fitting data using broken power-law in gnuplot

Question

I want to fit my data with broken power law function. I am using data from two different files and this is the code that I am using to fit my data in gnuplot tool

set term wxt
    p 'Data1.dat' u 1:($8*100):2:($9*100) w xyerr lt 7 t '', \
      'Data2.dat' u 1:($8*100):2:($9*100) w xyerr lt 8 t ''
    set log x 2
    f(x) = A*(x**p)*(x**(p-q))
    A = 1.0
    p = 1.5
    q = 0.1
    rep f(x) lt 8 lw 2 t ''
    fit f(x) '< cat Data1.dat Data2.dat' u 1:($8*100):($9*100) yerr via A,p,q
    replot`

But my fitting looks like this:

enter image description here

Is it the way I am using function wrong or something else?

With "broken power-law", do you mean [this thing](https://en.wikipedia.org/wiki/Power_law#Broken_power_law)? If yes, do you want to fit the "break point" as well? — maij, Jan 21 '19 at 19:44
Well, your starting parameters seem to be way off. That function doesn't look remotely like your data with them. Gnuplot can't do miracles, you have to supply reasonable values, except if you fit a pure polynome. — Karl, Jan 23 '19 at 17:55
@maij I was trying to find appropriate function that I can fit for my data and I found broken power to be the best one. As I am using data from two different files that makes entire plot with increasing X- axis value I am not concerned about "break point" as such — Arya Stark, Jan 24 '19 at 07:11
@Karl Function does fit when used in other plotting devices with pre-defined broken power-law. I am finding it difficult how to use the same function in gnuplot — Arya Stark, Jan 24 '19 at 07:15
I don't know if the function you defined is the one you wanted (you still haven't answered maij's question above !?! ), but if it is, please look at the plot before you do the fitting, and you see that it is _not at all_ describing you data. The parameter values A=1, p = 1.5, q=0.1 are definitively totally wrong. No fitting algorithm will succeed from such a starting point. — Karl, Jan 24 '19 at 21:09
@Arya Stark, did my answer solve your problem? Any kind of response would be reputable and appreciated. — theozh, Nov 20 '20 at 10:46
@theozh Yes! thank you so much. Extremely sorry for not responding for so long. — Arya Stark, May 16 '21 at 06:44
@AryaStark no problem. Thank you on StackOverflow = accept/upvote. — theozh, May 16 '21 at 07:19

theozh · Answer 1 · 2019-01-24T21:41:03.467

1

As I understand you basically want to fit different functions in different ranges. So, just use two functions in different ranges. Maybe something like this...

Edit: added a continuous function h(x). Data approximately from OP's graph.

# SO_data.dat
0.551   2.213
0.928   3.531
1.199   4.796
1.461   5.901
1.963   6.393
2.770   6.260
3.760   5.794
4.445   5.515
4.905   5.528
5.914   5.581
7.566   5.062
4.358   4.996
5.052   4.929
6.032   4.729
6.924   4.609
7.948   4.370
8.945   4.117
10.167  4.024
11.902  3.930
14.928  3.824
18.724  3.704
23.484  3.438
29.166  3.584
42.405  2.945

And the code:

### fitting two regions
reset session
set colorsequence classic

set logscale x 2
FILE = "SO_data.dat"
set xrange[0.25:64]
set yrange[2:9]

# some start values
A = C = 4
p = r = 0.8
B = D = 8
q = s = -0.3
d = 2
a = 3
f(x) = A*x**p 
g(x) = B*x**q
h(x) = C*x**r/(exp((x-d)*a)+1) + D*x**s/(exp((-x+d)*a)+1)

fit [:2] f(x) FILE u 1:2 via A,p
fit [2:] g(x) FILE u 1:2 via B,q
fit h(x) FILE u 1:2 via C,D,r,s,a,d

c = (B/A)**(1/(p-q))   # crossing point
print sprintf("A: %.3g, p: %.3g, B: %.3g, q: %.3g, c: %.3g",A,p,B,q,c)
print sprintf("C: %.3g, r: %.3g, D: %.3g, s: %.3g, a: %.3g, d: %.3g",C,r,D,s,a,d)
plot FILE u 1:2 w p ps 2,\
    f(x) noautoscale, g(x) noautoscale, h(x) noautoscale
### end of code

Output:

A: 3.96, p: 0.795, B: 8.1, q: -0.274, c: 1.95
C: 1.43, r: 0.046, D: 8.08, s: -0.272, a: 3.63, d: 1.15

edited Jan 24 '19 at 21:41

answered Jan 24 '19 at 09:43

theozh

22,244
5
28
72

You can make a piecewise definition in gnuplot `f(x) = x < t ? A(x-t)**2 : B*(x-t)**3`, or `r=2; f(x) = x – Karl Jan 24 '19 at 21:18
...haven't tested whether gnuplot also fits a piece-wise function. I added a single-piece and smooth function which looks nicer ;-) – theozh Jan 24 '19 at 21:43
Maybe, but it's probably unphysical. ;-) The Marquard-Levenberg algorithm (and gnuplot's implementation) doesn't care how your function is defined, as long as it returns function values. – Karl Jan 24 '19 at 21:50
Admittedly, the function h(x) you added looks like it actually *is* a good model. OP probably want's to use that one. – Karl Jan 24 '19 at 21:54
...well, the part `1/(exp((x-d)*a)+1)` is the same structure as Fermi-Dirac statistics, so I would say _very_ physical ;-) – theozh Jan 24 '19 at 21:55
You could still add a bit of Einstein or Heisenberg. :-)) – Karl Jan 24 '19 at 21:58

Fitting data using broken power-law in gnuplot

1 Answers1