Task

Given a set of data vectors in the following format:

   
  
    
      
        y
        =
        {
        
          y
          
            1
          
        
        ,
        
          y
          
            2
          
        
        ,
        .
        .
        .
        ,
        
          y
          
            n
          
        
        }
        
      
    
    {\displaystyle y=\{y_{1},y_{2},...,y_{n}\}\,}
  


   
  
    
      
        
          X
          
            i
          
        
        =
        {
        
          x
          
            i
            1
          
        
        ,
        
          x
          
            i
            2
          
        
        ,
        .
        .
        .
        ,
        
          x
          
            i
            n
          
        
        }
        ,
        i
        ∈
        1..
        k
        
      
    
    {\displaystyle X_{i}=\{x_{i1},x_{i2},...,x_{in}\},i\in 1..k\,}
  


Compute the vector 



β
=
{

β

1


,

β

2


,
.
.
.
,

β

k


}


{\displaystyle \beta =\{\beta _{1},\beta _{2},...,\beta _{k}\}}

 using ordinary least squares regression using the following equation:

   
  
    
      
        
          y
          
            j
          
        
        =
        
          Σ
          
            i
          
        
        
          β
          
            i
          
        
        ⋅
        
          x
          
            i
            j
          
        
        ,
        j
        ∈
        1..
        n
      
    
    {\displaystyle y_{j}=\Sigma _{i}\beta _{i}\cdot x_{ij},j\in 1..n}
  


You can assume  y  is given to you as a vector (a one-dimensional array), and  X  is given to you as a two-dimensional array (i.e. matrix).




#include <array>
#include <iostream>
 
void require(bool condition, const std::string &message) {
    if (condition) {
        return;
    }
    throw std::runtime_error(message);
}
 
template<typename T, size_t N>
std::ostream &operator<<(std::ostream &os, const std::array<T, N> &a) {
    auto it = a.cbegin();
    auto end = a.cend();
 
    os << '[';
    if (it != end) {
        os << *it;
        it = std::next(it);
    }
    while (it != end) {
        os << ", " << *it;
        it = std::next(it);
    }
    return os << ']';
}
 
template <size_t RC, size_t CC>
class Matrix {
    std::array<std::array<double, CC>, RC> data;
 
public:
    Matrix() : data{} {
        // empty
    }
 
    Matrix(std::initializer_list<std::initializer_list<double>> values) {
        size_t rp = 0;
        for (auto row : values) {
            size_t cp = 0;
            for (auto col : row) {
                data[rp][cp] = col;
                cp++;
            }
            rp++;
        }
    }
 
    double get(size_t row, size_t col) const {
        return data[row][col];
    }
 
    void set(size_t row, size_t col, double value) {
        data[row][col] = value;
    }
 
    std::array<double, CC> get(size_t row) {
        return data[row];
    }
 
    void set(size_t row, const std::array<double, CC> &values) {
        std::copy(values.begin(), values.end(), data[row].begin());
    }
 
    template <size_t D>
    Matrix<RC, D> operator*(const Matrix<CC, D> &rhs) const {
        Matrix<RC, D> result;
        for (size_t i = 0; i < RC; i++) {
            for (size_t j = 0; j < D; j++) {
                for (size_t k = 0; k < CC; k++) {
                    double prod = get(i, k) * rhs.get(k, j);
                    result.set(i, j, result.get(i, j) + prod);
                }
            }
        }
        return result;
    }
 
    Matrix<CC, RC> transpose() const {
        Matrix<CC, RC> trans;
        for (size_t i = 0; i < RC; i++) {
            for (size_t j = 0; j < CC; j++) {
                trans.set(j, i, data[i][j]);
            }
        }
        return trans;
    }
 
    void toReducedRowEchelonForm() {
        size_t lead = 0;
        for (size_t r = 0; r < RC; r++) {
            if (CC <= lead) {
                return;
            }
            auto i = r;
 
            while (get(i, lead) == 0.0) {
                i++;
                if (RC == i) {
                    i = r;
                    lead++;
                    if (CC == lead) {
                        return;
                    }
                }
            }
 
            auto temp = get(i);
            set(i, get(r));
            set(r, temp);
 
            if (get(r, lead) != 0.0) {
                auto div = get(r, lead);
                for (size_t j = 0; j < CC; j++) {
                    set(r, j, get(r, j) / div);
                }
            }
 
            for (size_t k = 0; k < RC; k++) {
                if (k != r) {
                    auto mult = get(k, lead);
                    for (size_t j = 0; j < CC; j++) {
                        auto prod = get(r, j) * mult;
                        set(k, j, get(k, j) - prod);
                    }
                }
            }
 
            lead++;
        }
    }
 
    Matrix<RC, RC> inverse() {
        require(RC == CC, "Not a square matrix");
 
        Matrix<RC, 2 * RC> aug;
        for (size_t i = 0; i < RC; i++) {
            for (size_t j = 0; j < RC; j++) {
                aug.set(i, j, get(i, j));
            }
            // augment identify matrix to right
            aug.set(i, i + RC, 1.0);
        }
 
        aug.toReducedRowEchelonForm();
 
        // remove identity matrix to left
        Matrix<RC, RC> inv;
        for (size_t i = 0; i < RC; i++) {
            for (size_t j = RC; j < 2 * RC; j++) {
                inv.set(i, j - RC, aug.get(i, j));
            }
        }
        return inv;
    }
 
    template <size_t RC, size_t CC>
    friend std::ostream &operator<<(std::ostream &, const Matrix<RC, CC> &);
};
 
template <size_t RC, size_t CC>
std::ostream &operator<<(std::ostream &os, const Matrix<RC, CC> &m) {
    for (size_t i = 0; i < RC; i++) {
        os << '[';
        for (size_t j = 0; j < CC; j++) {
            if (j > 0) {
                os << ", ";
            }
            os << m.get(i, j);
        }
        os << "]\n";
    }
 
    return os;
}
 
template <size_t RC, size_t CC>
std::array<double, RC> multiple_regression(const std::array<double, CC> &y, const Matrix<RC, CC> &x) {
    Matrix<1, CC> tm;
    tm.set(0, y);
 
    auto cy = tm.transpose();
    auto cx = x.transpose();
    return ((x * cx).inverse() * x * cy).transpose().get(0);
}
 
void case1() {
    std::array<double, 5> y{ 1.0, 2.0, 3.0, 4.0, 5.0 };
    Matrix<1, 5> x{ {2.0, 1.0, 3.0, 4.0, 5.0} };
    auto v = multiple_regression(y, x);
    std::cout << v << '\n';
}
 
void case2() {
    std::array<double, 3> y{ 3.0, 4.0, 5.0 };
    Matrix<2, 3> x{
        {1.0, 2.0, 1.0},
        {1.0, 1.0, 2.0}
    };
    auto v = multiple_regression(y, x);
    std::cout << v << '\n';
}
 
void case3() {
    std::array<double, 15> y{ 52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46 };
    std::array<double, 15> a{ 1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83 };
 
    Matrix<3, 15> x;
    for (size_t i = 0; i < 15; i++) {
        x.set(0, i, 1.0);
    }
    x.set(1, a);
    for (size_t i = 0; i < 15; i++) {
        x.set(2, i, a[i] * a[i]);
    }
 
    auto v = multiple_regression(y, x);
    std::cout << v << '\n';
}
 
int main() {
    case1();
    case2();
    case3();
 
    return 0;
}