Igor Code
The objective of this page is to serve a repository of useful Igor code that we keep needing for different applications. Feel free to add new code or to comment on existing code, but please don't delete old code w/o discussing with Jose.
Doing a Robust Linear Regression
Linear regression gets pretty complex when one starts thinking about the details. There are 3 important decisions, which depend on the data which you are trying to fit:
- (A) Error in one or both coordinates:
- (A1) Only one of the coordinates (Y) has error and the other does not (e.g. when fitting an aerosol measurement vs. time).
- (A2) Or do both coordinates have errors of the same relative order-of-magnitude? (e.g. as when comparing two instruments that measure the same quantity)
- (B) Assumptions about the relative magnitudes of the errors in different points
- (B1) All the points are assumed to have the same error
- (B2) The error varies point-by-point.
- (C) Standard vs. Robust regression
- (C1) The fit line is found by minimizing Chi^2, i.e. the sum of the squares of the deviations between each point and the regression line.
- (C2) The fit line is found by minimizing the sum of the absolute values of the deviations between each point and the regression line.
This gives rise to 2 x 2 x 2 = 8 options to choose from when fitting your data. Which to use depends on the nature of the data and the use of the fit, and should be chosen intelligently. In practice if there is a good relationship and the points are scattered around a line with few outliers, all the methods will produce similar results. However when there are outliers or significant scatter, there can be significant differences and some of the methods are plainly incorrect for some datasets.
To do regressions in Igor in each case:
- A1-B1-C1: this is the default assumption of the linear fit
- A2-B1-C1: an approximation to this fit can be done with the /ODR=2 option in Igor.
- A1-B1-C2: this can be done in Igor, but you need to define a new custom fit function (which can be just a line) and then set V_FitOptions = 2 in the command line. (Going back to V_FitOptions = 0 gets you back to the regular A1-B1-C1 fits)
- A2-B1-C2: this can be accomplished approximately by doing A1-B1-C2 for Y vs X, then for X vs Y, and averaging the fitting parameters appropriately (and the fitting parameter uncertainties in quadrature). There are some papers that show that this is approximately correct.
- For all the B2 options, choose the fitting weights in the Igor Curve Fitting tab
Regridding a Wave into the Time Axis of Another Wave
To regrid a wave into another framework, use the following functions. For example you might have one curve, defined by a certain x- and y-wave (xwave1, ywave1). Let's say there are 10 points in these waves, and the maximum is 20 (and extreme case of length discrepancy). Another curve is defined by two other x- and y-waves, also with a maximum of 20, but with 20 points each (xwave2, ywave2). It is important that they have approximately the same bound values, or it will have a hard time interpolating between them.
If I want to look at the first curve in terms of the second x-wave, you have to regrid the first y-wave into the x-wave of the second. This can be done in two ways, and that is dependent on which inintial set of waves has more data points. If the wave you want to regrid into has more points than the original, then use the first of these functions (Interp_1). If the opposite is true (wave you want to regrid to has less points) use the second function (Interp_2). They can be run as is, after changing the wave names appropriately. The [p] refers to each point, and should stay with this syntax (replacing it with a number refers only to that specified cell).
Function Interp_1 ()
- Wave xwave1, ywave1, xwave2 //ywave2 is not used
- duplicate/o xwave2 ywave_regrid
- ywave_regrid = ywave1 [binarysearchinterp(xwave1, xwave2[p])]
- appendtotable ywave_regrid
End Function
Function Interp_2 ()
- Wave xwave1, ywave2, xwave2
- duplicate/o xwave1 ywave_regrid
- ywave_regrid = ywave1 [fAverageXY(xwave2, xwave1, xwave2[p-1], xwave2[p])]
End Function
Compiled by Qi Zhang and Alex Huffman, 2004.
