Adding initial guesses to a curve fit procedure

Share icon

Wed, 05/14/2014 - 12:09 pm

I am fitting my data using a 2D Gaussian function and was wondering if there is a way to incorporate my initial guesses into my procedure. My current procedure sets a window to fit a 2D Gaussian, but it is much easier to enter initial guesses to fit a 2D Gaussian to where my peak is. I am trying to automate my procedure so it can fit 50 spectra at once to the same initial guesses.

Any help will be greatly appreciated.

Thank you,
-s0585353

Function Gaussian2Dfit(aavgNNN, w1, w3, xmin, xmax, ymin, ymax)
Wave aavgNNN
Wave w1
Wave w3
Variable xmin
Variable xmax
Variable ymin
Variable ymax
CurveFit/NTHR=0 Gauss2D aavgNNN[xmin,xmax][ymin,ymax] /X=w1[xmin,xmax] /Y=w3[ymin,ymax] /D
ModifyContour fit_aavgNNN labels=0
wave W_coef
End
Take a look at this previous topic.
http://www.igorexchange.com/node/4038
Since then Batch Curve Fit is in Igor 6.3 onwards.

You need to make a co-efficient wave containing your guesses. Relevant bit from that post is this
Make/D/N=6/O W_coef
        W_coef[0] = {0,1,0.0006,0.0155,600,0.999}
        Make/O/T/N=2 T_Constraints
        T_Constraints[0] = {"K5 > 0","K5 < 1"}
        FuncFit/H="000010"/NTHR=0 /N/Q ChapmanRichards W_coef aWave[5,44] /X=bWave /D/R /C=T_Constraints
        WAVE W_coef


The /H flag tells igor whether to hold any of the co-efficients and which ones.

May 14, 2014 at 01:28 pm - Permalink

sjr51 has it mostly right, but you are using CurveFit with the built-in Gauss2D and sjr51 gave you an answer using FuncFit.

You do, indeed, need to pre-make a coefficient wave and then use that with CurveFit:
Function Gaussian2Dfit(aavgNNN, w1, w3, xmin, xmax, ymin, ymax, amp)
    Wave aavgNNN
    Wave w1
    Wave w3
    Variable xmin
    Variable xmax
    Variable ymin
    Variable ymax
    Variable amp
    Make/D gcoefs={0, amp, (xmin+xmax)/2, xmax-xmin, (ymin+ymax)/2, ymax-ymin, 0}
    CurveFit/NTHR=0/G Gauss2D, kwCWave=gcoefs, aavgNNN[xmin,xmax][ymin,ymax] /X=w1[xmin,xmax] /Y=w3[ymin,ymax] /D
    ModifyContour fit_aavgNNN labels=0
    wave W_coef
End

I added an input parameter amp to give the initial guess of the peak height, I computed reasonable guesses for center and width from your window inputs (but you might want to pass those in directly), and I added the /G flag to tell CurveFit to use your coefficient wave to get the initial guesses.

And you might, in fact, like the convenience of Batch Curve Fit.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 14, 2014 at 04:29 pm - Permalink

the code almost works how I want it to, but there are a few glitches that need to be corrected. First one is when I run the procedure on a new spectrum (in the same experiment) I get the error message of the gcoefs wave already exists, how do I get rid of that? (All of my data will use the same gcoefs wave). Secondly, my fits are not coming out smooth, I have attached a picture to show you what I mean. The one on left is ran clicking on Analysis --> Curve Fitting --> and then typing in my initial guesses, the one on the right is ran using my procedure. I also provided a copy of my complete procedure below.
Thanks,
-s0585353

Function Gaussianfit2D(aavgNNN, w1, w3)
Wave aavgNNN
Wave w1
Wave w3
Make/D gcoefs={0, -0.05, 2049.3, 3.5, 2047.2, 3.5, 0.3}
CurveFit/NTHR=0/G Gauss2D, kwCWave=gcoefs, aavgNNN /X=w1 /Y=w3 /D
ModifyContour fit_aavgNNN labels=0
wave W_coef

if (waveExists (W_coef))
Variable integral
integral = W_coef[1]*2*pi* W_coef[3]*W_coef[5]*sqrt(1-W_coef[6]^2)
Print integral
endif
End
igorprob.png

May 15, 2014 at 12:24 am - Permalink

