Expanding PeakFunctions2 functionality
jgrinias
Fri, 03/23/2012 - 08:18 am
Long time user, very bad beginner programmer. I am currently trying to create a new fit function using a two-time constant exponentially modified Gaussian (see attached document). Initially, I thought it would be simple enough to build off some of the code already present for the EMG functionality in PeakFunctions2.ipf, but there are many more functions related to "Function EMGFunction(w,x)" then I originally thought.
It seems that the process shouldn't be too much more than creating a new fit function and then finding a way to expand the reported data table to include one more variable, but I'm not sure the best way to proceed. Any information would be greatly appreciated.
Thank you,
Jim
DisplayHelpTopic "Adding Your Own Peak Function"
It is a somewhat more involved process than I wish it was...
All those extra functions are required to give the Multipeak Fit 2 machinery information it needs about the peak function so that it will work with MPF2. If you don't need MPF2, you could just write the peak function as a simple user-defined fit function.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
March 23, 2012 at 02:22 pm - Permalink
In trying to fit a peak using the "New Fit Function" under curve fitting, the program seems to run in a loop without ever getting any data out. Is this what you meant or is there some other method to simply allow a new function to be one the pull-down options using the MPF (either 1.4 or 2)?
March 26, 2012 at 08:13 am - Permalink
DisplayHelpTopic "All-At-Once Fitting Functions"
Once you have successfully created a peak fitting function, then to use it with MPF2 you need to add the machinery documented in the help for MPF2.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
March 26, 2012 at 10:09 am - Permalink
Function DoubModGauss(pw, yw, xw) : FitFunc WAVE pw, yw, xw // pw[0] = x0 // pw[1] = w // pw[2] = h // pw[3] = s1 (tau1) // pw[4] = s2 (tau2) Variable C1=(((xw-pw[0])/pw[1])-(pw[1]/pw[3])) Variable C2=(((xw-pw[0])/pw[1])-(pw[1]/pw[4])) Variable D1=(exp(((pw[1]^2)/(2*(pw[3]^2)))-((xw-pw[0])/pw[3]))*(1-erf(-C1/sqrt(2))) Variable D2=(exp(((pw[1]^2)/(2*(pw[4]^2)))-((xw-pw[0])/pw[4]))*(1-erf(-C2/sqrt(2))) yw=(((pw[2]*pw[1])/(pw[3]-pw[4]))*(sqrt(3.14159/2))*(D1-D2)) endWhen I try running a curve fitting on my peak using this equation, it still seems to get stuck in a long loop. I think it's because it's initially making a long guess based on each data point in the initial wave. Do I need to create some sort of guess wave (use 'Make' to make a pw initial wave)?
March 26, 2012 at 07:15 pm - Permalink
http://www.igorexchange.com/node/2954
Essentially, the ExpModGauss function that is included in the MultiPeakFitting2 package doesn't quite get me a good fit because I know there are two time constants acting differently and the equation I describe in the above PDF is the best way I can find in the literature to discern the two time constants specific to my problem.
March 26, 2012 at 07:42 pm - Permalink
That's never going to work. Where in the calculation of yw is xw used (it's not)? D1 and D2 are variables, only containing a single value, not a whole range of xvalues. Here is an example for an all-at-once single gaussian:
Function mygauss(pw, yw, xw):fitfunc Wave pw, yw, xw yw = pw[0] + pw[1] * exp(-((xw - pw[2])/pw[3])^2) EndHere yw is calculated for every single point of xw.
March 27, 2012 at 01:00 am - Permalink
Whenever you develop a function like this, you should test it outside of curve fitting. Just use it directly to create model data sets like this:
SetScale/I x minX, maxX, junkx
junkx = x
MyAAOFunc(coefs, junk, junkx)
Display junk vs junkx
This assumes a pre-existing coefficient wave called "coefs". You can change /N=100 to any number you like. The SetScale and assignment of x is just a convenient way to fill a wave full of evenly-space X values. You will need to choose actual values for minX and maxX. Naturally, you only need the Display command once.
If you had done this previously, you would have found that the shape of the curve wasn't right.
An alternative to intermediate waves would be a loop that references xw[i] instead of xw. The loop would also have to assign a single value to yw[i] in each iteration. I would not recommend that as loops in Igor tend to be slow compared to wave assignments.
To learn more about wave assignments:
DisplayHelpTopic "Waveform Arithmetic and Assignment"
The most efficient method that would preserve your structure would probably be to use MatrixOP for the intermediate results.
DisplayHelpTopic "MatrixOP"
When you have digested all that, the next step would be to use free waves for any intermediate results, for an improvement in performance:
DisplayHelpTopic "Free Waves"
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
March 27, 2012 at 09:28 am - Permalink
Function DoubModGauss(pw, yw, xw) : FitFunc
WAVE pw, yw, xw
// pw[0] = x0
// pw[1] = w
// pw[2] = h
// pw[3] = s1 (tau1)
// pw[4] = s2 (tau2)
// C1=(((xw-pw[0])/pw[1])-(pw[1]/pw[3]))
// C2=(((xw-pw[0])/pw[1])-(pw[1]/pw[4]))
// D1=(exp(((pw[1]^2)/(2*(pw[3]^2)))-((xw-pw[0])/pw[3]))*(1-erf(-C1/sqrt(2)))
// D2=(exp(((pw[1]^2)/(2*(pw[4]^2)))-((xw-pw[0])/pw[4]))*(1-erf(-C2/sqrt(2)))
//yw=(((pw[2]*pw[1])/(pw[3]-pw[4]))*(sqrt(3.14159/2))*(D1-D2))
yw=(((pw[2]*pw[1])/(pw[3]-pw[4]))*(sqrt(3.14159/2))*((exp(((pw[1]^2)/(2*(pw[3]^2)))-((xw-pw[0])/pw[3]))*(1-erf(-(((xw-pw[0])/pw[1])-(pw[1]/pw[3]))/sqrt(2)))-(exp(((pw[1]^2)/(2*(pw[4]^2)))-((xw-pw[0])/pw[4]))*(1-erf(-(((xw-pw[0])/pw[1])-(pw[1]/pw[4]))/sqrt(2)))))))
end
Function Test()
Make/N=7025 junk,junkx
SetScale/I x 0, 7025, junkx
junkx = x
DoubModGauss(coef, junk, junkx)
Display junk vs junkx
End
The equation works to generate a basic curve. This generates the plot based on the values for 'pw' that I have selected in 'coefs'. How do I then apply the equation to take a y-wave that I input (and it's calculated x-wave for specific scaling) and get the 5 variables that I am looking for back out?
April 2, 2012 at 08:07 am - Permalink
You use it as a fit function with the FuncFit command. One way to do that is using the Curve Fit dialog. The FuncFit command takes care of generating appropriate inputs for your fit function and running the function.
To use it as a peak shape in Multipeak Fit 2, you need to add a fair bit of machinery as described in the Multipeak Fit 2 help file.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
April 2, 2012 at 01:22 pm - Permalink
I have attached an image of the comparison between the fit to generated data (giving coefficients that are essentially equal to the chosen parameters as expected) and then to some of the actual data needed to fit. The general shape of the fit is actually worse than what is used when applying the ExpModGauss function in MPF2. Does this fit program use the convolution between a Gaussian and an exponential rather than just a single fit equation as I am trying? Or, is it some relation to the initial guess generated using MPF2 being better than what I am choosing in my initial coefficient wave?
Additionally, are there issues related to wave scaling in using the curve fit? The data presented here was not scaled for simplicity sake, but most of the data that I will be characterizing is scaled to a recorded data acquisition rate when it is imported into Igor for analysis (thankfully written by graduate students past). When trying to run this peak fit on some of that data, I usually received even worse fits or a NaN return on some x-values.
April 3, 2012 at 11:46 am - Permalink
I am guessing that the fit that doesn't go through the right-hand tail is the one that used your fit function. It appears that your fit function simply isn't a good model for the data. If the fit is successful (that is, it converges to a solution without an error message) and the fit appears good to your eye (the peak itself is well fit, and the extreme tail is well fit) then it is not a failure of the initial guesses. In a successful fit, the initial guesses will have only a minor effect on the eventual solution.
Caveat: if the fit has multiple local minima, that is, multiple possible solutions, then different initial guesses can arrive at different solutions. But that should be evident from arriving at different values of the coefficients.
If you enter a data wave with wave scaling as the Y Data wave, and you don't use the /X= flag when you invoke FuncFit, then the X values will be pulled from the wave scaling of the Y wave automatically. You don't have to do anything to achieve that. FuncFit will then feed those X values to your fitting function via the input wave xw. FuncFit synthesizes waves to use as the inputs to your fitting function when it runs. They are not the same waves as your data waves.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
April 4, 2012 at 09:47 am - Permalink
An equation that had been developed empirically has been attached in the PDF as my next attempt at a fit. Now that I feel more comfortable with some of the main features of Igor programming (I sat down with the "programming training module" to try and learn more), I tried to put together a new fit function:
Function ExpGaussHybrid(pw, yw, xw) : FitFunc Wave pw,yw,xw //pw[0] = x0 //pw[1] = w //pw[2] = h //pw[3] = s (ExpTau) if (((2*(pw[1]^2))+(pw[3]*(xw-pw[0])))>0) yw=(pw[2]*exp((-(xw-pw[0])^2)/((2*(pw[1]^2))+(pw[3]*(xw-pw[0]))))) else yw=0 endif EndOf course, no results are able to be returned. What I can see happening here is something is wrong with the conditional, mainly because I seem to need quite a few more right parentheses than I had counted. I understand simple if-endif statements using dummy variables based on 0 and 1, but is it more complex when using values that are going to be determined by the fit function in one of the conditionals (should I be using some other functionality for this type of iterative loop)?
April 4, 2012 at 11:59 am - Permalink
for (i = 0; i < numpnts(yw); i += 1)
if ((((2*(pw[1]^2))+(pw[3]*(xw[i]-pw[0])))>0)
yw[i]=(pw[2]*exp((-(xw[i]-pw[0])^2)/((2*(pw[1]^2))+(pw[3]*(xw[i]-pw[0])))))
else
yw[i]=0
endif
endfor
I'm sure that this could be made faster using one or two matrixOP expressions. What I have shown here is conceptually easier to understand. Another option would be to make a fairly unreadable one-liner using the conditional operator, like
It is also likely that the expression in the conditional is also in the fit function expression, and should be pre-calculated and assigned to a variable. Then use the variable in the places where the expression appears.
April 4, 2012 at 05:41 pm - Permalink