Convolved Function: Curve Fitting

Share icon

Hello, I am somewhat new to IGOR programming and I am trying to make a program that will convolve a multi-exponential (2 or 3 exponential being added to each other) with a gaussian. The gaussian is known instrument response so this should be supplied by the user and then the multi-exponential component should be fit to the data set. When I use the curve fit the answers that come out are nonsense even though I tried to build in constraints for having non-negative time constants and amplitudes. Yet, it seems that IGOR will not find a solution that is even close to the real fit since I have a data set that I fit by manually and the error bars on the fit are enormous. Anyway here is the code along with the data set I am trying to fit. It is mostly modeled after the IGOR Help convolved fit program that is used as an example in the IGOR Help menu. Please let me know if you can help.

Here is the code:

#pragma rtGlobals=1 // Use modern global access method.
Function convfunc(pw, yw, xw) : FitFunc
WAVE pw, yw, xw
// pw[0] = c0 //vertical offset (Background noise) (Typically hold this parameter and use your own guess)
// pw[1] = Gaussian amplitude
// pw[2] = Gaussian Position
// pw[3] = Gaussian Width
// pw[4] = c1
// pw[5] = c2
// pw[6] = tau1
// pw[7] = tau2
// Gaussian parameters can set these to be held during the command line (/H="hhhhh..."
// where 1 denotes hold and 0 denotes search in the ith constraint on pw[i])

//Resolution sets the degree at which the exponential will be over-sampled with regard
//to fit parameters, increasing this number
Variable resolutionFactor = 10
Variable dT = min(pw[3]/resolutionFactor, pw[5]/resolutionFactor)
Variable nExpWavePnts = round((resolutionFactor*pw[5])/dT)*2 + 1

Make/D/O/N=(nExpWavePnts) expwave

Variable nYPnts = max(resolutionFactor*numpnts(yw), nExpWavePnts)
Make/D/O/N=(nYPnts) yWave

setscale/P x -dT*(nExpWavePnts/2),dT,expwave
setscale/P x xw[0],dT, yWave

//Heavside function of a multiexponential (x>=0) = 0 when x<0, and =1 when x>=0
expwave = (x>=0)*(pw[4]*exp(-xw/pw[5]))+pw[6]*exp(-xw/pw[7]))

Variable sumexp
sumexp = sum(expwave, -inf,inf)
expwave /= sumexp

ywave =pw[0] + pw[1]*exp(-(x-pw[2])^2/pw[3])
convolve/A expwave, ywave
yw=yWave(xw[p])
End

Also here is the command line that I have used to implement the function when analyzing the data set:

FuncFit/H="11110000"/NTHR=0 ConvFunc pw yw /X=xw /D /C=T_Constraint
9 iterations with no decrease in chi square
--Curve fit with constraints--
Active Constraint: Desired: K5>-0 Achieved: K5=-4.44089e-16
Duplicate/O fit_yw,WMCF_TempAutoXWave
ConvFunc(pw,fit_yw,WMCF_TempAutoXWave)
KillWaves/Z WMCF_TempAutoXWave
Your coefficient wave is single-precision. We recommend making it double-precision (use Redimension Waves from the Data menu).
pw={99,30000,125,180,0.13957,-4.4409e-16,1.417e+05,151.3}
V_chisq= 1.53992e+10;V_npnts= 3126;V_numNaNs= 0;V_numINFs= 0;
V_startRow= 0;V_endRow= 3125;
W_sigma={0,0,0,0,6.99e+07,5.65,4.3e+09,6.16}
Coefficient values ± one standard deviation
K0 =99 ± 0
K1 =30000 ± 0
K2 =125 ± 0
K3 =180 ± 0
K4 =0.13957 ± 6.99e+07
K5 =-4.4409e-16 ± 5.65
K6 =1.417e+05 ± 4.3e+09
K7 =151.3 ± 6.16
**** Singular matrix error during curve fitting ****
There may be no dependence on these parameters:
pw[5]
Fit converged properly

Here are the problems I can not seem to fix.
1.) Singular matrix error during the fit, I suspect I should not be using pw, xw, yw as the input waves. Should I be using a temporary wave as the input in the command line for the parameters and the two waves to be convolved.
2.) The constraints I have set seem to continually be violated no matter what I do. I especially do not want negative amplitudes and time constants, these are nonsense in the sense of defining components involved in the signal that needs to be analyzed. This also seems to be giving extremely large error bars.
3.) I can not seem to find a way with three exponential that define all of the parameters to be unique in the sense that they are mathematically distinguishable during the fit. It seems that even when using constraints to say that c1 4.) I additionally after this need to modify this to perform a second convolution, so I will need to open up some of the fit parameters. Right now parameters k0-k3 are constants that do not need to be fit. It seems as if there is a better way to enter these in than to use up space in the fit function, since I know there is a limit of 10 parameters that can be fit in parallel so I need a way to enter in the known gaussian. Right now I am entering the parameters in the parameter wave for the fit and then I am just holding these parameters during the curve fit.

Anything would help I am pretty lost right now, I have been attempting to use the debugger to find the problems, but it seems that the way I am inputting these waves may be fundamentally incorrect.
I also I attached the data set and constraints I am attempting to use as well as the parameters. The initial guesses I am using are the values I found to fit the graph manual by just iteratively searching and calculating residual values.

The relevant waves are named IRF and Decay. IRF is the simulated instrument response. Decay is the actual data set that I need to convolution to fit to. Of course pw is the parameter wave and expwave is the double exponential to be fit. (I choose a known (double exponential*gaussian) fit just to test out the program.) Please let me know if you need any more information I know the IGOR attachment I am sending you is a bit messy, so please let me know if there is any more information you need to decipher this mess. The graph shows mostly what I am attempting to do though.

Best,
Greg Fahs


 ConvolvedFitProcedure1.4.pxp (726.42 KB)
Out of curiosity, have you actually measured the system response to be a gaussian? In real equipment there can be a variety of other effects which skew the instrument response, or add another component (such as afterpulsing). Just for improved accuracy, I would suggest that you measure the instrument response and use that data for the convolution, instead of a generic gaussian curve.

March 26, 2012 at 11:02 am - Permalink

proland wrote:
Out of curiosity, have you actually measured the system response to be a gaussian? In real equipment there can be a variety of other effects which skew the instrument response, or add another component (such as afterpulsing). Just for improved accuracy, I would suggest that you measure the instrument response and use that data for the convolution, instead of a generic gaussian curve.



Generally we do actually measure the instrument response. In our case it is actually fairly Gaussian and for the purpose of getting the program working I am actually using a simulated instrument response that the analysis program Sym-Pho time gives me. This program actually already does convolution fits, but it does not perform a double convolution fit. The end result of this program should do a curve fit of
Fit(Convolve[Exp2*Convolve[Exp1*IRF]]). Where exp1 and exp2 are separate multiexponentials (usually of the form (x>=0)(c1*exp(-x/t1)+c2*exp(-x/t2)...)). It seems though that I can not even get this program working on a file where I supply guesses that are almost exact to the real functional form. Although in this case using my simulated IRF seems to be working well for the testing phase (we need to use this to fit a few hundred files so I am trying to use a decay curve that is well defined to test this.) Regardless it does not seem to be the inadequacy of the IRF I am using to fit the data, although you are right typically these do have some afterpulse effect, but it is minimized on our setup or at least for this data (at least so I believe). Anyway thank you.

March 26, 2012 at 11:20 am - Permalink

Fitting multiple exponentials can be very tricky. They will cause problems like what you describe if the time constants are not quite distinct. In my experience, even with clean data they time constants need to be different by a factor of three or more. Constraints won't help with that.

What are the initial guesses for the exponential terms?

Perhaps, since you say you have excellent estimates of the instrument parameters, you can verify correctness of your fitting function by holding the time constants of the exponential terms and doing the fit.

By the way, the experiment file you attached doesn't include the procedure file containing your fit function.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

March 26, 2012 at 01:02 pm - Permalink

johnweeks wrote:
Fitting multiple exponentials can be very tricky. They will cause problems like what you describe if the time constants are not quite distinct. In my experience, even with clean data they time constants need to be different by a factor of three or more. Constraints won't help with that.

What are the initial guesses for the exponential terms?

Perhaps, since you say you have excellent estimates of the instrument parameters, you can verify correctness of your fitting function by holding the time constants of the exponential terms and doing the fit.

By the way, the experiment file you attached doesn't include the procedure file containing your fit function.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com



Thanks for your quick reply John. I was wondering a few things further about this. It seems from what you are saying that if the exponentials are not extremely different the curve fit algorithm seems to fail to find good solutions (probably because there is a very large set of solutions where it will run into multiple minima on the chi squared surface). This would seem to be a problem with the search algorithm that IGOR's curve fit procedure algorithm uses. Now here is my real question. The end goal of this program is to be able to fit functions that are of the form
Convolve(Exp1*Convolve(Exp2*IRF))
where Exp1=a*exp(-x/t1)
Exp2=b*exp(-x/t2)+c*exp(-x/t3)
sometimes Exp2 is actually a triple exponential, but for simplicity sake I will present it as a double exponential linear combination.
The IRF is a known Gaussian function in this case 30000*exp(-(x-125)/180) (This is actually not the instrument response, but rather just a function that minimized a sum of residuals this is how I have been fitting these manually by having a graph of the residuals at each point and the minimizing each area locally since there is obviously long and short decay components in the double/triple exponential term.)

The coefficients I found manually were
a=0.002;b=0.0955;t1=682;t2=181;

The rest are constants. You specifically mentioned that the Curve fit algorithm has problems dealing with exponential constants that are this similar. As it would be these are actually the most dissimilar coefficients at least comparatively to most of the other files. Is this a somewhat lost cause if I keep attempting to do this by using constraints? It seems that even when I had this program working properly the curve fit would simply spit back out my initial guesses (which were good guesses in the first place) with enormous uncertainties (i.e 10^14 numbers in this range... in the range of machine precision to be honest which makes me even more nervous). Is there any other algorithm that you know of that would do this better, or should I attempt to continue to use the all at once curve fit, but in the procedure build in some if statements that are bit more strict as a regularizer than building a constraint wave (non-negativity constraints and inequalities on the exponential coefficients and time constants). I guess at this point I am just looking for advice since it seems like for a curve fit where the coefficients are so close that I will have to do something extra to get values that I can actually rely on for my analysis. Thank you to whoever decides to help me with this.

Also, as for your last piece of advice. The program we have actually works very well for the type of file I have in the Igor file I attached to this, but generally no the time-constants are not well defined before hand, this was just an example that was well-defined for testing the first trials of this program. The time constants actually the most important piece of information we are looking for. The program we currently have will not perform the double convolution that have I described in the beginning so we have no a priori information on the exact time coefficients for either Exp1 or Exp2, we really just have ranges on what is physically expected of these types of systems. We know for example that we have decay constants that around 650ish and 150ish and so on, we really do not know the exact time constant before hand and this is actually the most vital piece of information, actually the amplitude information is actually much less important if there was a reasonable way that you can think of to define these before hand. I honestly do not see a way to do that with any reliability since again we have no a priori information on this either. It seems that I would either have to find another optimization algorithm that is more sensitive for multi-exponentials or that I may need to cut my loses and fit all of these by hand. Any advice you have would be very much appreciated.

Best,
Greg
UMass Amherst Chemistry Department

March 27, 2012 at 11:57 am - Permalink

The problem with exponentials is not an Igor problem, it's mathematical. When the time constants are similar, there just isn't a lot of difference in shape between a single exponential and double, so the fit doesn't have much difference in chi-square to know which way to go. When you add a bit of noise, and probably when you blur everything with a convolution, the differences are even smaller.

Your time constants may be sufficiently different, but it's hard to say.

Singular matrix errors and the very large sigma are both signs of "identifiability problems". That is, two or more coefficients are strongly correlated such that if you change one, you can make a compensating change in another. This results in a chi-square surface with a very flat bottom that can't be optimized with any great certainty. You can diagnose identifiability problems by forming the correlation matrix as described here:

DisplayHelpTopic "Correlation Matrix"

Off-diagonal terms with absolute values around 0.9 are probably OK. Values around 0.99 are a sign of trouble but are sometimes OK. Values around 0.999 are very bad.

Have you tried your fitting function outside of a curve fit to make sure it is actually doing what you expect? If you do that, and play around with the time constants you may see that time constants with similar values are hard to distinguish from a single exponential.

Sampling can also influence success in fitting multiple exponentials. Sampling at the same rate all the time is actually not a good strategy, as it undersamples events with short lifetime. It is best to sample heavily during the lifetime of the short time constant, and more slowly later. Of course, if your time constants are similar, it is hard to define "early" and "late".

There are a couple of algorithms around that claim to be able to separate multiple exponentials better than a chi-square fit. I have never had the time to really investigate or to incorporate them into Igor. Try Googling "double exponential fit" or "multiple exponential fit". My recollection is that they didn't include convolution with a Gaussian or anything else.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

March 27, 2012 at 04:14 pm - Permalink

Gfahs wrote:
Hello, I am somewhat new to IGOR programming and I am trying to make a program that will convolve a multi-exponential (2 or 3 exponential being added to each other) with a gaussian. The gaussian is known instrument response so this should be supplied by the user and then the multi-exponential component should be fit to the data set. When I use the curve fit the answers that come out are nonsense even though I tried to build in constraints for having non-negative time constants and amplitudes. Yet, it seems that IGOR will not find a solution that is even close to the real fit since I have a data set that I fit by manually and the error bars on the fit are enormous. Anyway here is the code along with the data set I am trying to fit. It is mostly modeled after the IGOR Help convolved fit program that is used as an example in the IGOR Help menu. Please let me know if you can help.

Here is the code:

#pragma rtGlobals=1 // Use modern global access method.
Function convfunc(pw, yw, xw) : FitFunc
WAVE pw, yw, xw
// pw[0] = c0 //vertical offset (Background noise) (Typically hold this parameter and use your own guess)
// pw[1] = Gaussian amplitude
// pw[2] = Gaussian Position
// pw[3] = Gaussian Width
// pw[4] = c1
// pw[5] = c2
// pw[6] = tau1
// pw[7] = tau2
// Gaussian parameters can set these to be held during the command line (/H="hhhhh..."
// where 1 denotes hold and 0 denotes search in the ith constraint on pw[i])

