How do I get the error of a fitting function?

Question

JJH 2019-1-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/439682-how-do-i-get-the-error-of-a-fitting-function

评论： Star Strider 2019-1-16

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2019-1-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/439682-how-do-i-get-the-error-of-a-fitting-function#answer_356329

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.

5 个评论
显示 3更早的评论隐藏 3更早的评论

JJH 2019-1-15

编辑：JJH 2019-1-15

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 2019-1-15

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 2019-1-16