Fit.doFit with more than 6 parameters

Hi. So I aim to use the macro language to fit a FLIM curve with 8 parameters (which I cannot reduce), since I saw in the Curve Fitter that polynomial functions can work with up to 9 (a-i). Problem is, in a macro only 6 parameters (a-f) get recognized when trying to execute the code. Is there any way for me to incorporate more of them in an equation? Any help would be greatly appreciated.

This is working:

// This macro demonstrates how to fit a 
// polynomial to a set of data points.

// Create 25 coordinate pairs
n = 25;
xpoints = newArray(n);
ypoints = newArray(n);
xmin = -2;
xmax = 2;
xinc = (xmax - xmin) / (n - 1);
x = xmin;
for (i = 0; i < n; i++) {
    xpoints[i] = x + random / 6;
    ypoints[i] = x * x + random / 6;
    x += xinc;
}

Fit.doFit("8th Degree Polynomial", xpoints, ypoints);
Fit.plot;

Yes, thank you, but the problem is that I can’t use a polynomial function. I’m looking to use a custom function with 8 parameters, and the macro doesn’t seem to recognise more than 6 for any function that isn’t the built-in higher order polynomials. But since it works for those, I still hope there is a way.

Yes, it seems as noted here:

https://imagej.nih.gov/ij/macros/examples/UserDefinedCurveFits.txt

However, maybe @Wayne can answer that.

1 Like

Hi everyone,
as far as I can say, custom fits tend to be problematic already with six parameters; for more parameters I fear that the chances are quite slim that it would yield good results. ImageJ uses the simplex (Nelder-Mead algorithm) for minimizing. There are even cases with two or three parameters where the simplex method not converge to the correct values:
https://docs.mantidproject.org/v3.7.1/concepts/FittingMinimizersComparisonDetailed.html#minimizers-unweighted-comparison-in-terms-of-accuracy
Eight parameters might be possible in benign cases with Levenberg-Marquardt least-squares fitting, but at least for ImageJ 1 one would have find for a license-free Levenberg-Marquardt code, or write one (it should be a version that does not need analytic equations for the derivatives). In a very quick and superficial search I did not find any license-free code. If someone should know such code or at least detailed pseudo-code, it could be implemented in ImageJ.
But I don’t know whether Levenberg-Marquardt is as good as a simplex in case the parameters have very different orders of magnitude. E.g., the simplex minimizer in ImageJ still works if one parameter has a reasonable value between 0 and 1 and another parameter around 1e20 or 1e-20. E.g. I could get this one working with the standard ImageJ curve fitter data and suitable starting points:

y = 1e20*a*exp(-(x-b)*(x-b)/(1e-40*c*c))

It yields the correct result:
a = 2.51655E-19, b = 5.50923, c = 2.49040E20
(it is the “Gaussian no offset” with some factors to get large parameter values).
Simple implementations of Levenberg-Marquardt that use finite differences for the derivatives with fixed delta-x (or a minimum value of delta-x) would probably fail in this case.

Coming back to the 8-parameter fit: What is the function that you want to fit? Maybe it can be reduced to a 6-parameter fit when implemented in a plugin. ImageJ can eliminate up to two parameters, a fixed offset and either a slope or a prefactor; it calculates these via linear regression.
Eliminating two parameters is also the reason why ImageJ offers polynomial fits up to 8th order. Nevertheless, the high orders don’t always work correctly; for high orders the x values should straddle zero to ensure a correct result, and for 8th order the x values should be roughly symmetric around zero. The 8th-order polynomial fit in the above example macro fulfills this and works well. I have tried it and compared it with Mathematica - the fit parameters are essentially the same (some deviations in the 5th significant digit. Mathematica does it as a linear fit problem, so it is certainly more accurate; except for the elimination of two parameters ImageJ handles it as a general case which may be linear or nonlinear, which makes fitting more difficult).

–Michael

3 Likes

Hello, Michael, and thank you very much for the lengthy reply.

I am trying to fit a fluorescent signal to a biexponential decay convolved by a double Gaussian IRF as determined by the arrival times of photons for every pixel in an image time series (normally, those values would obey a normal Gaussian distribution, but an unknown artifacts causes them to be separated into two populations). In this case, to reduce the number of parameters, both the biexponential and the bigaussian are normalized:

Bigaussian as the sum of two Gaussians:

