CODEWITHSUNDEEP
gnuplotdata-fitting
mas
mas

Reputation: 35

Calculate area under curve in gnuplot

I have data (data can be downloaded here: gauss_data) and need to find the area of a particular peak. From my data set, the peak seems to have some contribution from another peak. I made the fit on my data with 3 Gaussians using this code:

# Gaussian fit                                              

reset 
set terminal wxt enhanced

# Set fitting function      

f(x)  = g1(x)+g2(x)+g3(x)
g1(x) = p1*exp(-(x-m1)**2/(2*s**2))
g2(x) = p2*exp(-(x-m2)**2/(2*s2**2))
g3(x) = p3*exp(-(x-m3)**2/(2*s3**2))


# Estimation of each parameter      

p1 = 100000
p2 = 2840
p3 = 28000
m1 = 70
m2 = 150
m3 = 350
s = 25
s2 = 100
s3 = 90

# Fitting & Plotting data                           

fit [0:480] f(x) 'spectrum_spl9.txt' via p1, m1, s, p2, m2, s2, p3, m3, s3
plot [0:550] 'spectrum_spl9.txt' lc rgb 'blue', f(x) ls 1, g1(x) lc rgb 'black', g2(x) lc rgb 'green' , g3(x) lc rgb 'orange'

and the result is shown in fig below

enter image description here

I need to calculate the area under the peak i.e. area f(x) - area g3(x). Is there any way to find/calculate the area of each function in Gnuplot?

Upvotes: 1

Views: 2778

Answers (2)

theozh
theozh

Reputation: 25714

Your data is equidistant in x-units with a step width of 1. So, you can simply sum up the intensity values multiplied by the width (which is 1). If you have irregular data then this would be a bit more complicated.

Code:

### determination of area below curve
reset session
FILE = "SO/spectrum_spl9.txt"

# fitting function      
f(x)  = g1(x)+g2(x)+g3(x)
g1(x) = p1*exp(-(x-m1)**2/(2*s1**2))
g2(x) = p2*exp(-(x-m2)**2/(2*s2**2))
g3(x) = p3*exp(-(x-m3)**2/(2*s3**2))

# Estimation of each parameter      
p1 = 100000
p2 = 2840
p3 = 28000
m1 = 70
m2 = 150
m3 = 350
s1 = 25
s2 = 100
s3 = 90

set fit quiet nolog
fit [0:480] f(x) FILE via p1, m1, s1, p2, m2, s2, p3, m3, s3

set table $Difference
    plot myIntegral=0 FILE u 1:(myIntegral=myIntegral+f($1)-g3($1),f($1)-g3($1)) w table
unset table

set samples 500    # set samples to plot the functions
plot [0:550] FILE u 1:2 w p lc 'blue' ti FILE noenhanced, \
      f(x) ls 1, \
      g1(x) lc rgb 'black', \
      g2(x) lc rgb 'green', \
      g3(x) lc rgb 'orange', \
      $Difference u 1:2 w filledcurves lc rgb 0xddff0000 ti sprintf("Area: %.3g",myIntegral)
### end of code

Result:

enter image description here

Upvotes: 1

maij
maij

Reputation: 4218

Can you use the analytic integral under a Gaussian function?

y(x) = 1/(s*sqrt(2*pi)) * exp(-(x-m1)**2/(2*s**2))
integral(y) [-inf:inf] = 1

This would mean that:

I1 = integral(g1) = p1 * s1 * sqrt(2.0*pi)
I2 = integral(g2) = p2 * s2 * sqrt(2.0*pi)

area f(x) - area g3(x) = I1 + I2

Please double check the math :)

Upvotes: 1

Related Questions