Advice on use of Line fit statistic Pr and ChiSqr

If the time values are the independent variable, and they have errors, you should be using the /ODR=2 "errors in variables" method:

DisplayHelpTopic "Errors in Variables: Orthogonal Distance Regression"

In that case, the fit is inherently nonlinear even when fitting a line, and V_Pr is not generated. Instead, you get the usual coefficient sigmas and if you ask for it, a covariance matrix.

I have been through the help topic and examples given and not having used this aspect of curve fitting before, found them straight forward and easy to follow (so I have probably made mistakes:-) 

I applied ODR to two sets of XY pairs of the same event from two different instruments with the SD from WaveStat (V_sdev)

The first set is 19 pairs and the second set 14 pairs.

The first set fit converged properly whereas the second threw an error -"singular matrix of other numeric error"  When I arbitrarily removed the last points of the second set the fit converged.  The removed data looks quite acceptable.  Any clues?

In practice there are some 60 or more XY pairs with different lengths and slopes for analysis and curve fitting and an 'automated' approach is essential.

Here is the ouput with the reduced second set.

 I have also read the commentry from John Weeks in a thread from 01/10/2013  - "Orthogonal Distance Regression - weight with two SDs".  How would I go ablout determining if the slopes of the lines within each others error bands ?

I read the John's comment  "One complication of all this stuff is that the ODR fitting technique is inherently nonlinear, so normal statistics aren't exact"  and I am not sure how to move forward with analysis of the line fit....

I may be misunderstanding something here, but it looks as though you are filling the weighting waves with the standard deviation of the whole dataset, whereas I would have thought the weighting should be the (estimate of) the standard deviation for each data point, on a point by point basis?

If you had multiple independent measurements of each data point then you could use the standard deviation of these replicates to estimate the uncertainty of a given data point, and fill the weighting wave with these.

Have you tried adding these weighting (st dev) as error bars on a x-y plot so that you can see whether they look reasonable?

Hallo Kurt,

Thnaks for the relpy. You may well be right but I do not understand (and cant find it in the help files - but it may be there) how to generate the SD for weighting needed when starting with raw data of unknown errors. 

I can look at the data on a one off basis but I really need to operate in a batch mode and the convergence/non-convergence when I remove a data pair doesn't seem right

Hi Mike,

The problem here is your errors are unknown, but you need to somehow be able to estimate these 'unknown' errors because these are the weighting values that you need.

One approach is if you can a priori assume that a straight line is a good fit, then fit this line and calculate the residuals (this is an option in the built-in CurveFit. Then, further assuming the errors are, on average, uniform across your dataset, calculate the sd of the residuals and use this value.

Hope this helps!

1) I presume "DP" and "SP" are double precision and single precision. Basically, you should never use a single-precision coefficient wave. Igor copies the other waves into double-precision arrays, but the coefficient wave actually participates in the fit.

2) in ODR fitting, the X and Y weights do two things: as in non-ODR fitting, the sigmas are computed using those weights, which makes the chi-square into something that is useful. But in ODR fitting it also establishes a geometric scale factor. If the slope is extreme, like 10^3 or 10^-3, the perpendiculars from the solution to your data points are going to be either nearly vertical or nearly horizontal, and the fit loses any meaningful information in the horizontal or vertical direction. But the perpendicular distance is computed using the weighted data, not the raw data. That means that if your slopes aren't near 1, you really need some sort of weighting to establish the 2D geometrical scaling. The span of the X and Y values probably is a reasonable ad-hoc starting point for situations where you don't have a good estimate of actual errors. It's been a while since I looked at this in detail... With normal least-squares, you can assume that the (chi-square / N) should be 1, and then use (chi-square / N) as a factor to adjust sigma's. That's common when you need weighting (you know that there is a measurement error that is proportional to log(Y) for instance) but you only know it to within a constant multiplier. I would have to do some research to figure out the proper measure for ODR fitting, but something similar should be possible.

Thanks Kurt. I'm thinking about this!

1) I presume "DP" and "SP" are double precision and single precision. Basically, you should never use a single-precision coefficient wave. Igor copies the other waves into double-precision arrays, but the coefficient wave actually participates in the fit.

The prefix D1 and SP refer to the instruments that logged the data.  The X data are Igor Date/Time and of course in double precision. However, the Y data are in single precision. Thanks for pointing this out; I'll correct this

If the slope is extreme, like 10^3 or 10^-3, the perpendiculars from the solution to your data points are going to be either nearly vertical or nearly horizontal, and the fit loses any meaningful information in the horizontal or vertical direction.

For this event with D1, the Y data spans, from ~ 1300 to 40 Hz and the X data, from 3749504289.114 to 3749504289.308 (25/10/2022 00:58:09.114 to 25/10/2022 00:58:09.308). The slope is ~ 6400 Hz/s.  

