Confidence bands after a fit
The function at the top of the procedure file, vfit(), implements the fitting function. It is used only for calculating estimated derivatives by numerical methods.
The function CPInterval, as you are getting it, uses anaDerivs() for the derivatives, a function that implements the derivative of the function using an analytic expression. The one included in the procedure file is for a particular fitting function (the one in vfit, of course).
You will need to replace either vfit or anaDerivs with one appropriate to your fitting function. This was written before we had FUNCREF in Igor, so it was difficult to write it to use an arbitrary fitting function.
To use this, call CPInterval with the value of X where you want the confidence band, plus the covariance matrix from the fit, the coefficient wave, desired confidence level expressed as a fraction (0.95 for 95 per cent) and the number of degrees of freedom (number of input data points minus number of fit coefficients). Returns the confidence interval.
You might use it this way:
Duplicate fitWave, upperConfWave, lowerConfWave
upperConfWave = fit_inputData + CPInterval(x, M_covar, W_coef, 0.95, V_npnts-numpnts(W_coef))
lowerConfWave = fit_inputData - CPInterval(x, M_covar, W_coef, 0.95, V_npnts-numpnts(W_coef))
Function/D vfit(cw, xx) Wave/D cw;Variable/D xx return ((cw*xx)/(xx+cw)) End Function CPInterval(xx, covar, params, conflevel, DegFree) Variable xx Wave covar Wave params Variable conflevel Variable DegFree Variable i Variable j Variable LpEnd = numpnts(params) Variable temp, Yvar Variable tP // Uncomment these three lines and the one below the Duplicate, in order to // use numerical estimates of derivatives. You might do that if your fitting function // is very complicated to compute, doesn't have a good analytic derivative, or if you lost // your installation of Mathematica... // Duplicate/O params, epsilon // epsilon = 1e-6 // epsilon = 1e-18 Duplicate/O params, dyda // calcDerivs(xx, params, dyda, epsilon) // You must comment out this line in order to use numerical derivatives. anaDerivs(xx, params, dyda) i = 0 YVar = 0 do temp = 0 j = 0 do temp += covar[j][i]*dyda[j] j += 1 while (j < LpEnd) YVar += temp*dyda[i] i += 1 while (i < LpEnd) tP = StudentT(confLevel, DegFree) return tP*sqrt(YVar) end Function calcDerivs(xx, params, dyda, epsilon) Variable xx Wave params Wave dyda Wave epsilon Duplicate/O params, theP Variable yhat = vfit(params, xx) Variable i = 0 Variable LpEnd = numpnts(params) do theP = params theP[i] = params[i]-epsilon[i] yhat = vfit(theP, xx) theP[i] = params[i]+epsilon[i] dyda[i] = (yhat - vfit(theP, xx))/(2*epsilon[i]) i += 1 while (i < LpEnd) end Function anaDerivs(xx, w, dyda) Variable xx Wave w Wave dyda Variable denom = xx+w dyda = xx/denom dyda = -w*xx/(denom*denom) end