Double gaussian fit using FuncFit

Share icon

Thu, 01/22/2015 - 01:17 pm

Hello everyone..

I am trying to fit a doublet using the FuncFit where I need to tell Igor the initial guesses as Location, Amplitude and Width. The displayhelptopic about FuncFit tells that we can inform the initial guesses in a wave using this: kwCWave=wGuessesWave , but it dosn't tell about which information on which position. Searching in the internet I have found an exemple and I am adopting that the positions in the wave for amplitude, width and location are correct!

However.. I have one error that I don't know how to fix it: "singular matrix or other numeric error", and in the command window this massage is printed: There may be no dependence on these parameters: W_coef[4],W_coef[5],W_coef[6]. Probably it is a silly mistake done by me.

By the way, is there any way to have the FWHM and Chi Square in the W_coefs result wave?

Thanks so much

Andre



Function doublefit()
    wave wSpectrum=$("FR_Y")
    wave wWavenumber=$("FR_X")
   
      variable StartP, EndP
             // Searching for Point Numbers on X Wave
             FindValue /V=690 /T=0.5 wWavenumber
             StartP=v_value    
             FindValue /V=720 /T=0.5 wWavenumber
             EndP=v_value    

             Make /O/N=7 root:W_coef
             Wave W_coef=root:W_coef
            W_coef[0]=0 //offset
            W_coef[1]=1 // peak height
            W_coef[2]=700.5 // centre
            W_coef[3]=2.042 //width
            W_coef[4]=1 // peak height
            W_coef[5]=705 // centre
            W_coef[6]=2.462 //width            
             
             Print W_coef

        // fit a double gaussian, hold the centre values (W_coef[2] and W_coef[5])
        FuncFit /H="0000000"/NTHR=0/N/W=0 DoubleGauss  kwCWave=W_coef, wSpectrum[StartP,EndP] /X=wWavenumber
        // fit results stored in W_coef, which is effectively returned by reference
End
Hello

My first guess would be a "typo" in the fit function. Could you please post it?
Otherwise it could be that your two peaks "merge": One peak takes it all and the other one vanishes. In this case you might have to use constraints (->manual).

BTW: FuncFit officially does not take the "kwCWave" keyword (CurveFit does!):
FuncFit [ flags ] fitFuncName, cwaveName, waveName [ flag parameters ]

Is there any reason why you use an x-wave? It's labeled "wavenumber" which makes me think it could have constant spacing? In this case apply "setscale" to your data wave and do not use an x wave at all. You might also specify a fit range just by /R=(StartX,EndX) in this case ( note: "()" vs "[]" ).

The fwhm should be correlated to your width coefficients and your fit function (maybe the widths are the same?).
Chi square comes from the "StatsChiTest" applied on your data and the fit result, please see the manual for details.

Your error message tells you, that you have a linear dependency on your coefficients, i.e., the temporary matrix does not have full rank (but it need to have it).
y= a + bx +cx would be an example for this case; only (a+b) matters for y; you can give any value for b (or a), and a (or b) will adjust. But adjusting both will not work.

HJ

January 22, 2015 at 03:11 pm - Permalink

Hello HJDrescher, first of all thanks for your considerations!

HJDrescher wrote:

My first guess would be a "typo" in the fit function. Could you please post it?
Otherwise it could be that your two peaks "merge": One peak takes it all and the other one vanishes. In this case you might have to use constraints (->manual).


What do you mean by "typo"? I am pretty sure that those peaks are not convoluted! I can fit them in the manual mode using the Multipeak fitting 2.0 option, and those values I am suggesting in the wave are pretty close from the results I got from manual fitting!

HJDrescher wrote:

Is there any reason why you use an x-wave? It's labeled "wavenumber" which makes me think it could have constant spacing? In this case apply "setscale" to your data wave and do not use an x wave at all. You might also specify a fit range just by /R=(StartX,EndX) in this case ( note: "()" vs "[]" ).


Yes.. that's the information we use for Raman Spectroscopy, and we need the peak info (FWHM and Position) related to the X axis. About the scale, it suppose to be constant, but sometimes we have light changes on it, that's why I am keeping the X-axis.

HJDrescher wrote:

Chi square comes from the "StatsChiTest" applied on your data and the fit result, please see the manual for details.


I am using the Curvefit function to fit the other peaks on my graph, and at the end the function plots a few results in the command window (including the Chi Square). But they are not in the result's wave. My next goal is to do batch fitting, so that's why I asked if it is possible to include the Chi square in the wave!

Can I use the Curvefit function to fit double peaks or it is just for single peaks?

Thanks again!

Andre

January 23, 2015 at 01:56 am - Permalink

Hello Andre

a typical "typo" would be typing the function for the first peak and then copy the text for the second one and forgetting to change "w[2]" to "w[5]". The widths might depend on the use of "[]" or "()" in the fit function (x-waves...) - They have a different "unit": wave numbers or data points. This is why I would like to have a look at the function - just post it.

