MLRModel.hpp

Go to the documentation of this file.

/* 
 * MLRModel.hpp 
 *
 * This file is part of the Chemical Data Processing Toolkit
 *
 * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; see the file COPYING. If not, write to
 * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */
 
#ifndef CDPL_MATH_MLRMODEL_HPP
#define CDPL_MATH_MLRMODEL_HPP
 
#include <limits>
 
#include "CDPL/Math/Vector.hpp"
#include "CDPL/Math/Matrix.hpp"
#include "CDPL/Math/CommonType.hpp"
#include "CDPL/Math/TypeTraits.hpp"
#include "CDPL/Math/SVDecomposition.hpp"
#include "CDPL/Math/SpecialFunctions.hpp"
#include "CDPL/Base/Exceptions.hpp"
 
 
namespace CDPL
{
 
    namespace Math
    {
 
        template <typename T = double>
        class MLRModel
        {
 
          public:
            typedef typename CommonType<typename Vector<T>::SizeType, typename Matrix<T>::SizeType>::Type SizeType;
            typedef T                                                                                     ValueType;
            typedef Matrix<T>                                                                             MatrixType;
            typedef Vector<T>                                                                             VectorType;
 
            MLRModel():
                chiSquare(0), Q(0), r(0), stdDeviation(0) {}
 
            void resizeDataSet(SizeType num_points, SizeType num_vars);
 
            void clearDataSet();
 
            template <typename V>
            void setXYData(SizeType i, const VectorExpression<V>& x_vars, ValueType y);
 
            template <typename V>
            void addXYData(const VectorExpression<V>& x_vars, ValueType y);
 
            MatrixType& getXMatrix();
 
            const MatrixType& getXMatrix() const;
 
            VectorType& getYValues();
 
            const VectorType& getYValues() const;
 
            void buildModel();
 
            template <typename V>
            ValueType calcYValue(const VectorExpression<V>& x_vars) const;
 
            template <typename V>
            ValueType operator()(const VectorExpression<V>& x_vars) const;
 
            const VectorType& getCoefficients() const;
 
            ValueType getChiSquare() const;
 
            ValueType getGoodnessOfFit() const;
 
            ValueType getCorrelationCoefficient() const;
 
            ValueType getStandardDeviation() const;
 
            void calcStatistics();
 
          private:
            MatrixType xMatrix;
            VectorType yValues;
            VectorType calcYValues;
            VectorType coefficients;
            MatrixType svdU;
            MatrixType svdV;
            VectorType svdW;
            ValueType  chiSquare;
            ValueType  Q;
            ValueType  r;
            ValueType  stdDeviation;
        };
    } // namespace Math
} // namespace CDPL
 
 
// Implementation
 
template <typename T>
void CDPL::Math::MLRModel<T>::resizeDataSet(SizeType num_points, SizeType num_vars)
{
    if (num_points == xMatrix.getSize1() && num_vars == xMatrix.getSize2())
        return;
 
    xMatrix.resize(num_points, num_vars, true, ValueType());
    yValues.resize(num_points, ValueType());
}
 
template <typename T>
void CDPL::Math::MLRModel<T>::clearDataSet()
{
    yValues.resize(0);
    xMatrix.resize(0, 0, false);
}
 
template <typename T>
template <typename V>
void CDPL::Math::MLRModel<T>::setXYData(SizeType i, const VectorExpression<V>& x_vars, ValueType y)
{
    SizeType x_mtx_size1 = xMatrix.getSize1();
    SizeType x_mtx_size2 = xMatrix.getSize2();
    SizeType x_vars_size = x_vars().getSize();
    SizeType y_vals_size = yValues.getSize();
 
    resizeDataSet(std::max(i + 1, std::max(x_mtx_size1, y_vals_size)), std::max(x_vars_size, x_mtx_size2));
 
    for (SizeType j = 0; j < x_vars_size; j++)
        xMatrix(i, j) = x_vars()(j);
 
    if (x_vars_size < x_mtx_size2)
        for (SizeType j = x_vars_size; j < x_mtx_size2; j++)
            xMatrix(i, j) = ValueType();
 
    yValues(i) = y;
}
 
template <typename T>
template <typename V>
void CDPL::Math::MLRModel<T>::addXYData(const VectorExpression<V>& x_vars, ValueType y)
{
    SizeType i           = xMatrix.getSize1();
    SizeType x_mtx_size2 = xMatrix.getSize2();
    SizeType x_vars_size = x_vars().getSize();
 
    resizeDataSet(i + 1, std::max(x_mtx_size2, x_vars_size));
 
    for (SizeType j = 0; j < x_vars_size; j++)
        xMatrix(i, j) = x_vars()(j);
 
    if (x_vars_size < x_mtx_size2)
        for (SizeType j = x_vars_size; j < x_mtx_size2; j++)
            xMatrix(i, j) = ValueType();
 
    yValues(i) = y;
}
 
template <typename T>
CDPL::Math::Matrix<typename CDPL::Math::MLRModel<T>::ValueType>&
CDPL::Math::MLRModel<T>::getXMatrix()
{
    return xMatrix;
}
 
template <typename T>
const CDPL::Math::Matrix<typename CDPL::Math::MLRModel<T>::ValueType>&
CDPL::Math::MLRModel<T>::getXMatrix() const
{
    return xMatrix;
}
 