s0585353 wrote:
... there are a few glitches ... one is when I run the procedure on a new spectrum (in the same experiment) I get the error message of the gcoefs wave already exists, how do I get rid of that?


You have three options:

* use make/O/D gcoeffs
* use make/D/FREE gcoeffs
* use killwaves gcoeffs

See the Help information for each command to decide which will be best for your needs.

s0585353 wrote:
... Secondly, my fits are not coming out smooth, I have attached a picture to show you what I mean. ...


The resultant waves to the left fit have more points than those for the right fit. Increase the number of points in your resultant.

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAHuntsville

May 15, 2014 at 05:00 am - Permalink

jjweimer wrote:


The resultant waves to the left fit have more points than those for the right fit. Increase the number of points in your resultant.



How do I increase the number of points?

May 15, 2014 at 09:10 am - Permalink

Thanks everybody for your help the procedure seems to be working correctly. I have a few other questions, 1) I set my contour fit labels to be off but they are still appearing on the fit, am I missing something in my code? 2) I just tried to run this procedure on a spectrum that has multiple peaks, where I would only want to fit one peak at a time (if there is away to fit multiple 2D Gaussians at once on one spectrum, I would be interested in knowing how), so I input my initial guesses for one of the peaks and the result is one large 2D Gaussian over the entire spectrum and the same thing happened when I held the center of the peak (xo and yo) constant. The next thing tried was setting an area to fit the Gaussian in (I combined the 1st code I posted with the last code I posted) and this worked for the most part. If the set area was small enough, the Gaussian would be centered on my peak, but if I increased to a larger area the Gaussian would no longer be centered on my peak, the only way to have larger area is to hold the center of the peak constant, which I really don't want to do. So if you have any suggestions it would be greatly appreciated. I attached an image of the type of spectrum I am trying to fit.

Thanks,
-s0585353
igorprobmultipeaks.png

May 15, 2014 at 10:26 am - Permalink

Clearly you have overall a data set that isn't Gaussian, it's a sum of Gaussians, possibly superposed on a gradient. The built-in Gauss2D fit fits a single peak sitting on a horizontal plane. That is, it fits a Z offset but not a gradient. If you broaden the area of the fit so that it includes places where two peak overlap, you will get a fit that tries hard to interpret the whole thing as a single peak. That usually means broadening the peak to fit tails that are elevated by overlap with a neighboring peak. If the peaks sit on top of a gradient, the same thing happens- the Z offset will fit more or less the average background, and the tails may broaden somewhat to try to fit some portion of the gradient.

You can fit more than one peak if you write a user-defined function to fit a sum of Gaussian peaks. You may also need to incorporate a gradient or polynomial background.

You will have to supply initial guesses for each of the peaks in the sum, but you already have that- something is generating the windows for your original attempt.

Here are a pair of suitable fitting functions:
Function Fit2DGaussians(pw, xx, yy) : FitFunc
    Wave pw
    Variable xx, yy
   
    Variable result = pw[0]     // w[0] is the Z offset
    Make/D/N=7/FREE onePeakCoefs = 0
   
    Variable npeaks = (numpnts(pw)-1)/6
    Variable i
    for (i = 0; i < npeaks; i += 1)
        onePeakCoefs[1,] = pw[6*i+p]        // leave the Z offset at zero for the individual peaks
        result += Gauss2d(onepeakCoefs, xx, yy)
    endfor
   
    return result
end

Function Fit2DGaussiansAAO(pw, zw, xw, yw) : FitFunc
    Wave pw, zw, xw, yw
   
    zw = pw[0]      // w[0] is the Z offset
    Make/D/N=7/FREE onePeakCoefs = 0
   
    Variable npeaks = (numpnts(pw)-1)/6
    Variable i
    for (i = 0; i < npeaks; i += 1)
        onePeakCoefs[1,] = pw[6*i+p]        // leave the Z offset at zero for the individual peaks
        zw += Gauss2d(onepeakCoefs, xw[p], yw[p])
    endfor
end

Both do the same thing, the first is easier to use outside of curve fitting, for instance to fill a test matrix fill of values. The second is much faster than the first which is good for curve fitting. Both take a coefficient wave that has 6*npeaks+1 points. Point zero is the z offset. The next six points are coefficients for the first peak: amplitude, X center, X width, Y center, Y width and correlation coefficient. To fit two peaks you would need 13 points.

Hope this helps.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 16, 2014 at 04:31 pm - Permalink