curve fitting and constraints

Share icon

I have a several thousand datasets from which I want to subtract a quadratic background. (single example attached)

I think that the best way to do so is to fit a curve subject to an additional constraint that the value of the curve must be less than or equal to the data at every point.

The constraints make it easy to constrain the value of the coefficients, but I can't see how to further constrain based on the actual data values.

In Excel, I could add a column for the result of (data - calculated background) and have a single cell equal to the minimum of that column, and constrain the model such that the cell value is greater than or equal to zero.

.

Question: How can I constrain the model such that the model is always less than the data? curveFit.pxp (43.63 KB)
You could use an approach similar to the one I posted here:
http://www.igorexchange.com/node/6006

In outline the process is as follows:
1. Fit your quadratic (or whatever function you desire) and calculate residuals.
2. Determine the data points that have the largest (positive) values for their residuals and mask these from subsequent fits.
3. Repeat 1 and 2 a few times.
The parameters are the number/fraction of points you exclude at each iteration, and the number of iterations.

The fit will then end up being on the background only (if the data is well behaved).
Note: there is a chance that a quadratic fit actually fits the peak initially - this is more likely to happen if there is more peak than background in your data. In this case, I suggest you make the first iteration a straight line fit.

HTH,
Kurt

March 17, 2015 at 12:21 am - Permalink

I made it work, just not very elegantly. (apologies for making your eyes bleed)

Code to remove background from a 1 dimensional wave:
//in-place removal of BG.
function removeBG(withBG)
    wave withBG
 
    Make/D/O/N=(3) W_coef  // the same name as for CurveFit
    W_coef[0]=0.5 // offset
    W_coef[1]=0 // scale
    W_coef[2]=0 // scale
    
    Duplicate/O withBG, fitme, waveFit, diff  //these are the waves I need to fit the function
    
    variable i,j
    
    //this loop goes though and does a quadratic curve fit to the data
    //  after this, it finds the point that is the largest positive distance away
    //  from the curve and removes it from the calculation. It does this until 
    //  there are 4 points left, being one more than the degree of the curve.
    for(i = 0; i < numpnts(withBG) - 4; i+=1)
        CurveFit/Q poly 3, fitMe /D //fit the wave 
        for(j = 0; j < numpnts(waveFit); j += 1) //calc the value of the fit at each row value. there is probably an "igor way" of doing it, but I don't know it
            waveFit[j] = W_coef[0] + W_coef[1]*j + W_coef[2]*j^2
        endfor
        diff = fitme - waveFit //calc the difference plot
        WaveStats/Q diff //find the location of the max values
        fitme[V_maxRowLoc] = NaN //remove the points with the biggest positive difference
        diff[V_maxRowLoc] = NaN 
    endfor
    
    //One more curve fit with the last few points
    CurveFit/Q poly 3, fitMe /D 
    for(j = 0; j < numpnts(waveFit); j += 1)
        waveFit[j] = W_coef[0] + W_coef[1]*j + W_coef[2]*j^2
    endfor
 
    withBG = withBG - waveFit //withBG now has no BG
    
    //I don't need these, so I might as well clean them up
    KillWaves/Z fitme, waveFit, diff, W_coef, fit_fitme, W_sigma, W_ParamConfidenceInterval
end


code to apply the 1D case to a 2D wave of 1D waves that I need to fit:
function removeBGsurface(surf)
    wave surf
    variable ii,jj
    
    Make/O/N=(DimSize(surf,1)) wi
    
    for(ii = 0; ii < DimSize(surf, 0); ii += 1)
        print ii
        wi = surf[ii][p]
        removeBG(wi)
        
        //I can't figure out how to do the assignment in igor-speak
        for(jj=0; jj < DimSize(surf,1); jj += 1)
            surf[ii][jj] = wi[jj]
        endfor
    endfor
end

March 17, 2015 at 10:13 pm - Permalink

masheroz wrote: 

        //I can't figure out how to do the assignment in igor-speak
        for(jj=0; jj < DimSize(surf,1); jj += 1)
            surf[ii][jj] = wi[jj]
        endfor


I think the following should work (not tested):
surf[ii][] = wi[q]

March 18, 2015 at 12:58 am - Permalink

masheroz wrote: 

        for(j = 0; j < numpnts(waveFit); j += 1) //calc the value of the fit at each row value. there is probably an "igor way" of doing it, but I don't know it
            waveFit[j] = W_coef[0] + W_coef[1]*j + W_coef[2]*j^2
        endfor


Oh, and also this should work (not tested):
waveFit[] = W_coef[0] + W_coef[1]*p + W_coef[2]*p^2

March 18, 2015 at 01:06 am - Permalink