# 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[0]*xx)/(xx+cw[1]))

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

//  Duplicate/O params, epsilon
//  epsilon[0] = 1e-6
//  epsilon[1] = 1e-18

Duplicate/O params, dyda

//  calcDerivs(xx, params, dyda, epsilon)

//  You must comment out this line in order to use numerical derivatives.
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

Variable xx
Wave w
Wave dyda

Variable denom = xx+w[1]
dyda[0] = xx/denom
dyda[1] = -w[0]*xx/(denom*denom)
end

Forum

Support

Gallery