[rfc] Inexact quantities arithmetics
Hello to everyone, I'm developing a library for working with inexact quantities, like (0.3162 +- 0.0079), in other words, for auto-estimating the propagation of errors in calculations. I plan to release it under Boost license when ready. Any comments (including sensibility of the whole project) are welcome. I. Some basic principles and limitations I established so far 1. Inexact quantity is a pair of 'value' and 'error' components. 2. Error component is always non-negative, but can be infinite or NaN. 3. Value and error are of the same type (at least for the user). For example, value can not be 'double' while error is 'float', nor value can be 'complex<double>' while error is 'double'. 4. Value result of calculation does not depend on error and coincides with result of the same calculation for bare value. Correspondingly it is biased from math expectation for any non-linear function. For example, cos(0 +- 0.1).value() gives 1, not exp(-0.05). 5. Error calculation is limited to the first order. This means f(mx +- dx) gives (f(mx) +- abs(dx * f'(mx))), f'' and subsequent terms proportional to greater powers of dx are ignored. For example: cos(0 +- 0.1) gives exactly 1; sqrt(0 +- 0.1) gives (0 +- inf). Later, template parameter for altering behavior described in the previous two items can be added. 6. Correlations are not considered; instead, worst case is always assumed. For example: (mx +- dx) + (my +- dy) gives ((mx + my) +- (dx + dy), not ((mx + my) +- sqrt(dx*dx + dy*dy); x *= x (where x is of inexact type) produces correct result, while x -= x will not produce 0 but, incorrectly, (0 +- 2*x.error()). In the distant future additional operations with correlations (and distributions) explicitly specified by user can probably be added. NOTE: Items 4-6 make specific distribution (normal, lognormal or whatever) irrelevant. II. Proposed IO format 1. Currently IO format is specified by template parameter. I really need your comments here. I want IO format to be both tweakable and replaceable by user with minimum headache. Ideally I'd like to see my specific manipulators similar to std::scientific and std::precision(int), but I'm not sure it's possible. 2. Already, the following format for inexact<T> is implemented: a) When T is simple floating-point type, the format is: '(' value dummy error ')' [('e' | 'E') common_exponent] Here dummy is any whitespace-terminated string. For wchar_t streams sensible value is L" \xB1 ". For char streams default value is " +- ", but one might decide to use " \pm ", or " ± ", or UTF-8 representation of L" \xB1 ", or probably something else as long as it ends on space. The point is to read parenthesized part with: in >> value >> ws >> dummy >> ws >> error; In output, presence of common exponent and specific format of value and error (fixed or scientific, precision) depend on stream flags and values of value and error in a complex way I will not elaborate on here. Results look like this: Default: (-0.0010676 +- 0.0000010) (-1.16767e+06 +- 0.0011) Fixed: (-0.001068 +- 0.000001) (-1167672.1657 +- 0.0011) Scientific: (-1.0676 +- 0.0010)e-3 (-1.167672e+06 +- 1.1e-03) b) Alternatively, just the bare value (without parentheses) can be used in place of (value +- 0). c) When T is collection C<E> of elements E, the format is by definition the same as for C<inexact<E> >. Corresponding converters are provided. For example, output of both inexact<complex<double> >( complex<double>(1, 2), complex<double>(3, 4)) and complex<inexact<double> >( inexact<double>(1, 3), inexact<double>(2, 4)) is "((1 +- 2),(3 +- 4))". 3. Format where error is stored in and retrieved from the number of value digits can be useful as well. III. Some notes about implementation 1. In implementation the type of error is wrapped in order to guarantee non-negativeness. I'm working on this right now. a) When value type T is simple floating-point type, error will be of type absolute<T>, that trivially converts to T but calls abs() when converting from T. b) When T is collection C<E> of elements E, error will be of type C<A>, where A is error type of E. For example, in inexact<vector<complex<double> > > the value type will be vector<complex<double> >, and error type will be vector<complex<absolute<double> > >. Hopefully all this stuff will get optimized out. I tried to do without this wrapping first, but code is much more transparent with. 2. I'll need derivatives for all cmath functions at least. It is probable that some template-based derivation-finder will be created as a by-product. Or, does one exist already? With Best Regards, Marat P.S. Pardon for possible double-post. I'm really lost how boost list + gmane combination works.
On Sun, Jul 08, 2007 at 10:07:29PM +0400, Marat Khalili <0x8.0p15@rol.ru> wrote:
Hello to everyone,
I'm developing a library for working with inexact quantities, like (0.3162 +- 0.0079), in other words, for auto-estimating the propagation of errors in calculations. I plan to release it under Boost license when ready. Any comments (including sensibility of the whole project) are welcome.
This looks like interval arithmetic: http://www.boost.org/libs/numeric/interval/doc/interval.htm and: http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2006/n2067.pdf Regards Andreas Pokorny
Andreas Pokorny wrote:
This looks like interval arithmetic:
http://www.boost.org/libs/numeric/interval/doc/interval.htm and: http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2006/n2067.pdf
Yes, this point needs to be explained. I did know this library exists, but it was not directly useful, there were differences: 1) First of all if one uses intervals instead if inexact quantities result for value (what's it? middle of the interval?) is affected. One might not want to estimate error this way. 2) Then, it only supports 'totally ordered' base types, no complex, vectors and so on. Adding support for these is most part of the work. 3) And it doesn't use derivatives, but simply calculates each function twice. It can be faster or slower, but the result is different. 4) Moreover, dyadic functions, mostly pow(inexact, inexact), are not implemented there yet. Implementing them up to the second order is the most part of the fun. :) 5) On the other hand, there are some interval class features, like rounding, I really don't need. What is true is interval can be treated as a distribution (uniformly between specified values), so results have to be similar sometimes, especially at the first order. So, Could it be wrapper for the interval class? Probably it could, but it looked like headache for me. Will inexact class support converting to and from intervals? Most likely, but to confidence intervals with all this quantile machinery. With Best Regards, Marat
Marat Khalili wrote:
I. Some basic principles and limitations I established so far
<snip>
6. Correlations are not considered; instead, worst case is always assumed. For example: (mx +- dx) + (my +- dy) gives ((mx + my) +- (dx + dy), not ((mx + my) +- sqrt(dx*dx + dy*dy); x *= x (where x is of inexact type) produces correct result, while x -= x will not produce 0 but, incorrectly, (0 +- 2*x.error()).
In the distant future additional operations with correlations (and distributions) explicitly specified by user can probably be added.
Is there a use case where this behavior is useful and meaningful? Without this feature, I can't come up with any reason for me to use it. I certainly can't see a use for it in publication quality basic research or any engineering work where you need to rely on the errors it calculates. Do you have an intended use case? I don't want to sound too critical of your work (although it is hard to NOT sound critical, I suppose ... sorry). I know how hard this area is to do right. I've tried to build my own class to handle this type of work and failed miserably. I've come to the conclusion that it isn't even possible to do in pure C++ without substantial run time performance penalties (although I would love to be proved wrong). The problem as I see it is that proper error analysis requires a level of global knowledge to which the language just doesn't give you access at compile time. Now, maybe you could do it with sufficient compile-time reflection/introspection support that interacts well with current template metaprogramming facilities, but even then I'm not sure it can be done. Perhaps you could write a compile-time DSL that does the right thing, but then you aren't programming transparently in C++ anymore. There was a brief exchange concerning this topic during the review of the Quantitative Units Library back in the first half of the year. Matthias Schabel was persuaded not to include his version of this class in the library at the time since it didn't "do the right thing". All that said, I would love to see an efficient library solution to this problem ... I'd love to be able to use real C++ and not the current poorly performing, hacky solutions that I have access to. -- ------------------------------------------------------------------------------- Kevin Lynch voice: (617) 353-6025 Physics Department Fax: (617) 353-9393 Boston University office: PRB-361 590 Commonwealth Ave. e-mail: krlynch@bu.edu Boston, MA 02215 USA http://budoe.bu.edu/~krlynch -------------------------------------------------------------------------------
Hello, At least someone else needs this feature if not implementation, that's great. Kevin Lynch wrote:
6. Correlations are not considered; instead, worst case is always assumed. Is there a use case where this behavior is useful and meaningful?
I do have this kind of case, but I don't know if it's typical or not. It depends on: (a) how many independent sources of error do you have, and (b) is it fine to have just upper bound for error. In my case I have few numerical integrations (each producing inexact value) and various operations on the results. I just wanted to be sure errors do not grow too large. Every time I add uncorrelated values as correlated I'm wrong by factor of two at most, and who can tell for sure if my integrations are really uncorrelated or not.
Without this feature, I can't come up with any reason for me to use it.
You can always write: z = inexact( x.value() + y.value(), sqrt(x.error()*x.error() + y.error()*y.error())) in place of z = x + y where I really need it. Feature is just formalizing it. Problem is, one might want to add x and z later. Zeroth principle for me was, it is not opaque library, but just some help with arithmetics. It will not allow person knowing nothing about errors to calculate errors.
I certainly can't see a use for it in publication quality basic research or any engineering work where you need to rely on the errors it calculates.
It correctly produces upper bounds almost everywhere. Where it doesn't (like cos(0.01 +- 0.1)) unexperienced physicist will also make a mistake :).
Do you have an intended use case? I don't want to sound too critical of your work (although it is hard to NOT sound critical, I suppose ... sorry).
No, I did need it myself, another question is I already spent more (like x10 :)) time on the library than just doing it manually could take. So I hope it will be of a use for someone else, that's why I need some criticism, really.
I know how hard this area is to do right. I've tried to build my own class to handle this type of work and failed miserably.
I'd be happy to hear more about your experience.
I've come to the conclusion that it isn't even possible to do in pure C++ without substantial run time performance penalties (although I would love to be proved wrong).
Yes, at the very least we would have to keep lots of correlations in memory. N^2 correlations for N floating-point variables in the worse case.
The problem as I see it is that proper error analysis requires a level of global knowledge to which the language just doesn't give you access at compile time. Now, maybe you could do it with sufficient compile-time reflection/introspection support that interacts well with current template metaprogramming facilities, but even then I'm not sure it can be done. Perhaps you could write a compile-time DSL that does the right thing, but then you aren't programming transparently in C++ anymore.
Well, it's an interesting idea - this can be done in Java with reflection and runtime disassembling (Soot anyone here?). If only I had some time... :(
There was a brief exchange concerning this topic during the review of the Quantitative Units Library back in the first half of the year. Matthias Schabel was persuaded not to include his version of this class in the library at the time since it didn't "do the right thing".
Was it in this list? Where can I take a look at the library?
All that said, I would love to see an efficient library solution to this problem ... I'd love to be able to use real C++ and not the current poorly performing, hacky solutions that I have access to.
One day temptation to write a helper function or two becomes too strong :) With Best Regards, Marat
6. Correlations are not considered; instead, worst case is always assumed. For example: (mx +- dx) + (my +- dy) gives ((mx + my) +- (dx + dy), not ((mx + my) +- sqrt(dx*dx + dy*dy);
If x and y are uncorrelated normally distributed variables, then the variance of x+y is the sum of the variances of x and y. This implies the latter error law, which you say you don't implement. So setting the error to (dx+dy) has nothing to do with ignoring correlations. It is hard for me to imagine a case when simply adding the errors like this is in any way meaningful. Note, in particular, that this error law would imply no benefit whatsoever to averaging the results from many observations together... you would get the same error as you would from any one trial. This is the same as assuming the quantities mx and my are Cauchy-distributed, which is a very odd assumption for the general case. -Lewis
Hello, Thank you for response. Lewis Hyatt wrote:
6. Correlations are not considered; instead, worst case is always assumed. For example: (mx +- dx) + (my +- dy) gives ((mx + my) +- (dx + dy), not ((mx + my) +- sqrt(dx*dx + dy*dy);
It is hard for me to imagine a case when simply adding the errors like this is in any way meaningful.
Isn't it the case when numbers are fully correlated? This gives largest error, this is the way I go. Well, _actual_ correlations are not considered, user has to do it manually when needed. With Best Regards, Marat
Lewis Hyatt <lhyatt@princeton.edu> writes:
6. Correlations are not considered; instead, worst case is always assumed. For example: (mx +- dx) + (my +- dy) gives ((mx + my) +- (dx + dy), not ((mx + my) +- sqrt(dx*dx + dy*dy);
If x and y are uncorrelated normally distributed variables, then the variance of x+y is the sum of the variances of x and y. This implies the latter error law, which you say you don't implement. So setting the error to (dx+dy) has nothing to do with ignoring correlations. It is hard for me to imagine a case when simply adding the errors like this is in any way meaningful. Note, in particular, that this error law would imply no benefit whatsoever to averaging the results from many observations together... you would get the same error as you would from any one trial. This is the same as assuming the quantities mx and my are Cauchy-distributed, which is a very odd assumption for the general case.
It seems that the error term may be intended to mean an exact upper/lower bound, not a variance. -- Jeremy Maitin-Shepard
Lewis Hyatt wrote:
Note, in particular, that this error law would imply no benefit whatsoever to averaging the results from many observations together... you would get the same error as you would from any one trial.
On a positive note, is this what you need for your work? Do you need function: estimation_from_statistics(list<T> &statistics) returning average and deviation? I bet one exists already, but I can add one to the library. What do you do next with the value you receive? Do you substitute it somewhere? With Best Regards, Marat
Marat Khalili wrote:
Lewis Hyatt wrote:
Note, in particular, that this error law would imply no benefit whatsoever to averaging the results from many observations together... you would get the same error as you would from any one trial.
On a positive note, is this what you need for your work? Do you need function: estimation_from_statistics(list<T> &statistics) returning average and deviation? I bet one exists already, but I can add one to the library. What do you do next with the value you receive? Do you substitute it somewhere?
Well, the boost accumulators library offers a lot of things to do this kind of calculation, but being able to keep track of the errors automatically would be nice too. Like others have pointed out, it's very hard to do this in complete generality, allowing for arbitrary correlations between all numbers. I just wanted to point out that when you see something like x +- dx, that carries a specific meaning to most people, which is that dx is roughly the 1-sigma error on x, or the standard deviation. It is true that if x and y are correlated 100%, then the error on their sum is the sum of their errors, but that is a pretty unusual case.
Lewis Hyatt wrote:
I just wanted to point out that when you see something like x +- dx, that carries a specific meaning to most people, which is that dx is roughly the 1-sigma error on x, or the standard deviation. OK, I'll try to make this distinction clear.
I'll announce when the code is ready to be seen by public. With Best Regards, Marat
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Marat Khalili Sent: 08 July 2007 19:07 To: boost@lists.boost.org Subject: [boost] [rfc] Inexact quantities arithmetics
I'm developing a library for working with inexact quantities, like (0.3162 +- 0.0079), in other words, for auto-estimating the propagation of errors in calculations. I plan to release it under Boost license when ready. Any comments (including sensibility of the whole project) are welcome. <big snip>
This all looks most interesting, and would be a big step forward in handling uncertainty - but use of the word 'error' is now so-last-millenium ;-) Joking apart, this is a serious point http://www.rsc.org/delivery/_ArticleLinking/DisplayArticleForFree.cfm?doi=a7... http://www.teknologiportalen.dk/EN/Teknologi/Events/Workshop.htm The word error gives off some very wrong vibes, and I would strongly urge you to change to using uncertainty instead of error, despite the historial use of the word error. It *really, really* isn't an error - it is a interval of (un-)certainty. You should also see http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/29189/1/95-0368.pdf // Evan Manning (manning@alumni.caltech.edu) // Evan Marshal Manning, C/C++ Users Journal, March 1996 page 29 to 38. // original downloaded from ftp://beowulf.jpl.nasa.gov/pub/manning // This is a simple model of uncertainties, designed to // accompany an article published in C/C++ Users Journal March 1996. // A fuller collection of even fancier classes also given in UReal.h. You will see from the full code that he templates two versions, one (unusual) where the uncertainties are correlated (must add up to 100%) - as well as the normal uncorrelated situation. Both are needed I think. and FWIW some of my comments some years ago http://lists.boost.org/Archives/boost/2003/04/46828.php I built my attempt without using Boost interval library, but you might be able to use it. I also felt it was vital to store the degrees of freedom with the value and uncertainty, and I also added some other bits of information like quantization (from A/D conversion), exact values, where the uncertainty estimate came from ... Paul --- Paul A Bristow Prizet Farmhouse, Kendal, Cumbria UK LA8 8AB +44 1539561830 & SMS, Mobile +44 7714 330204 & SMS pbristow@hetp.u-net.com
Hello, Paul A Bristow wrote:
The word error gives off some very wrong vibes, and I would strongly urge you to change to using uncertainty instead of error, despite the historial use of the word error. It *really, really* isn't an error - it is a interval of (un-)certainty.
Not an interval, but it's half-width, right? What do you think about 'uncertainty_interval_half_width'? Seriously I wanted to call it 'deviation' first, but opted for shorter identifier; same number of letters as in 'value' is additional bonus, makes code more symmetrical. BTW shouldn't 'value' be called average, mean or expectancy? I thought it is still colloquially being called 'error', but I'll think about it. Some help from English-speaking folks will be appreciated.
You should also see
http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/29189/1/95-0368.pdf
Yes, this one looks very much like one I'm thinking about. As a work of a government employee his code must be free, shouldn't it?
You will see from the full code that he templates two versions, one (unusual) where the uncertainties are correlated (must add up to 100%) - as well as the normal uncorrelated situation. Both are needed I think.
I'm thinking about providing a class of uncorrelated values, but return value from any operation on them will be correlated. Even this allows to underestimate the er... uncertainty if someone uses a value twice. I don't know how to prohibit this (compile time).
and FWIW some of my comments some years ago
Is it open-source?
My uncertain class has a value, standard deviation, degrees of freedom, and 16 uncertainty type flags - square, triangular, gaussian, exact etc
I also felt it was vital to store the degrees of freedom with the value and uncertainty, and I also added some other bits of information like quantization (from A/D conversion), exact values, where the uncertainty estimate came from ...
This looks too heavy for me, it was tuned for some specific task, right? Generally, adding or substracting two lognormal values produces god-knows-what. Either you keep the whole (joint!) distribution or nothing.
I built my attempt without using Boost interval library, but you might be able to use it.
I wanted to keep the main value intact. Probably value and interval is a way to go... Also, uBLAS vectors have problems even with wrapped double, I should test intervals first. With Best Regards, Marat
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Marat Khalili Sent: 15 July 2007 07:36 To: boost@lists.boost.org Subject: Re: [boost] [rfc] Inexact quantities arithmetics
Paul A Bristow wrote:
The word error gives off some very wrong vibes, and I would strongly urge you to change to using uncertainty instead of error, despite the historial use of the word error. It *really, really* isn't an error - it is a interval of (un-)certainty.
What do you think about 'uncertainty_interval_half_width'?
It could be - but that makes all sorts of assumptions about distribution, normality? 95% confidence? ... so I think you want a more vague (and shorter!) name.
Seriously I wanted to call it 'deviation' first, but opted for shorter identifier; same number of letters as in 'value' is additional bonus, makes code more symmetrical. BTW shouldn't 'value' be called average, mean or expectancy?
Well it might be, but if there is only one exact value 2 +- 0, and you might consider it 'most likely' so again I eventually favoured a vague term like value.
I thought it is still colloquially being called 'error',
You are quite right - but I strongly advise that this is out-dated.
but I'll think about it. Some help from English-speaking folks will be appreciated.
I *really* am English-speaking! - unlike Boost's many American-English speakers ;-)
You should also see
http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/29189/1/95-0368.pdf
Yes, this one looks very much like one I'm thinking about. As a work of a government employee his code must be free, shouldn't it?
I have assumed so.
and FWIW some of my comments some years ago
Is it open-source?
It isn't in a state where I would want anyone else to see it with my name on it ;-((
My uncertain class has a value, standard deviation, degrees of freedom, and 16 uncertainty type flags - square, triangular, gaussian, exact etc
This fits into 128 bits, so only doubling the space required to store a 'value'.
I also felt it was vital to store the degrees of freedom with the value and uncertainty, and I also added some other bits of information like quantization (from A/D conversion), exact values, where the uncertainty estimate came from ...
This looks too heavy for me, it was tuned for some specific task, right?
Not really, I was trying to cater for ALL types of physical measurement, from astronomy to social science.
Generally, adding or substracting two lognormal values produces god-knows-what. Either you keep the whole (joint!) distribution or nothing.
True - but I felt there was value in recording where the uncertainty estimates came from - repeated (with how many) measurements, some other estimating method, A/D conversion quantisation, round-off from number of significant digits... But I also wanted to achieve the compile-time checking of units, now a Boost library, and get output nice, including the 'right' number of significant digits (from the uncertainty). This is more that I can chew at present ;-) Though it still seems a reasonable enough requirement. Paul PS I recall some New Zealand programmers work on this, but can't find a reference immediately. --- Paul A Bristow Prizet Farmhouse, Kendal, Cumbria UK LA8 8AB +44 1539561830 & SMS, Mobile +44 7714 330204 & SMS pbristow@hetp.u-net.com
Paul A Bristow wrote:
What do you think about 'uncertainty_interval_half_width'?
It could be - but that makes all sorts of assumptions about distribution, normality? 95% confidence? ...
Not really, the library can just crunch whatever it was given by user, 95% or 66%, as long as it is consistent. Problems only arise when it is necessary to estimate 95% interval from the 66% one.
so I think you want a more vague (and shorter!) name.
Shorter! IMO 'uncertainty' is not short, and, at the same time, is neither strict nor it is immediately clear if it's width, half-width or what.
but I'll think about it. Some help from English-speaking folks will be appreciated. I *really* am English-speaking!
I just wanted to point that I'm not, and hence not that good in finding synonyms.
- unlike Boost's many American-English speakers ;-)
:)))
It isn't in a state where I would want anyone else to see it with my name on it ;-((
But so is mine. :)
My uncertain class has a value, standard deviation, degrees of freedom, and 16 uncertainty type flags - square, triangular, gaussian, exact etc
This fits into 128 bits, so only doubling the space required to store a 'value'.
I'm sure it can be very useful for storing measurement results, but AFAIU applying various arithmetic operations to a limited number of uncertainty types will create unlimited (probably even uncountable) number of more types, so implementing corresponding algebra in C++ is hardly realistic. On the other hand, AFAIU the error in estimating uncertainty if we disregard distributions completely is some constant power number of operations.
True - but I felt there was value in recording where the uncertainty estimates came from - repeated (with how many) measurements, some other estimating method, A/D conversion quantisation, round-off from number of significant digits...
Indeed it is valuable to know this for the raw data, it will help to aggregate your statistics correctly, but as soon as you aggregated it (for number of measurements >~= 30) the distribution of the estimate should resemble normal. So IMHO you are talking about a different task or stage. BTW shouldn't you have large volumes of data sharing same uncertainty type?
But I also wanted to achieve the compile-time checking of units, now a Boost library, and get output nice, including the 'right' number of significant digits (from the uncertainty).
Compile-time checking of units using templates is completely different task, and it is solved at least once AFAIK (but I don't remember the link). As for the output, I also spent some time on it, and finally more or less satisfied with the results.
PS I recall some New Zealand programmers work on this, but can't find a reference immediately.
Would be nice to see their considerations. With Best Regards, Marat
participants (6)
-
Andreas Pokorny -
Jeremy Maitin-Shepard -
Kevin Lynch -
Lewis Hyatt -
Marat Khalili -
Paul A Bristow