Fit Curve to a user defined function

Hello,

I have 2 defined variables, S and u, both with experimental data shown below:

S =

   0.403000000000000
   0.331000000000000
   0.270400000000000
   0.222200000000000
   0.181600000000000
   0.148900000000000
   0.122100000000000
   0.100500000000000
   0.080800000000000
   0.067000000000000
   0.054000000000000
   0.045200000000000
   0.036300000000000
   0.029500000000000
   0.024800000000000
   0.020100000000000
   0.016400000000000
   0.014200000000000
   0.011900000000000
   0.008800000000000
   0.008000000000000
   0.005100000000000
   0.005100000000000
u =

   0.957033515000000
   0.903278685000000
   0.815542061000000
   0.732249907000000
   0.656767397000000
   0.578471565000000
   0.506853196000000
   0.440901480000000
   0.379658753000000
   0.324856464000000
   0.276170883000000
   0.233324443000000
   0.196245474000000
   0.164407417000000
   0.137095397000000
   0.113974220000000
   0.094601760000000
   0.078298766000000
   0.064655532000000
   0.053307707000000
   0.043933989000000
   0.036147733000000
   0.029677689000000

I want to fit them to the Monod Equation using programmatic curve fitting, where the monod equation is given as u = (umax*S)/(Ks+S), where umax and Ks are constants.

However, I do not know how to define the Monod Equation as a function in MATLAB so that I can fit u and S to it, because I need to do so to obtain the 95% prediction intervals, instead I end up with the following error messages instead:

>>[fitresult,gof] = fit(S,u,'Monod');

Error using fittype>iCreateFromLibrary (line 414)
Library function Monod not found.

Error in fittype>iCreateFittype (line 345)
    obj = iCreateFromLibrary( obj, varargin{:} );

Error in fittype (line 330)
                obj = iCreateFittype( obj, varargin{:} );

Error in fit>iFit (line 165)
    model = fittype( fittypeobj, 'numindep', size( xdatain, 2 ) );

Error in fit (line 116)
[fitobj, goodness, output, convmsg] = iFit( xdatain, ydatain, fittypeobj, ...

and

>> [fitresult,gof] = fit(S,u,'u = (umax*S)/(Ks+S)')
Error using fittype>iDeduceCoefficients (line 621)
The independent variable x does not appear in the equation expression.
Use x in the expression or indicate another variable as the independent variable.

Error in fittype>iCreateCustomFittype (line 477)
    obj = iDeduceCoefficients(obj);

Error in fittype>iCreateFittype (line 353)
    obj = iCreateCustomFittype( obj, varargin{:} );

Error in fittype (line 330)
                obj = iCreateFittype( obj, varargin{:} );

Error in fit>iFit (line 165)
    model = fittype( fittypeobj, 'numindep', size( xdatain, 2 ) );

Error in fit (line 116)
[fitobj, goodness, output, convmsg] = iFit( xdatain, ydatain, fittypeobj, ...

How would you all solve this problem?

Answer 1

[John D'Errico]

You cannot assume that fit will know the name of every possible nonlinear model form. monod is apparently not one of the predefined types.

S = [0.403000000000000   0.331000000000000   0.270400000000000   0.222200000000000   0.181600000000000   0.148900000000000   0.122100000000000   0.100500000000000   0.080800000000000   0.067000000000000   0.054000000000000   0.045200000000000   0.036300000000000   0.029500000000000   0.024800000000000   0.020100000000000   0.016400000000000   0.014200000000000   0.011900000000000   0.008800000000000   0.008000000000000   0.005100000000000   0.005100000000000]';
u = [0.957033515000000 0.903278685000000   0.815542061000000    0.732249907000000 0.656767397000000   0.578471565000000   0.506853196000000   0.440901480000000   0.379658753000000   0.324856464000000   0.276170883000000   0.233324443000000   0.196245474000000   0.164407417000000  0.137095397000000   0.113974220000000   0.094601760000000   0.078298766000000   0.064655532000000   0.053307707000000   0.043933989000000   0.036147733000000   0.029677689000000]';
plot(S,u,'o')

mdl = fittype('(umax*S)/(Ks+S)','indep','S')

mdl = 
     General model:
     mdl(Ks,umax,S) = (umax*S)/(Ks+S)

fittedmdl = fit(S,u,mdl)

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

fittedmdl = 
     General model:
     fittedmdl(S) = (umax*S)/(Ks+S)
     Coefficients (with 95% confidence bounds):
       Ks =      0.2603  (0.2524, 0.2683)
       umax =       1.593  (1.566, 1.621)

plot(fittedmdl,S,u)

The fit seems quite reasonable overall, even where the legend covers up the last data point. (Sorry about that. Blame legend, and I am too lazy now to force it to be better.)

Fit would have been happier if I had provided initial guesses, but even though it complained, it was ok at the end. This was not a difficult model to fit, and a random start is adequate.

[Gavin Leong]

This works perfectly! How would you plot the 95% prediction bands in this case?

[John D'Errico]

Make sure you understand the somewhat subtle difference. But then try it. My guess is you were asking for confint.

help confint

--- help for cfit/confint ---

 CONFINT  Confidence intervals for the coefficients of a fit result object.
    CI = CONFINT(FITRESULT) returns 95% confidence intervals for the
    coefficient estimates.
 
    CI = CONFINT(FITRESULT,LEVEL) specifies the confidence level.
    LEVEL must be between 0 and 1, and has a default value of 0.95.

    Documentation for cfit/confint
       doc cfit/confint

    Other uses of confint

       sfit/confint

help predint

--- help for cfit/predint ---

 PREDINT  Prediction intervals for a fit result object or new observations.
    CI = PREDINT(FITRESULT,X,LEVEL) returns prediction intervals
    for a new Y value at the specified X value.  LEVEL is the
    confidence level and has a default value of 0.95.
 
    CI = PREDINT(FITRESULT,X,LEVEL,'INTOPT','SIMOPT') specifies the
    type of interval to compute.  'INTOPT' can be either 'observation'
    (the default) to compute bounds for a new observation, or
    'functional' to compute bounds for the curve evaluated at X.
    'SIMOPT' can be 'on' to compute simultaneous confidence bounds
    or 'off' to compute non-simultaneous bounds.
 
    If 'INTOPT' is 'functional', the bounds measure the uncertainty
    in estimating the curve.  If 'INTOPT' is 'observation', the bounds
    are wider to represent the addition uncertainty in predicting a
    new Y value (the curve plus random noise).
 
    Suppose the confidence level is 95% and 'INTOPT' is 'functional'.
    If 'SIMOPT' is 'off' (the default), then given a single pre-determined
    X value you have 95% confidence that the true curve lies between
    the confidence bounds.  If 'SIMOPT' is 'on', then you have 95%
    confidence that the entire curve (at all X values) lies between
    the bounds.

    Documentation for cfit/predint
       doc cfit/predint

    Other uses of predint

       sfit/predint

[Gavin Leong]

Oh I was asking for Prediction Intervals, not the Confidence Intervals. Basically I want to know how to determine and plot those such intervals on a graph (I believe it is represented by a dotted curve as shown when you use the curve fitting toolbox.)

[Gavin Leong]

Never mind, I have figured it out. Thank you so much for the help!

For reference I used the following code:

p11 = predint(fittedmdl,S,0.95,'observation','off')

plot(S,u, '.', 'markersize', 8), hold on, plot(fittedmdl), hold on, plot(S,p11,'m--'),  xlim([0 0.5]), ylim([0 1.2])
title('Nonsimultaneous Prediction Bands','FontSize',9) 
legend({'Data','Fitted curve', 'Prediction intervals'}, 'FontSize',8,'Location', 'southeast')