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