gnuplotdistributiondata-fitting
How to fit a bimodal distribution function via Gnuplot?
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 ?
Solution
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:
