How do I get the error of a fitting function?

I've written a function that does a Lorentzian fit on a set of data. I want to calculate an error bar for the accuray of this fit with respect to the x coordinates (the parameter Wavelength). Is this a built in command or is it more difficult?

function [params, Wavelength, Intensity, wavpeak]=SpectrumPeakFitCSV(filename)
    
    SpectrumData=csvread(filename);
    Wavelength=SpectrumData(2:end-1,1);
    Intensity=SpectrumData(2:end-1,2);
    Max = max(Intensity); 
    peak = find(Intensity==Max(1));
    wavpeak = Wavelength(peak(1));
    fun=@(x,Wav)Lorentzian(x(1),x(2),x(3),Wav)+x(4);
    dark=Intensity(1);
    opts = optimset('Display','off');
    x=lsqcurvefit(fun,[Max(1),1,wavpeak,dark],Wavelength,Intensity,[0.5*Max(1) 0 wavpeak-0.2 0.9*dark],[2*Max(1) 3 wavpeak+0.2 1.5*dark],opts);
    params=x; %scale, width, centre, dark
end

Answer 1

[Star Strider]

If I understand correctly what you want to do, the Statistics and Machine Learning Toolbox function nlpredci (link) is likely what you are looking for. Use the ‘jacobian’ output (last output) from lsqcurvefit with it.

[JJH]

I've been looking into this more and I'm not sure exactly what you mean by use the Jacobian output from lsqcurvefit with nlpredci. I can get lsqcurvefit to print a Jacobian but I'm confused about how to use this in conjunction with nlpredci. Could you explain a bit more?

[JJH]

To clarify, I've modified my code to look like this so far

function [params, Wavelength, Intensity, wavpeak]=SpectrumPeakFitCSV(filename)
    
    SpectrumData=csvread(filename);
    Wavelength=SpectrumData(2:end-1,1);
    Intensity=SpectrumData(2:end-1,2);
    Max = max(Intensity); 
    peak = find(Intensity==Max(1));
    wavpeak = Wavelength(peak(1));
    fun=@(x,Wav)Lorentzian(x(1),x(2),x(3),Wav)+x(4);
    dark=Intensity(1);
    opts = optimset('Display','off');
    [x,~,R,~,~,~,J]=lsqcurvefit(fun,[Max(1),1,wavpeak,dark],Wavelength,Intensity,[0.5*Max(1) 0 wavpeak-0.2 0.9*dark],[2*Max(1) 3 wavpeak+0.2 1.5*dark],opts)
    params=x; %scale, width, centre, dark
    [Ypred,delta] = nlpredci(fun,Wavelength,x,R,'Jacobian',J)
end

but I get a few errors doing it this way:

Error using svd
SVD does not support sparse matrices. Use SVDS to compute a subset of the singular values and vectors of a sparse matrix.

Error in internal.stats.isEstimable (line 108)
            [V,Sx,U] = svd(X,0);

Error in nlpredci (line 306)
        [estimable,rankJ,pinvJTJ] = internal.stats.isEstimable(delta, 'DesignMatrix', J,'TolSVD',TolSVD);

Error in SpectrumPeakFitCSV (line 19)
    [Ypred,delta] = nlpredci(fun,Wavelength,x,R,'Jacobian',J)

Error in VCSELCharacLIVPlusJunctionTempLowRes (line 111)
        peakfit = SpectrumPeakFitCSV([FolderName '\' Direc(kk).name]);

[Star Strider]

Your nlpredci call uses the 'Jacobian', as I intended in my Answer (I’ve used it similarly in the past with lsqcurvefit).

With respect to the error nlpredci is throwing, it seems that something (possibly the Jacobian) is sparse (i.e. has a significant number of zero elements). I don’t know what to suggest, since using the svds function would involve hacking into the nlpredci code (if that’s even possible) to change that one call.

I would begin by looking at the ‘J’ output, and possibly plotting it (surf or contourf) to see what it looks like if it is too large to examine easily when you display it. (You can also see if using svds on your Jacobian produces acceptable results. This will not change anything with respect to nlpredci, although it could provide you with some information about the Jacobian matrix.)

Alternatively, see if using the nlinfit function and the covariance matrix with nlpredci works. This also assumes that the covariance matrix will have a relatively low condition number (not close to being singular). You appear to be regressing two vectors (as opposed to the dependent variable being a matrix), so nlinfit should work with your data.

[Janna Hinchliff]

Thanks, nlinfit seems to be a bit easier!

[Star Strider]

@Janna Hinchliff — Please see my Comment to your Question.