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
//  your installation of Mathematica...

//  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.
    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[1]
    dyda[0] = xx/denom
    dyda[1] = -w[0]*xx/(denom*denom)
end

Forum

Support

Gallery

Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More