Maybe your FindValue construct does not work as intended? You should check whether it returns -1 and if StartP and EndP are different. (You might end up fitting ONE data point outside your spectrum).
For debugging a print StartP,EndP will do.

Regarding the x wave: If most of the data points are equidistant you might consider to use interpolate to get a "normal" data wave. It makes life much easier.

Finally the chi²: I cannot find the option for automatic calculation in the manual. But you can add the /D=FittedWave flag at the end of the FuncFit-line and you get one of the waves you need. The other one is your (ranged) source wave.
StatsChiTest /S wSpectrum[StartP,EndP], FittedWave

Should provide a wave with the results of the chi square test. Have a look at the StatsChiTest-Demo.

For sums of fit functions and user defined functions you need FuncFit. Please also have a look at "Fitting Sums of Fit Functions" mentioned in the FuncFit help.

HJ

January 23, 2015 at 03:55 am - Permalink

Hello again HJ

I still don't understand what I suppose to do with this "typo" thing... is it a function that may return any information about my waves or something like that?

I did the test printing the StartP and EndP values, and they are correct ! They are supposed to be, because in my other peaks (the single ones) the procedure works fine, but for them I am using the CurveFit.

I have uploaded a figure where shows the doublet I am trying to fit via procedure. Sometimes the auto-guess (for the manual mode) finds really good parameters like in the figure, but sometimes not and then I have to suggest values to have the perfect fitting. So, that's why I have to inform to the IGOR where my peaks are!

I included in the same figure a small part of my X-wave... the step between lines are around 0.71, having differences in the third decimal case. The uncertainty of my equipment is around 0.7 or 0.8, so do think it is a problem discarding the information in the third decimal case. But what is the advantage of doing that?

For now I will discard the Chi Square info.. first of all, I have to find a way to solve this fitting problem, and create a looping for batch fitting!

Thanks again for your time.. and as you can see.. I am a beginner with Igor!

Andre
fit.jpg

January 23, 2015 at 06:30 am - Permalink

Hello Andre,

you should have a function called "DoubleGauss". Please post it.

The fact that your StartP and EndP values are correct in this case is good news but your code should test it and print an error message is something is wrong (including StartP>EndP).

Your figure suggest that a fit might fail with a poor initial guess (especially if the width parameter is too large) and it might be necessary to use constraints (one peak position <702 the other one >703).
FuncFit always need initial values. In some cases they can be off by a lot, sometimes not; it depends on the function you are using for fitting (including constraints).

There is a chapter in the Igor manual about the waveform model (Chapter II-5):
Quote:
Some people who have uniformly spaced data still use the XY model because it is what they are accustomed to. This is a mistake. If your data is uniformly spaced, it will be well worth your while to learn and use the waveform model. It greatly simplifies graphing and analysis and makes it easier to write Igor procedures.


One your fitting routine works the chi square should be easy.
A hint for batch fitting: If your peaks change only a little between two spectra use the results from the previous fit as guesses for the next one.

Things will work out just fine
HJ

January 23, 2015 at 07:20 am - Permalink

Yes.. there is a procedure called DoubleGauss().. now that you mentioned it, I have found it in the original procedure that I used as template to build mine. Here it goes:

Function DoubleGauss(w,x) : FitFunc
    Wave w
    Variable x

    //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
    //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
    //CurveFitDialog/ Equation:
    //CurveFitDialog/ f(x) = y0+A_1*exp(-((x-x0_1)/width_1)^2)+A_2*exp(-((x-x0_2)/width_2)^2)
    //CurveFitDialog/ End of Equation
    //CurveFitDialog/ Independent Variables 1
    //CurveFitDialog/ x
    //CurveFitDialog/ Coefficients 7
    //CurveFitDialog/ w[0] = y0
    //CurveFitDialog/ w[1] = A_1
    //CurveFitDialog/ w[2] = x0_1
    //CurveFitDialog/ w[3] = width_1
    //CurveFitDialog/ w[4] = A_2
    //CurveFitDialog/ w[5] = x0_2
    //CurveFitDialog/ w[6] = width_2

    return w[0]+w[1]*exp(-((x-w[2])/w[3])^2)+w[4]*exp(-((x-w[5])/w[6])^2)
End

January 23, 2015 at 07:33 am - Permalink

Yes.. it worked.. and printed this:

W_coef[0]= {0,1,700.5,2.042,1,705,2.462}
Fit converged properly
Curve fit with data subrange:
FR_Y[1046,1096]
y= DoubleGauss(W_coef,x)
Your coefficient wave is single-precision. We recommend making it double-precision (use Redimension Waves from the Data menu).
W_coef={49.312,14.414,700.5,2.0654,17.355,705,3.1229}
V_chisq= 40.5692;V_npnts= 51;V_numNaNs= 0;V_numINFs= 0;
V_startRow= 1046;V_endRow= 1096;
W_sigma={0.167,0.651,0,0.1,0.496,0,0.121}
Coefficient values ± one standard deviation
y0 =49.312 ± 0.167
A_1 =14.414 ± 0.651
x0_1 =700.5 ± 0
width_1 =2.0654 ± 0.1
A_2 =17.355 ± 0.496
x0_2 =705 ± 0
width_2 =3.1229 ± 0.121

