curve fitting and constraints
masheroz
Mon, 03/16/2015 - 08:22 pm
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?
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'll have a look.
March 17, 2015 at 12:42 am - Permalink
Code to remove background from a 1 dimensional wave:
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:
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
I think the following should work (not tested):
March 18, 2015 at 12:58 am - Permalink
Oh, and also this should work (not tested):
March 18, 2015 at 01:06 am - Permalink