Arya Stark
Arya Stark

Reputation: 11

Fitting data using broken power-law in gnuplot

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?

Upvotes: 0

Views: 1226

Answers (1)

theozh
theozh

Reputation: 26088

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

enter image description here

Upvotes: 1

Related Questions