January 23, 2015 at 09:37 am - Permalink

Now I am pretty sure that one peak eats the other.
Changing the /H="0000000" pinned the peak positions to the original values. I am afraid you will need constraints ...
Please read the manual about fitting with constraints.

Otherwise your fitting seems to work...

HJ

January 23, 2015 at 10:33 am - Permalink


Ok.. I will read tomorrow! In the meanwhile..

Is there a way to turn automatic guesses for this FuncFit procedure?

Is there a way to suggest the peaks positions and widths for the Curvefit, in the same way I am trying with the Funcfit ?

I was checking about wavescalling... my X axis is not linear, and if I try to use the auto-scale I have different positions for my peaks!


Tks a lot...


Andre

January 23, 2015 at 01:10 pm - Permalink

andreponzoni wrote:
Is there a way to turn automatic guesses for this FuncFit procedure?

No. Automatic guesses require a great deal of knowledge of the function, and in the case of a user-defined function, Igor knows nothing about it. There have been lots of requests for the ability to define a user-defined autoguess function, but I've never had the time to figure out how to make it work smoothly.
Quote:
Is there a way to suggest the peaks positions and widths for the Curvefit, in the same way I am trying with the Funcfit ?

Make/D coefwave={...guesses...}
CurveFit fitname, kwCWave = coefwave, ...

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 23, 2015 at 01:29 pm - Permalink

In your case the initial guesses might be "simple":
The width of your lines is around 2 to 4 wave numbers, so use 3 as a start.
If your "doublet" is centered in your spectrum you can use startp+(end-startp)/3 as position1 and startp+(end-startp)*2/3 as position2.
My experience is that a too high amplitude is tolerated better than a too small one so you might try wavemax(wSpectrum[StartP,EndP]) as an initial value.

You might wonder why I suggest this rough guess if your good guesses didn't work: On one hand I think your amplitude guess was too small. On the other hand the fit routine calculates "gradients" and uses the steepest one assume it will lead to the best solution. If your guesses are very close in some parameters and off by others the steepest gradient might lead away from the optimal fit. In your case it increases amplitude of one peak, thinks its a good way and shifts the position to the center "killing" the other peak.

Play around a little bit with different starting values and see how fast the fit converges and in which parameter radius it runs or crashes.

Happy Weekend
HJ

January 24, 2015 at 01:11 am - Permalink

HJDrescher wrote:
In your case the initial guesses might be "simple":
The width of your lines is around 2 to 4 wave numbers, so use 3 as a start.
If your "doublet" is centered in your spectrum you can use startp+(end-startp)/3 as position1 and startp+(end-startp)*2/3 as position2.
My experience is that a too high amplitude is tolerated better than a too small one so you might try wavemax(wSpectrum[StartP,EndP]) as an initial value.

You might wonder why I suggest this rough guess if your good guesses didn't work: On one hand I think your amplitude guess was too small. On the other hand the fit routine calculates "gradients" and uses the steepest one assume it will lead to the best solution. If your guesses are very close in some parameters and off by others the steepest gradient might lead away from the optimal fit. In your case it increases amplitude of one peak, thinks its a good way and shifts the position to the center "killing" the other peak.

Play around a little bit with different starting values and see how fast the fit converges and in which parameter radius it runs or crashes.

Happy Weekend
HJ


Thanks so much HJ! If I increase the amplitude, it fits very well both peaks... similar to the values I got from the manual mode! Once all my waves are very similar, the initial guesses should work for them! The next step is doing the looping for batch fitting, but I shoudn't have problems with it!

I would like to leave you one more question, this time related with the wavemax function you suggested to collect the maximum intensity from the wave.

According with what I read on displayhelptopic "wavemax", the correct way should be

variable a

a= wavemax(Y_waveName[,100,200])


That should return the maximum value between 100 and 200 (i am not concerned now if is data points or X axis, I can figure that out later). But in the end I always got the following error:

Expected right parenthesis...

If I try the same without defining the range ( a=wavemax(Y_wavename) ) works pretty fine!

I also tryed without the first comma; parenthesis instead of square brackets... but this simple function doesn't work if I fix the range! Any suggestions?

Thanks again and have a nice weekend!

Andre

January 24, 2015 at 05:32 am - Permalink

You want to use

a= wavemax(Y_waveName, 100,200)

In the command reference the square brackets are indicating that the enclosed parameters are optional. The brackets are not meant to be included when the function is used.

January 25, 2015 at 12:59 pm - Permalink