fitiing with integrals

Share icon

Fri, 11/17/2023 - 09:12 am

hello, i need to fit a data set to the integral of a gaussian (which looks like a sigmoid) but for biological reasons, the fitting parameters are those of the gaussian. I don't know how to use either integrate1D or integrate in a FuncFit (user defined fitting function). Does anyone have an example of doing such things?

thanks a million

diego

 

I think you can simply use the function StatsNormalCDF instead of doing the integration within the fit function.

e.g.:

Function MakeData()
    Make/O/D/N=50 wData
    SetScale x -5,5,"", wData
    //wData = StatsNormalCDF(x,0,1) + gnoise(00.05)
    wData = StatsNormalPDF(x,0,1) + gnoise(0.1)
    Integrate wData/D=W_INT
    Display W_INT
    AppendToGraph/R wData
    ModifyGraph mode=3,rgb(W_INT)=(0,0,65280)
    ModifyGraph rgb(wData)=(0,52224,0) 
End


Function Fit_NormalCDF(w,x) : FitFunc
    Wave w
    Variable x

    //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
    //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
    //CurveFitDialog/ Equation:
    //CurveFitDialog/ f(x) = StatsNormalCDF(x,centre,width)
    //CurveFitDialog/ End of Equation
    //CurveFitDialog/ Independent Variables 1
    //CurveFitDialog/ x
    //CurveFitDialog/ Coefficients 4
    //CurveFitDialog/ w[0] = centre
    //CurveFitDialog/ w[1] = width
    //CurveFitDialog/ w[2] = height
    //CurveFitDialog/ w[3] = offset

    return w[3]+w[2]*StatsNormalCDF(x,w[0],w[1])
End

Make/D/N=4/O W_coef
W_coef[0] = {0,1,1,0}
FuncFit/NTHR=0/TBOX=768 Fit_NormalCDF W_coef  W_INT /D

EDIT - I have updated the code to explicitly make the test data using integrate, and added height and offset to the fit.

November 17, 2023 at 09:37 am - Permalink

But since there is a well-known analytic integral for a Gaussian, use that instead. In fact, that will be available in either the StatsNormalCDF() function, or the erfc() function. Naturally, you will need to figure out how to transform the result to account for parameters you use in the Gaussian function.

November 20, 2023 at 11:03 am - Permalink

Here is an example of two functions I have used for fitting XPS peaks as a Gaussian peak with a Shirley background, which is the integral of the peak, calculated using erf().

ThreadSafe Static Function XPSFitAssGauss(x, Position, FWHM)
//  Returns a Gaussian peak normalized to an area of one
Variable x, Position, FWHM

    //  f(x) = sqrt( 4 * ln(2) / pi ) / FWHM * exp( -4 * ln(2) * (( x - Position ) / FWHM )^2 )
    //  sqrt( 4 ln(2) / pi ) = 0.939437
    //  4*ln(2) = 2.77259

    Return 0.939437/FWHM*exp(-2.77259*((x-Position)/FWHM)^2)
end



ThreadSafe Static Function XPSFitAssGaussShirley(x, Position, FWHM)
//  Calculates the Shirley background (integrated area) of a Gaussian peak with an area of one
//  D.A. Shirley, Phys. Rev. 55 (1972) 4709
Variable x, Position, FWHM

    //  f(x) = 0.5 * ( Erf(sqrt( 4 * ln(2) ) * ( x - Position ) / FWHM ) + 1 )
    //  sqrt( 4 ln(2) ) = 1.66511

    Return 0.5*(Erf(1.66511*(x-Position)/FWHM)+1)
end

 

November 20, 2023 at 12:34 pm - Permalink

i am very thankful to your answers, the fit_normalCDF does work but i still think i need to integrate a gaussian. i apologize for my poor coding skills, is there a flagrant error in this code? its trying to fit all at once. again i apologize if what i am writing is absurd.

function gaussintall(pw,yw,xw) : FitFunc
    WAVE pw,yw,xw

    variable/g a=pw[0],b=pw[1],c=pw[2],d=pw[3]

    variable i
    for(i=0;i<(numpnts(yw));i+=1)
        yw[i] = integrate1D(gaussE,0,i,1)// integrates the gaussian one point at a time    
    endfor
end

function gausse(inx)
    variable inx
   
    NVAR a,b,c,d
    return a + b*exp(-((inx-c)/d)^2)
end

 

November 23, 2023 at 10:30 am - Permalink