errorprop.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_ERRORPROP_H
#define THEORETICA_ERRORPROP_H
 
#include "../core/core_traits.h"
#include "../algebra/vec.h"
#include "../autodiff/autodiff.h"
#include "../statistics/statistics.h"
#include "../pseudorandom/montecarlo.h"
#include "../pseudorandom/sampling.h"
 
 
namespace theoretica {
 
 
    namespace stats {
 
 
        template <
            typename Matrix = mat<real>, typename Dataset = vec<real>,
            enable_vector<Dataset> = true
        >
        inline Matrix covar_mat(const std::vector<Dataset>& v) {
 
            Matrix cm;
            cm.resize(v.size(), v.size());
 
            for (unsigned int i = 0; i < cm.rows(); ++i)
                for (unsigned int j = 0; j < cm.cols(); ++j)
                    cm(i, j) = stats::covariance(v[i], v[j]);
 
            return cm;
        }
 
 
        template<
            unsigned int N = 0,
            typename MultiDualFunction = autodiff::dreal_t<N>(*)(autodiff::dvec_t<N>)>
        inline real error_propagation(
            MultiDualFunction f,
            const vec<real, N>& x_best, const vec<real, N>& delta_x) {
 
            real err_sqr = 0;
            const multidual<N> df = f(multidual<N>::make_argument(x_best));
 
            for (unsigned int i = 0; i < x_best.size(); ++i)
                err_sqr += square(df.Dual().get(i) * delta_x.get(i));
 
            return sqrt(err_sqr);
        }
 
 
        template <
            unsigned int N = 0,
            typename Matrix, enable_matrix<Matrix> = true,
            typename MultiDualFunction = autodiff::dreal_t<N>(*)(autodiff::dvec_t<N>)
        >
        inline real error_propagation(
            MultiDualFunction f, const vec<real, N>& x_best, const Matrix& cm) {
 
 
            if(cm.rows() != x_best.size()) {
                TH_MATH_ERROR("error_propagation", cm.rows(), INVALID_ARGUMENT);
                return nan();
            }
 
            if(cm.cols() != x_best.size()) {
                TH_MATH_ERROR("error_propagation", cm.cols(), INVALID_ARGUMENT);
                return nan();
            }
 
            real err_sqr = 0;
            const multidual<N> df = f(multidual<N>::make_argument(x_best));
 
            for (unsigned int i = 0; i < cm.rows(); ++i)
                for (unsigned int j = 0; j < cm.cols(); ++j)
                    err_sqr += df.Dual().get(i) * df.Dual().get(j) * cm(i, j);
 
            return sqrt(err_sqr);
        }
 
 
        template<
            unsigned int N = 0,
            typename MultiDualFunction = multidual<N>(*)(autodiff::dvec_t<N>),
            typename Dataset = vec<real, N>>
        inline real error_propagation(
            MultiDualFunction f,
            const std::vector<Dataset>& v) {
 
            vec<real, N> x_mean;
            x_mean.resize(v.size());
 
            for (unsigned int i = 0; i < v.size(); ++i)
                x_mean[i] = stats::mean(v[i]);
 
            return error_propagation(
                f, x_mean, covar_mat<mat<real, N, N>, Dataset>(v));
        }
 
 
        template<typename Function>
        real error_propagation_mc(
            Function f, std::vector<pdf_sampler>& rv, unsigned int N = 1E+6) {
 
            return stats::stdev(sample_mc(f, rv, N));
        }
    }
}
 
 
#endif