gnuplotdata-fitting

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?


Solution

  • 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.

    Script:

    ### 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 script
    

    Result:

    enter image description here

    Addition: (irregular spacing in x)

    In case your x-data is not spaced by 1 and not even equidistant you would have to adjust the procedure to calculate the area under a curve.

    The area under the curve between two datapoints (x0,y0) and (x1,y1 is simply (x1-x0)*(y0+y1)*0.5. So, you can calculate the area under the curve during plotting the original data. In the example below, the plotting line with impulses is just there to better illustrate the different spacings of the datapoints.

    Script:

    ### calculate area under curve for irregular x-spacings
    reset session
    
    # create some random test data
    set table $Data
        N = 15
        set samples N
        plot x0=0 '+' u (x0=x0+rand(0)*10):(100*exp(-((x0-N*2)/N)**2)) w table
    unset table
    
    set style fill solid 0.3 border
    set key noautotitle
    
    plot auc=0 $Data u ($0==0?(x1=$1,y1=$2):0, y0=y1, y1=$2, x0=x1, x1=$1): \
         (auc=auc+(x1-x0)*(y0+y1)*0.5,y1) w lp pt 7 ps 0.5 lw 2, \
         '' u 1:2 w impulses lc "red", \
         '' u 1:2 w filledcurves y=0 lc rgb 0xccff0000 ti sprintf("Area: %g",auc)
    ### end of script
    

    Result:

    enter image description here