Create a Code Snippet by Logging In or Creating an Account

Confidence bands after a fit

Fri, 03/14/2008 - 09:04 am

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

Confidence bands after a fit

Igor Pro 9

Igor XOP Toolkit

Igor NIDAQ Tools MX