How to fit a bimodal distribution function via Gnuplot?

Question

I have the following data which I like to fit with a Bimodal distribution

     <0       0.000000e+00         0          0.000000e+00    0.000000e+00         0          0.000000e+00   
     <=0.3       3.000000e-01         1          8.333333e-02    8.333333e+00         1          8.333333e+00   
     <=0.6       6.000000e-01         0          0.000000e+00    0.000000e+00         1          8.333333e+00   
     <=0.9       9.000000e-01         2          1.666667e-01    1.666667e+01         3          2.500000e+01   
     <=1.2       1.200000e+00         6          5.000000e-01    5.000000e+01         9          7.500000e+01   
     <=1.5       1.500000e+00         3          2.500000e-01    2.500000e+01        12          1.000000e+02   
     <=1.8       1.800000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
     <=2.1       2.100000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
     <=2.4       2.400000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
     <=2.7       2.700000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
       <=3       3.000000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
        >3       3.000000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   

pl 'data.txt' u 2:($5/100) w lp
gauss(x)=a/(sqrt(2*pi)*sigma)*exp(-(x-mean)**2/(2*sigma**2))
fit gauss(x) 'data.txt' u 2:($5/100) via a,sigma,mean

This is the result and with the fit I just get 1 peak

Can anyone suggest how to proceed? Is it possible to fit this data via Gnuplot ?

Answer 1

Simply define your fitting function consisting of two Gaussian peaks. And help the fitting algorithm with some reasonable starting values.

Code:

### fitting two gauss peaks
reset session

$Data <<EOD
     <  0       0.000000e+00         0          0.000000e+00    0.000000e+00         0          0.000000e+00   
     <=0.3       3.000000e-01         1          8.333333e-02    8.333333e+00         1          8.333333e+00   
     <=0.6       6.000000e-01         0          0.000000e+00    0.000000e+00         1          8.333333e+00   
     <=0.9       9.000000e-01         2          1.666667e-01    1.666667e+01         3          2.500000e+01   
     <=1.2       1.200000e+00         6          5.000000e-01    5.000000e+01         9          7.500000e+01   
     <=1.5       1.500000e+00         3          2.500000e-01    2.500000e+01        12          1.000000e+02   
     <=1.8       1.800000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
     <=2.1       2.100000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
     <=2.4       2.400000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
     <=2.7       2.700000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
       <=3       3.000000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
        >3       3.000000e+00         0          0.000000e+00    0.000000e+00        12          1.000000e+02   
EOD

Gauss(x,a,mean,sigma) = a/(sqrt(2*pi)*sigma)*exp(-(x-mean)**2/(2*sigma**2))
g1(x) = Gauss(x,a1,mean1,sigma1)
g2(x) = Gauss(x,a2,mean2,sigma2)
f(x) = g1(x) + g2(x)

# set some initial values to help the fitting algorithm
mean1 = 0.3
mean2 = 1.2

set fit results nolog
fit f(x) $Data u 2:($5/100) via a1,mean1,sigma1,a2,mean2,sigma2

set samples 200

myResults = sprintf("a1 = %g\nmean1 = %g\nsigma1 = %g\na2 = %g\nmean2 = %g\nsigma2 = %g\n",a1,mean1,sigma1,a2,mean2,sigma2)
set label 1 at graph 0.05, 0.97 myResults

plot $Data u 2:($5/100)  w lp pt 7 dt 3 ti "Data", \
     f(x) w l lc "red"
### end of code

Result: