[boost] Automated error analysis with correlation

13 Apr 2008

      AMDG

The topic of treating error analysis has come up on this
list several times in the past.  The main difficulty is that
the uncertainties of variables are not always independent.
I recently had an idea of how to overcome this weakness.
The first assumption is that all the variables that we start with
are independent.  Then, each variable keeps track of how
it depends on each of the initial variables.

The simplest implementation is to create a unique integer id
for every new variable constructed.  The
format of a variable in the follwing discussion is
(value, [(id1, alternate_value1), (id2, alternate_value2)])
where the alternate_values are the values that we would get
if the value were increased by their uncertainties.

value_with_uncertainty x(1.2, 0.2); // holds (1.2, [(0, 1.4)])
value_with_uncertainty y(1.2, 0.2); // holds (1.2, [(1, 1.4)])

Now, adding x + x, we add the values to get (2.4, ...).  We then see
that both arguments are dependent on the input with id 0.  So,
we add these two alternate values to get (2.4, [(0, 2.8)]).  The
total uncertainty is then 0.4.

If on the other hand, we add x + y, the value will still
be 2.4, but only the first element of the list of uncertainties
depends on id 0, so we use the value of y and the adjusted value
of x yielding (2.4, [(0, 2.6), ...]).  The same rule gives us the
dependence on id 1, for a final result of (2.4, [(0, 2.6), (1, 2.6)])
Now, the combined uncertainty is sqrt(0.2*0.2 + 0.2*0.2) = 0.28.

There are several ways this can be made more accurate.
1) instead of taking x + uncertainty_of_x we can take several
    sample points inside the range x +- uncertainty_of_x
2) it should be possible to track dependence on groups
    of variables in addition to dependence on single variables
    independently.

All of this adds overhead of course, but I think that a reliable
mechanism for error propagation is often going to be worth
that cost.

The code is attached.

Comments?

In Christ,
Steven Watanabe

// value_and_uncertainty.hpp
//
// Copyright (c) 2008
// Steven Watanabe
//
// Distributed under the Boost Software License, Version 10.0 (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef VALUE_AND_UNCERTAINTY_HPP_INCLUDED
#define VALUE_AND_UNCERTAINTY_HPP_INCLUDED

#include <vector>
#include <utility>
#include <numeric>
#include <cmath>

#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>

template<class T>
class value_and_uncertainty {
public:
    value_and_uncertainty(const T& val, const T& unc) : value_(val) {
        if(unc != T()) {
            alternate_values.push_back(std::make_pair(new_id(), val + unc));
        }
    }
    template<class F>
    value_and_uncertainty(F f, const value_and_uncertainty& arg0) : value_(f(arg0.value_)) {
        typename alternate_t::const_iterator
            begin = arg0.alternate_values.begin(),
            end = arg0.alternate_values.end();
        for(; begin != end; ++begin) {
            alternate_values.push_back(std::make_pair(begin->first, f(begin->second)));
        }
    }
    template<class F>
    value_and_uncertainty(F f, const value_and_uncertainty& arg0, const value_and_uncertainty& arg1) : value_(f(arg0.value_, arg1.value_)) {
        typename alternate_t::const_iterator
            begin0 = arg0.alternate_values.begin(),
            end0 = arg0.alternate_values.end();

        typename alternate_t::const_iterator
            begin1 = arg1.alternate_values.begin(),
            end1 = arg1.alternate_values.end();

        while(true) {
            if(begin0 == end0) {
                for(; begin1 != end1; ++begin1) {
                    alternate_values.push_back(std::make_pair(begin1->first, f(arg0.value_, begin1->second)));
                }
                break;
            } else if(begin1 == end1) {
                for(; begin0 != end0; ++begin0) {
                    alternate_values.push_back(std::make_pair(begin0->first, f(begin0->second, arg1.value_)));
                }
                break;
            } else {
                if(begin0->first < begin1->first) {
                    alternate_values.push_back(std::make_pair(begin0->first, f(begin0->second, arg1.value_)));
                    ++begin0;
                } else if(begin1->first < begin0->first) {
                    alternate_values.push_back(std::make_pair(begin1->first, f(arg0.value_, begin1->second)));
                    ++begin1;
                } else {
                    alternate_values.push_back(std::make_pair(begin0->first, f(begin0->second, begin1->second)));
                    ++begin0;
                    ++begin1;
                }
            }
        }
    }
    const T& value() const {
        return(value_);
    }
    T uncertainty() const {
        namespace bll = boost::lambda;
        using namespace bll;
        using std::sqrt;
        return(sqrt(std::accumulate(alternate_values.begin(), alternate_values.end(), T(),
            _1 + bll::bind(bll::protect(_1 * _1), (bll::bind(&std::pair<id_type, T>::second, _2) - value_)))));
    }
private:
    T value_;

    typedef int id_type;
    static id_type new_id() {
        static id_type result = 0;
        return(result++);
    }
    typedef std::vector<std::pair<id_type, T> > alternate_t;
    std::vector<std::pair<id_type, T> > alternate_values;
};

template<class T>
value_and_uncertainty<T> operator+(const value_and_uncertainty<T>& arg0, const value_and_uncertainty<T>& arg1) {
    using namespace boost::lambda;
    return(value_and_uncertainty<T>(_1 + _2, arg0, arg1));
}

template<class T>
value_and_uncertainty<T> operator*(const value_and_uncertainty<T>& arg0, const value_and_uncertainty<T>& arg1) {
    using namespace boost::lambda;
    return(value_and_uncertainty<T>(_1 * _2, arg0, arg1));
}

template<class T>
std::ostream& operator<<(std::ostream& os, const value_and_uncertainty<T>& arg) {
    return(os << arg.value() << " +- " << arg.uncertainty());
}

#endif