How does one know that the correct weighting has been used?  Or perhaps more importantly the other way round, how does one know with a given weighting that the best straight line fit has been achieved.

Mike German wrote:

1) I presume "DP" and "SP" are double precision and single precision. Basically, you should never use a single-precision coefficient wave. Igor copies the other waves into double-precision arrays, but the coefficient wave actually participates in the fit.

The prefix D1 and SP refer to the instruments that logged the data.  The X data are Igor Date/Time and of course in double precision. However, the Y data are in single precision. Thanks for pointing this out; I'll correct this

Oops. I read "SP" and then my eyes automatically registered "DP" instead of "D1" or "D2". I must be a programmer :)

Quote:

If the slope is extreme, like 10^3 or 10^-3, the perpendiculars from the solution to your data points are going to be either nearly vertical or nearly horizontal, and the fit loses any meaningful information in the horizontal or vertical direction.

For this event with D1, the Y data spans, from ~ 1300 to 40 Hz and the X data, from 3749504289.114 to 3749504289.308 (25/10/2022 00:58:09.114 to 25/10/2022 00:58:09.308). The slope is ~ 6400 Hz/s.  

How does one know that the correct weighting has been used?  Or perhaps more importantly the other way round, how does one know with a given weighting that the best straight line fit has been achieved.

Without an independent way to estimate actual errors, you can't really tell if you have "correct" weighting. In some cases, it is reasonable to assume that the errors are a proportion of the value span. You will get the same fit with the X and Y weights multiplied by the same factor, but the estimated coefficient errors will scale with that multiplier. For an ordinary least squares fit, if you have the correct weights, the chi-square as Igor reports it should be about N (number of data points). I don't know how that compares to ODR fitting, or how to combine the two sets of weights, but I note that if you multiply the weights by 10, the solution is the same, and the sigmas go up by a factor of 10, and chi-square by a factor of 100.

I would say that it would be reasonable to set the X weights to 1300 - 40 = 1260 and the Y weights to 3749504289.308 - 3749504289.114 = 0.194. Be sure to use the /I=1 flag- that's Igor's historical oddity that says that the weights are SD, not 1/SD. You can get outputs of both X and Y residuals from an ODR fit, that might give you some indication of the actual distribution of errors.

OK I am now using double precision for both X and Y data and the spans of X and Y for the weighting waves. D2_Fdp and SP_Fdp are the two Y waves.  From history area:-

As mentioned elsewhere the sigmas are clearly not right. Using the above weighting it was necessary to remove an outlier on the SP_Fdp data (point 6) on the previous "XY data with non-ODR fit" atatchment to get convergence!  I am not sure how I will overcome this editing without eye-balling which isn't practical for batch processing.  The fits and residuals are at the attachment

For comparison the same data using non-ODR curvefit and no weighting. The error estiamtes (sigmas) are sensible

I find it difficult to see the benefit of using ODR other than it is teh right way to do it.:-)

Can you post your Igor experiment file? Or send it to support@wavemetrics.com?

The large sigmas are probably because using the data spans as weight is a considerable over-estimation of measurement errors. It *probably* gives you X and Y weights that are nearly proportional.

Mike German wrote:

Ideally I would like a number that told me how much I can depend on the paramater based on the line fit. 

Hi Mike,

Here is a thought - not sure about the robustness of the stats though - can you fit you data to both a least squares on y and a least squares on x, and then use the difference in the fitted slopes to determine some measure of the uncertainty in the (fitted) slope of the data (which could be estimated as the average of these two slopes)?

In case it helps to visualize this, below is some code I wrote many years ago to generate some example plots when writing a curve fitting tutorial for my work colleagues. In the generated graph, green is least-squares y, blue is least-squares x, and red is ODR.

Nice demo, KurtB!

Please note that the demo has data with a slope that is nearly 1, which avoids all this bother with the weighting.

There is a demonstration of the importance of the weighting in our documentation: DisplayHelpTopic "ODR Fitting Examples" I used error bars to represent the XY residuals in those examples. I like KurtB's use of lines to the fit line much better; perhaps I can revise my documentation.

I draw your attention to the examples of exponential fits with and without weighting. I think that makes a pretty good visualization of the importance of weighting.

As I mentioned previously, at this point it may be more useful for me to have your actual data that I can work with myself.

Yes John, agreed it is a nice demo from Kurt.  Thanks for sharing Kurt it does nicely show the construction.

I have sent a sample of my real data to John as requested.

Not withstanding the recommendations on how to do robust regression fitting, is your system amenable to an analytical approach to determine an uncertainty budget? The method uses (linear) propagation of uncertainty.

https://en.wikipedia.org/wiki/Propagation_of_uncertainty

By example, when your parameter P depends on slope m as a function P(m), I would imagine doing a derivation to obtain the uncertainty DP as a function of the uncertainty Dm in a form (DP/P)^2 = f((Dm/m)^2).

I am going to look that up as I am not familiar with it. Thanks