Is it possible to select the parameters in the approximation equation using MATLAB?

I have a set of empirical data that is displayed on a graph (figure). There is an equation describing the dependence y=a*exp(b*x+c*x^2), you need to select the constants a, b, c. Is it possible to do this using MATLAB tools?

Answer 1

[John D'Errico]

Yes. Use the curve fitting toolbox.

help fittype

 fittype - Fit type for curve and surface fitting
    This MATLAB function creates the fittype object aFittype for the model
    specified by libraryModelName.

    Syntax
      aFittype = fittype(libraryModelName)

      aFittype = fittype(expression)
      aFittype = fittype(expression,Name,Value)

      aFittype = fittype(linearModelTerms)
      aFittype = fittype(linearModelTerms,Name,Value)

      aFittype = fittype(anonymousFunction)
      aFittype = fittype(anonymousFunction,Name,Value)

    Input Arguments
      libraryModelName - Library model to fit
        character vector | string scalar
      expression - Model to fit
        character vector | string scalar
      linearModelTerms - Model to fit
        cell array of character vectors | string array
      anonymousFunction - Model to fit
        anonymous function

    Name-Value Arguments
      coefficients - Coefficient names
        character vector | string scalar | cell array of character vectors |
        string array
      dependent - Dependent variable name
        y (default) | character vector | string scalar
      independent - Independent variable names
        x (default) | character vector | string scalar |
        cell array of character vectors | string array
      options - Fit options
        fitoptions
      problem - Problem-dependent (fixed) parameter names
        character vector | string scalar | cell array of character vectors |
        string array

    Output Arguments
      aFittype - Model to fit
        fittype object

    Examples
      openExample('curvefit/CreateFitTypesforLibraryModelsExample')
      openExample('curvefit/CreateCustomNonlinearModelsandSpecifyProblemParametersanExample')
      openExample('curvefit/CreateCustomLinearModelExample')
      openExample('curvefit/FitACurveDefinedByAFileExample')
      web /MATLAB/help/curvefit/fittype.html#btpaend-11
      web /MATLAB/help/curvefit/fittype.html#btpaend-14
      web /MATLAB/help/curvefit/fittype.html#btpaend-15

    See also fit, fitoptions, Curve Fitter

    Introduced in Curve Fitting Toolbox before R2006a
    Documentation for fittype
       doc fittype

help fit

 fit - Fit curve or surface to data
    This MATLAB function creates the fit to the data in x and y with the
    model specified by fitType.

    Syntax
      fitobject = fit(x,y,fitType)
      fitobject = fit([x,y],z,fitType)
      fitobject = fit(x,y,fitType,fitOptions)
      fitobject = fit(x,y,fitType,Name=Value)
      [fitobject,gof] = fit(x,y,fitType)
      [fitobject,gof,output] = fit(x,y,fitType)

    Input Arguments
      x - Data to fit
        matrix
      y - Data to fit
        vector
      z - Data to fit
        vector
      fitType - Model type to fit
        character vector | string scalar | string array |
        cell array of character vectors | anonymous function | fittype
      fitOptions - Algorithm options
        fitoptions

    Name-Value Arguments
      Options for All Fitting Methods
        Normalize - Option to center and scale data
          'off' (default) | 'on'
        Exclude - Points to exclude from fit
          expression | index vector | logical vector | empty
        Weights - Weights for fit
          [ ] (default) | vector
        problem - Values to assign to problem-dependent constants
          cell array | double

      Smoothing Options
        SmoothingParam - Smoothing parameter
          scalar value in the range (0,1)
        Span - Proportion of data points to use in local regressions
          0.25 (default) | scalar value in the range (0,1)

      Interpolation Options
        ExtrapolationMethod - Extrapolation method
          "auto" (default) | "none" | "linear" | "nearest" | "thinplate" |
          "biharmonic" | "pchip" | "cubic"

      Linear and Nonlinear Least-Squares Options
        Robust - Robust linear least-squares fitting method
          'off' (default) | LAR | Bisquare
        Lower - Lower bounds on coefficients to be fitted
          [ ] (default) | vector
        Upper - Upper bounds on coefficients to be fitted
          [ ] (default) | vector

      Nonlinear Least-Squares Options
        StartPoint - Initial values for the coefficients
          [ ] (default) | vector
        Algorithm - Algorithm to use for fitting procedure
          'Trust-Region' (default) | 'Levenberg-Marquardt'
        DiffMaxChange - Maximum change in coefficients for finite difference gradients
          0.1 (default)
        DiffMinChange - Minimum change in coefficients for finite difference gradients
          10^{–8} (default)
        Display - Display option in Command Window
          'notify' (default) | 'final' | 'iter' | 'off'
        MaxFunEvals - Maximum number of evaluations of model allowed
          600 (default)
        MaxIter - Maximum number of iterations allowed for fit
          400 (default)
        TolFun - Termination tolerance on model value
          10^{–6} (default)
        TolX - Termination tolerance on coefficient values
          10^{–6} (default)

    Output Arguments
      fitobject - Fit result
        cfit | sfit
      gof - Goodness-of-fit statistics
        gof structure
      output - Fitting algorithm information
        output structure

    Examples
      openExample('curvefit/FitaQuadraticCurveExample')
      openExample('curvefit/FitaPolynomialSurfaceExample')
      openExample('curvefit/FitaSurfaceUsingVariablesinaMATLABTableExample')
      openExample('curvefit/CreateFitOptionsandFitTypeBeforeFittingExample')
      openExample('curvefit/ComparePolynomialFitsforFitFunctionExample')
      openExample('curvefit/FitaCubicPolynomialSpecifyingNormalizeandRobustOptionsExample')
      openExample('curvefit/FitACurveDefinedByAFileExample')
      openExample('curvefit/ExcludePointsfromFitExample')
      openExample('curvefit/ExcludePointsandPlotFitShowingExcludedDataExample')
      openExample('curvefit/CompareExtrapolationMethodsExample')
      openExample('curvefit/FitaSmoothingSplineCurveandReturnGoodnessofFitInformatioExample')
      openExample('curvefit/FitaSingleTermExponentialExample')
      web /MATLAB/help/curvefit/fit.html#bto2vuv-11
      openExample('curvefit/FindCoefficientOrdertoSetStartPointsandBoundsExample')

    See also Curve Fitter, fittype, fitoptions, prepareCurveData,
      prepareSurfaceData, feval, plot, confint

    Introduced in Curve Fitting Toolbox before R2006a
    Documentation for fit
       doc fit

    Other uses of fit

       conbase/fit
       conswitch/fit
       contwostage/fit
       fitproblem/fit
       gmdistribution/fit
       incrementalClassificationECOC/fit
       incrementalClassificationKernel/fit
       incrementalClassificationNaiveBayes/fit
       incrementalDriftAwareLearner/fit
       incrementalOneClassSVM/fit
       incrementalPCA/fit
       incrementalRegressionLinear/fit
       incrementalRobustRandomCutForest/fit
       lime/fit
       linearDegradationModel/fit
       localmod/fit
       localmulti/fit
       mbcmodel.model/fit
       shapley/fit
       tireModel/fit
       xregmodel/fit

Note that the estimates of these coefficients will be very uncertain, due to the highly noisy data.

[Darya]

Could you please advise me where the error is. I am getting the coefficients, but they do not fit my data at all.

x=(data(:,1))'; 
y=(data(:,2))';
myfittype = fittype('F*exp(a*x-b*x.^2)', 'dependent',{'y'},'independent',{'x'}, 'coefficients',{'a','b','F'});
myfit = fit(x',y',myfittype)
plot(myfit,x,y)

[Walter Roberson]

Note that that is not the equation of a guassian. A guassian would involve exp(-(x-B)^2/(2*C^2)) . If you expand that then you get a coefficient for x^2 and a coefficient for x -- but you also get a constant coefficient. In order for that constant coefficient to vanish, B would have to be 0.. in which case there would be no a*x coefficient.

[John D'Errico]

@Walter Roberson

Well, actually, that IS an equation of a Gaussian. Just not in the standard form you expect to see. That is, start with this expression:

F*exp(a*x-b*x.^2)

Complete the square inside the parens. You will find that indeed, it does fall into a Gaussian form.

A constant term inside the parens can be absorbed into the coefficient outside the exponential, given the identity:

exp(-b*x^2 + a) == exp(a)*exp(-b*x^2)

So even though it does not look like a Gaussian model, it is one.

@Darya

I don't have your data. I don't see you attached it anywhere. And, while a picture may say a thousand words, it lacks a little when we start to talk about numbers.

Regardless, the problem is, you did not provide starting values for the parameters. If you don't do so, then fit uses random numbers. And those random numbers are NEVER going to be very good. You need to provide good starting values. At least something in the right ballpark.

If you do provide the data, I might still suggest you have a problem, because that curve just does not look like a gaussian bell curve. And it looks like you have seriously noisy data, so any result you get will be a poor predictor. If you want me to go into more depth, I would need to have your data.

[Darya]

This is my data.

data = xlsread('data.xlsx');
x=(data(:,1))'; 
y=(data(:,2))';
myfittype = fittype('F*exp(a*x-b*x.^2)', 'dependent',{'y'},'independent',{'x'}, 'coefficients',{'a','b','F'});
myfit = fit(x',y',myfittype)

Warning: Start point not provided, choosing random start point.

myfit = 
     General model:
     myfit(x) = F*exp(a*x-b*x.^2)
     Coefficients (with 95% confidence bounds):
       a =      0.2661
       b =      0.4182
       F =      0.4057

[Walter Roberson]

@John D'Errico

Oh, good point about the term being movable to the outside!

[John D'Errico]

data = xlsread('data.xlsx');
x=(data(:,1))'; 
y=(data(:,2))';

ALWAYS PLOT YOUR DATA! Actually, plot everything. (Well, you did plot it, so well done there.)

plot(x,y,'o')

I would note this does not show the classical bell shaped curve. Not even remotely close. It has a peak near x=26, but it is not at all symmetric. And that is a problem.

Regardless, I want you to see the warning you got.

"Warning: Start point not provided, choosing random start point."

And that should tell you to at least try to do better.

In that model, b will be related to the variance, but it will be roughly inversely proportional to that parameter. (With an extra factor of 2 thrown in.) If the standard deviation might be around 3 or so, then an estimate for b might be around 0.05. This is a parameter we will care about.

But what should a be? If the peak is around 26, then a classical Gaussian form would look like this:

exp(-((X-26)/sigma)^2/2)

If we expand that form with sigma around 3, then we would see this:

exp(-X^2/18 + (26*X)/9 - 338/9)

That suggests a should be roughly 3, thus 26/9. The problem is, that constant term will kill you, because it will look like exp(-338/9)

exp(-338/9)

ans = 4.8959e-17

And that means we will have serious variable scaling problems. However, we MIGHT just survive.

myfittype = fittype('F*exp(a*x-b*x.^2)', 'dependent',{'y'},'independent',{'x'}, 'coefficients',{'a','b','F'});
myfit = fit(x',y',myfittype,'start',[3 0.05 1e-17])

myfit = 
     General model:
     myfit(x) = F*exp(a*x-b*x.^2)
     Coefficients (with 95% confidence bounds):
       a =           3  (-0.3965, 6.396)
       b =     0.05399  (-0.006018, 0.114)
       F =   2.477e-18  (-1.163e-16, 1.213e-16)

plot(myfit,x,y)
grid on
ylim([0,4])

And we see a fit now. Looks like complete absolute crapola on a shingle. You should see it tries to go down to zero at the right.

The problem is, your model simply does not fit the data. A true Gaussian model will go to zero in the tails. And it will be symmetric around the peak. Your data does not show anything of the sort. If you allow me to change the model of course, we can do something a little better. I would add an offset in y, allowing the entire model to translate up and down in y.

myfittype2 = fittype('F*exp(a*x-b*x.^2) + c', 'dependent',{'y'},'independent',{'x'}, 'coefficients',{'a','b','F','c'});
myfit2 = fit(x',y',myfittype2,'start',[3 0.05 1e-17 2])

myfit2 = 
     General model:
     myfit2(x) = F*exp(a*x-b*x.^2) + c
     Coefficients (with 95% confidence bounds):
       a =           3  (-20.08, 26.08)
       b =     0.05777  (-0.3729, 0.4885)
       F =   1.132e-17  (-3.504e-15, 3.527e-15)
       c =           2  (0.06586, 3.934)

plot(myfit2,x,y)

Still complete crap. Though I do observe that the parameters c and a have not changed at all from their start points. That gets into the serious issues we have with parameter scaling.

Sorry. I'm not at all sure why you think anything of the sort should fit your data.

Answer 2

[Walter Roberson]

If you do not have the curve fitting toolbox, then you can try

x = %...
y = %...
objfun = @(abc) sum( (abc(1).*exp(abc(2).*x+abc(3).*x.^2) - y).^2 );
abc0 = [.1 .3 .4];
ABC = fminunc(objfun, abc0);
a = ABC(1); b = ABC(2); c = ABC(3);

Do not expect the fit to be all that good.

[Walter Roberson]

data = readmatrix('data.xlsx');
x = data(:,1);
y = data(:,2);
objfun = @(abc) sum( (abc(1).*exp(abc(2).*x+abc(3).*x.^2) - y).^2 );
abc0 = [.001 .003 .004];
objfun(abc0)

ans = 70.2212

ABC = fminunc(objfun, abc0);

Solver stopped prematurely.

fminunc stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+02.

a = ABC(1); b = ABC(2); c = ABC(3);
[a, b, c]

ans = 1×3

0.2550 0.2212 -0.0049

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

F = @(x) a.*exp(b.*x + c.*x.^2);
[xs, xorder] = sort(x);
ys = y(xorder);
Y = F(xs);
plot(xs, Y, '-', xs, ys, 'b.')

[Walter Roberson]

Fitting the first two points would require a guassian that rises sharply -- and consequently also falls sharply. It would fall to near the level of the first point, and stay low from there. The sum of squares of distance from that low value to the rest of the points would add up to a fair bit.

Whereas, the fit that is being generated here looks somewhat reasonable. The first two points are roughly the same distance from the line, so they balance out. The rest is similar to a linear fit in effect. You can see how this curve is a plausible fit.