BiG(x;a,e,f,g,h) = G1(x;e,g) + a*G2(x;f;h)

where e and f (e>f) are the modes, s1 and s2 are the standard deviations and p is the ratio of the second to the first peak (in a normalized form a fraction).

The biexponential takes the form:

BiE(x;b,c,d) = E1(x;b,c) + E2(x;b,d) = b*exp(-(x/c))+(1-b)*exp(-(x/d))

where b is the fractional coefficient of one decay (dependent on stoichiometry of the two fluorescent compounds, 0<b<1) and c and d are the exponential decay constants which I am looking to calculate.

If I’m not mistaken, the convolution of the sum of two exponential decays by the sum of two Gaussians is the sum of the convolutions of each exponential with each Gaussian, since it’s the closed form solution of the convolution integral thereof:

Exp(x;k) X Gauss(x;m,s) = 1/2*exp(((s*s)/(2*k*k))-((x-m)/k))*(1+erf(((s/k)-((x-m)/s))/sqrt(2)))

yielding:

y= 1/2*b*exp(((e*e)/(2*c*c))-((x-g)/c))*(1+erf(((e/c)-((x-g)/e))/sqrt(2)))+1/2*(1-b)*exp(((e*e)/(2*d*d))-((x-g)/d))*(1+erf(((e/d)-((x-g)/e))/sqrt(2))) + a*(1/2*b*exp(((f*f)/(2*c*c))-((x-h)/c))*(1+erf(((f/c)-((x-h)/f))/sqrt(2)))+1/2*(1-b)*exp(((f*f)/(2*d*d))-((x-h)/d))*(1+erf(((e/d)-((x-h)/f))/sqrt(2))))

I can at least reduce the parameters to 7 by locating the peak of the signal curve, cutting it there and setting h to 0:

y= 1/2*b*exp(((e*e)/(2*c*c))-((x-g)/c))*(1+erf(((e/c)-((x-g)/e))/sqrt(2)))+1/2*(1-b)*exp(((e*e)/(2*d*d))-((x-g)/d))*(1+erf(((e/d)-((x-g)/e))/sqrt(2))) + a*(1/2*b*exp(((f*f)/(2*c*c))-(x/c))*(1+erf(((f/c)-(x/f))/sqrt(2)))+1/2*(1-b)*exp(((f*f)/(2*d*d))-(x/d))*(1+erf(((e/d)-(x/f))/sqrt(2))))

Am I on the right track and is there any chance to get ImageJ to eliminate at least one parameter here? Thank you very much in advance as well as for the help you’ve already provided.

Hi Anatas,
as far as I can see, your functions have no constant offset. The equations are quite lengthy, but if you can rewrite them to have a constant multiplier, i.e., they are of the form y = b*f(x, other parameters), you can eliminate that multiplier ‘b’ in a plugin, which will make fitting easier.
The key is the CurveFitter’s setOffsetMultiplySlopeParams(int offsetParam, int multiplyParam, int slopeParam) method.
https://github.com/imagej/imagej1/blob/7445d833ff2814215e1803c685f94617dcc6d901/ij/measure/CurveFitter.java#L340
Your code might look roughly like the following:

public class ... implements UserFunction {
  ...
    CurveFitter cf = new CurveFitter(xData, yData);
    cf.setOffsetMultiplySlopeParams(-1, 1, -1); //if parameter 1 ('b') is a multiplier for the whole function
    cf. doCustomFit(this, /*numParams=*/6, "y=... for display only",
        initialParams, initialParamVariations, /*showSettings=*/false);
    if (cf.getStatus() != Minimizer.SUCCESS) {
        ... error message with cf.getStatusString()
    }
    double[] fitParameters = cf.getParams();
 ... 
    
  /** Define your fit function here. Fit parameters: a = params[0], b = params[1], etc. */
  public final double userFunction(double[] params, double x) {
    return params[1]*Math.exp(...)*...; //params[1] must be multiplier for everything
  }
}

The initialParams array should be initial guesses for the fit parameters, and the initialParamVariations may hold numbers for roughly how large the initial steps for these parameters are. By default, ImageJ uses 10% of the initial parameter value, or a step of 0.01 if the initial value is zero. You can use null for the initialParamVariations to use these defaults.

It won’t help if you have no parameter that is a common multiplier to the whole fit function (i.e., if you do not fit the initial magnitude of the decaying function, because it is already normalized).

–Michael