//Resolution sets the degree at which the exponential will be over-sampled with regard
//to fit parameters, increasing this number
Variable resolutionFactor = 10
Variable dT = min(pw[3]/resolutionFactor, pw[5]/resolutionFactor)
Variable nExpWavePnts = round((resolutionFactor*pw[5])/dT)*2 + 1

Make/D/O/N=(nExpWavePnts) expwave

Variable nYPnts = max(resolutionFactor*numpnts(yw), nExpWavePnts)
Make/D/O/N=(nYPnts) yWave

setscale/P x -dT*(nExpWavePnts/2),dT,expwave
setscale/P x xw[0],dT, yWave

//Heavside function of a multiexponential (x>=0) = 0 when x<0, and =1 when x>=0
expwave = (x>=0)*(pw[4]*exp(-xw/pw[5]))+pw[6]*exp(-xw/pw[7]))

Variable sumexp
sumexp = sum(expwave, -inf,inf)
expwave /= sumexp

ywave =pw[0] + pw[1]*exp(-(x-pw[2])^2/pw[3])
convolve/A expwave, ywave
yw=yWave(xw[p])
End

Also here is the command line that I have used to implement the function when analyzing the data set:

FuncFit/H="11110000"/NTHR=0 ConvFunc pw yw /X=xw /D /C=T_Constraint
9 iterations with no decrease in chi square
--Curve fit with constraints--
Active Constraint: Desired: K5>-0 Achieved: K5=-4.44089e-16
Duplicate/O fit_yw,WMCF_TempAutoXWave
ConvFunc(pw,fit_yw,WMCF_TempAutoXWave)
KillWaves/Z WMCF_TempAutoXWave
Your coefficient wave is single-precision. We recommend making it double-precision (use Redimension Waves from the Data menu).
pw={99,30000,125,180,0.13957,-4.4409e-16,1.417e+05,151.3}
V_chisq= 1.53992e+10;V_npnts= 3126;V_numNaNs= 0;V_numINFs= 0;
V_startRow= 0;V_endRow= 3125;
W_sigma={0,0,0,0,6.99e+07,5.65,4.3e+09,6.16}
Coefficient values ± one standard deviation
K0 =99 ± 0
K1 =30000 ± 0
K2 =125 ± 0
K3 =180 ± 0
K4 =0.13957 ± 6.99e+07
K5 =-4.4409e-16 ± 5.65
K6 =1.417e+05 ± 4.3e+09
K7 =151.3 ± 6.16
**** Singular matrix error during curve fitting ****
There may be no dependence on these parameters:
pw[5]
Fit converged properly

Here are the problems I can not seem to fix.
1.) Singular matrix error during the fit, I suspect I should not be using pw, xw, yw as the input waves. Should I be using a temporary wave as the input in the command line for the parameters and the two waves to be convolved.
2.) The constraints I have set seem to continually be violated no matter what I do. I especially do not want negative amplitudes and time constants, these are nonsense in the sense of defining components involved in the signal that needs to be analyzed. This also seems to be giving extremely large error bars.
3.) I can not seem to find a way with three exponential that define all of the parameters to be unique in the sense that they are mathematically distinguishable during the fit. It seems that even when using constraints to say that c1 4.) I additionally after this need to modify this to perform a second convolution, so I will need to open up some of the fit parameters. Right now parameters k0-k3 are constants that do not need to be fit. It seems as if there is a better way to enter these in than to use up space in the fit function, since I know there is a limit of 10 parameters that can be fit in parallel so I need a way to enter in the known gaussian. Right now I am entering the parameters in the parameter wave for the fit and then I am just holding these parameters during the curve fit.

Anything would help I am pretty lost right now, I have been attempting to use the debugger to find the problems, but it seems that the way I am inputting these waves may be fundamentally incorrect.
I also I attached the data set and constraints I am attempting to use as well as the parameters. The initial guesses I am using are the values I found to fit the graph manual by just iteratively searching and calculating residual values.

The relevant waves are named IRF and Decay. IRF is the simulated instrument response. Decay is the actual data set that I need to convolution to fit to. Of course pw is the parameter wave and expwave is the double exponential to be fit. (I choose a known (double exponential*gaussian) fit just to test out the program.) Please let me know if you need any more information I know the IGOR attachment I am sending you is a bit messy, so please let me know if there is any more information you need to decipher this mess. The graph shows mostly what I am attempting to do though.

Best,
Greg Fahs


Hello,

I came across the same convolution problem. Did you fix the problem and can you upload the new convolute fit program if possible.

Thanks.

January 24, 2018 at 11:30 am - Permalink