template <typename T>
CDPL::Math::Vector<typename CDPL::Math::MLRModel<T>::ValueType>&
CDPL::Math::MLRModel<T>::getYValues()
{
    return yValues;
}
 
template <typename T>
const CDPL::Math::Vector<typename CDPL::Math::MLRModel<T>::ValueType>&
CDPL::Math::MLRModel<T>::getYValues() const
{
    return yValues;
}
 
template <typename T>
void CDPL::Math::MLRModel<T>::buildModel()
{
    static const ValueType TOLERANCE(1.0e-6);
 
    SizeType n = xMatrix.getSize1();
    SizeType m = xMatrix.getSize2();
 
    if (m == 0 || n == 0)
        throw Base::CalculationFailed("MLRModel: empty data set");
 
    if (n != SizeType(yValues.getSize()))
        resizeDataSet(std::max(SizeType(yValues.getSize()), n), m);
 
    svdU.resize(n, m, false);
    svdV.resize(m, m, false);
    svdW.resize(m, false);
 
    svdU = xMatrix;
 
    if (!svDecompose(svdU, svdW, svdV))
        throw Base::CalculationFailed("MLRModel: singular value decomposition failed");
 
    ValueType w_max = svdW(0);
 
    for (SizeType i = 1; i < m; i++)
        if (svdW(i) > w_max)
            w_max = svdW(i);
 
    ValueType tresh = TOLERANCE * w_max;
 
    for (SizeType i = 0; i < m; i++)
        if (svdW(i) < tresh)
            svdW(i) = 0;
 
    coefficients.resize(m, false);
 
    svSubstitute(svdU, svdW, svdV, yValues, coefficients);
}
 
template <typename T>
template <typename V>
typename CDPL::Math::MLRModel<T>::ValueType
CDPL::Math::MLRModel<T>::calcYValue(const VectorExpression<V>& x) const
{
    if (SizeType(x().getSize()) != SizeType(coefficients.getSize()))
        throw Base::CalculationFailed("MLRModel: number of regression coefficients does not match number of independent variables");
 
    return innerProd(coefficients, x);
}
 
template <typename T>
template <typename V>
typename CDPL::Math::MLRModel<T>::ValueType
CDPL::Math::MLRModel<T>::operator()(const VectorExpression<V>& x) const
{
    return calcYValue(x);
}
 
template <typename T>
const CDPL::Math::Vector<typename CDPL::Math::MLRModel<T>::ValueType>&
CDPL::Math::MLRModel<T>::getCoefficients() const
{
    return coefficients;
}
 
template <typename T>
typename CDPL::Math::MLRModel<T>::ValueType
CDPL::Math::MLRModel<T>::getChiSquare() const
{
    return chiSquare;
}
 
template <typename T>
typename CDPL::Math::MLRModel<T>::ValueType
CDPL::Math::MLRModel<T>::getGoodnessOfFit() const
{
    return Q;
}
 
template <typename T>
typename CDPL::Math::MLRModel<T>::ValueType
CDPL::Math::MLRModel<T>::getCorrelationCoefficient() const
{
    return r;
}
 
template <typename T>
typename CDPL::Math::MLRModel<T>::ValueType
CDPL::Math::MLRModel<T>::getStandardDeviation() const
{
    return stdDeviation;
}
 
template <typename T>
void CDPL::Math::MLRModel<T>::calcStatistics()
{
    ValueType mean_data_y = 0;
    ValueType mean_calc_y = 0;
 
    chiSquare = ValueType();
 
    SizeType m = xMatrix.getSize2();
    SizeType n = xMatrix.getSize1();
 
    if (m != SizeType(coefficients.getSize()))
        throw Base::CalculationFailed("MLRModel: number of independent variables does not match number of regression coefficients");
 
    if (n != SizeType(yValues.getSize()))
        throw Base::CalculationFailed("MLRModel: number of dependent variables does not match number of vectors with independent variables");
 
    calcYValues.resize(n);
 
    for (SizeType i = 0; i < n; i++) {
        ValueType data_y = yValues(i);
        ValueType calc_y = innerProd(row(xMatrix, i), coefficients);
        ValueType y_diff = data_y - calc_y;
 
        mean_data_y += data_y;
        mean_calc_y += calc_y;
        chiSquare += y_diff * y_diff;
 
        calcYValues(i) = calc_y;
    }
 
    mean_data_y /= n;
    mean_calc_y /= n;
 
    ValueType sxx = ValueType();
    ValueType syy = ValueType();
    ValueType sxy = ValueType();
 
    for (SizeType i = 0; i < n; i++) {
        ValueType xt = yValues(i) - mean_data_y;
        ValueType yt = calcYValues(i) - mean_calc_y;
 
        sxx += xt * xt;
        syy += yt * yt;
        sxy += xt * yt;
    }
 
    r = sxy / (TypeTraits<ValueType>::sqrt(sxx * syy) + std::numeric_limits<ValueType>::epsilon());
 
    stdDeviation = TypeTraits<ValueType>::sqrt(chiSquare / (n - m));
 
    Q = gammaQ(ValueType(n - 2) / 2, chiSquare / 2);
}
 
#endif // CDPL_MATH_MLRMODEL_HPP