diff --git a/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp b/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp index 3ab35610a..f50ccf1b1 100644 --- a/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp +++ b/computer_oriented_statistical_methods/ordinary_least_squares_regressor.cpp @@ -1,280 +1,279 @@ -#include +/** + * Program that gets the number of data samples and number of features per + * sample along with output per sample. It applies OLS regression to compute + * the regression output for additional test data samples. + **/ #include #include -using namespace std; +#include template -ostream &operator<<(ostream &out, vector> const &v) -{ - const int width = 10; - const char separator = ' '; +std::ostream &operator<<(std::ostream &out, + std::vector> const &v) { + const int width = 10; + const char separator = ' '; - for (size_t row = 0; row < v.size(); row++) - { - for (size_t col = 0; col < v[row].size(); col++) - out - << left << setw(width) << setfill(separator) << v[row][col]; - out << endl; - } + for (size_t row = 0; row < v.size(); row++) { + for (size_t col = 0; col < v[row].size(); col++) + out << std::left << std::setw(width) << std::setfill(separator) + << v[row][col]; + out << std::endl; + } - return out; + return out; } template -ostream &operator<<(ostream &out, vector const &v) -{ - const int width = 15; - const char separator = ' '; +std::ostream &operator<<(std::ostream &out, std::vector const &v) { + const int width = 15; + const char separator = ' '; - for (size_t row = 0; row < v.size(); row++) - out - << left << setw(width) << setfill(separator) << v[row]; + for (size_t row = 0; row < v.size(); row++) + out << std::left << std::setw(width) << std::setfill(separator) << v[row]; - return out; + return out; } template -inline bool is_square(vector> const &A) -{ - size_t N = A.size(); // Assuming A is square matrix - for (size_t i = 0; i < N; i++) - if (A[i].size() != N) - return false; - return true; +inline bool is_square(std::vector> const &A) { + // Assuming A is square matrix + size_t N = A.size(); + for (size_t i = 0; i < N; i++) + if (A[i].size() != N) + return false; + return true; } /** * matrix multiplication **/ template -vector> operator*(vector> const &A, vector> const &B) -{ - size_t N_A = A.size(); // Number of rows in A - size_t N_B = B[0].size(); // Number of columns in B +std::vector> operator*(std::vector> const &A, + std::vector> const &B) { + // Number of rows in A + size_t N_A = A.size(); + // Number of columns in B + size_t N_B = B[0].size(); - vector> result(N_A); - - if (A[0].size() != B.size()) - { - cerr << "Number of columns in A != Number of rows in B (" << A[0].size() << ", " << B.size() << ")" << endl; - return result; - } - - for (size_t row = 0; row < N_A; row++) - { - vector v(N_B); - for (size_t col = 0; col < N_B; col++) - { - v[col] = static_cast(0); - for (size_t j = 0; j < B.size(); j++) - v[col] += A[row][j] * B[j][col]; - } - result[row] = v; - } + std::vector> result(N_A); + if (A[0].size() != B.size()) { + std::cerr << "Number of columns in A != Number of rows in B (" + << A[0].size() << ", " << B.size() << ")" << std::endl; return result; + } + + for (size_t row = 0; row < N_A; row++) { + std::vector v(N_B); + for (size_t col = 0; col < N_B; col++) { + v[col] = static_cast(0); + for (size_t j = 0; j < B.size(); j++) + v[col] += A[row][j] * B[j][col]; + } + result[row] = v; + } + + return result; } template -vector operator*(vector> const &A, vector const &B) -{ - size_t N_A = A.size(); // Number of rows in A +std::vector operator*(std::vector> const &A, + std::vector const &B) { + // Number of rows in A + size_t N_A = A.size(); - vector result(N_A); - - if (A[0].size() != B.size()) - { - cerr << "Number of columns in A != Number of rows in B (" << A[0].size() << ", " << B.size() << ")" << endl; - return result; - } - - for (size_t row = 0; row < N_A; row++) - { - result[row] = static_cast(0); - for (size_t j = 0; j < B.size(); j++) - result[row] += A[row][j] * B[j]; - } + std::vector result(N_A); + if (A[0].size() != B.size()) { + std::cerr << "Number of columns in A != Number of rows in B (" + << A[0].size() << ", " << B.size() << ")" << std::endl; return result; + } + + for (size_t row = 0; row < N_A; row++) { + result[row] = static_cast(0); + for (size_t j = 0; j < B.size(); j++) + result[row] += A[row][j] * B[j]; + } + + return result; } template -vector operator*(float const scalar, vector const &A) -{ - size_t N_A = A.size(); // Number of rows in A +std::vector operator*(float const scalar, std::vector const &A) { + // Number of rows in A + size_t N_A = A.size(); - vector result(N_A); + std::vector result(N_A); - for (size_t row = 0; row < N_A; row++) - { - result[row] += A[row] * static_cast(scalar); - } + for (size_t row = 0; row < N_A; row++) { + result[row] += A[row] * static_cast(scalar); + } - return result; + return result; } template -vector operator*(vector const &A, float const scalar) -{ - size_t N_A = A.size(); // Number of rows in A +std::vector operator*(std::vector const &A, float const scalar) { + // Number of rows in A + size_t N_A = A.size(); - vector result(N_A); + std::vector result(N_A); - for (size_t row = 0; row < N_A; row++) - result[row] = A[row] * static_cast(scalar); + for (size_t row = 0; row < N_A; row++) + result[row] = A[row] * static_cast(scalar); - return result; + return result; } template -vector operator/(vector const &A, float const scalar) -{ - return (1.f / scalar) * A; +std::vector operator/(std::vector const &A, float const scalar) { + return (1.f / scalar) * A; } template -vector operator-(vector const &A, vector const &B) -{ - size_t N = A.size(); // Number of rows in A +std::vector operator-(std::vector const &A, std::vector const &B) { + // Number of rows in A + size_t N = A.size(); - vector result(N); + std::vector result(N); - if (B.size() != N) - { - cerr << "Vector dimensions shouldbe identical!" << endl; - return A; - } + if (B.size() != N) { + std::cerr << "Vector dimensions shouldbe identical!" << std::endl; + return A; + } - for (size_t row = 0; row < N; row++) - result[row] = A[row] - B[row]; + for (size_t row = 0; row < N; row++) + result[row] = A[row] - B[row]; - return result; + return result; } template -vector operator+(vector const &A, vector const &B) -{ - size_t N = A.size(); // Number of rows in A +std::vector operator+(std::vector const &A, std::vector const &B) { + // Number of rows in A + size_t N = A.size(); - vector result(N); + std::vector result(N); - if (B.size() != N) - { - cerr << "Vector dimensions shouldbe identical!" << endl; - return A; - } + if (B.size() != N) { + std::cerr << "Vector dimensions shouldbe identical!" << std::endl; + return A; + } - for (size_t row = 0; row < N; row++) - result[row] = A[row] + B[row]; + for (size_t row = 0; row < N; row++) + result[row] = A[row] + B[row]; - return result; + return result; } /** * Get matrix inverse using Row-trasnformations **/ template -vector> get_inverse(vector> const &A) -{ - size_t N = A.size(); // Assuming A is square matrix +std::vector> +get_inverse(std::vector> const &A) { + // Assuming A is square matrix + size_t N = A.size(); - vector> inverse(N); - for (size_t row = 0; row < N; row++) // preallocatae a resultant identity matrix - { - inverse[row] = vector(N); - for (size_t col = 0; col < N; col++) - inverse[row][col] = (row == col) ? 1.f : 0.f; - } - - if (!is_square(A)) - { - cerr << "A must be a square matrix!" << endl; - return inverse; - } - - vector> temp(N); // preallocatae a temporary matrix identical to A - for (size_t row = 0; row < N; row++) - { - vector v(N); - for (size_t col = 0; col < N; col++) - v[col] = static_cast(A[row][col]); - temp[row] = v; - } - - // start transformations - for (size_t row = 0; row < N; row++) - { - for (size_t row2 = row; row2 < N && temp[row][row] == 0; row2++) // this to ensure diagonal elements are not 0 - { - temp[row] = temp[row] + temp[row2]; - inverse[row] = inverse[row] + inverse[row2]; - } - - for (size_t col2 = row; col2 < N && temp[row][row] == 0; col2++) // this to further ensure diagonal elements are not 0 - { - for (size_t row2 = 0; row2 < N; row2++) - { - temp[row2][row] = temp[row2][row] + temp[row2][col2]; - inverse[row2][row] = inverse[row2][row] + inverse[row2][col2]; - } - } - - if (temp[row][row] == 0) - { - // Probably a low-rank matrix and hence singular - cerr << "Low-rank matrix, no inverse!" << endl; - return inverse; - } - - // set diagonal to 1 - float divisor = static_cast(temp[row][row]); - temp[row] = temp[row] / divisor; - inverse[row] = inverse[row] / divisor; - // Row transformations - for (size_t row2 = 0; row2 < N; row2++) - { - if (row2 == row) - continue; - float factor = temp[row2][row]; - temp[row2] = temp[row2] - factor * temp[row]; - inverse[row2] = inverse[row2] - factor * inverse[row]; - } - } + std::vector> inverse(N); + for (size_t row = 0; row < N; row++) { + // preallocatae a resultant identity matrix + inverse[row] = std::vector(N); + for (size_t col = 0; col < N; col++) + inverse[row][col] = (row == col) ? 1.f : 0.f; + } + if (!is_square(A)) { + std::cerr << "A must be a square matrix!" << std::endl; return inverse; + } + + // preallocatae a temporary matrix identical to A + std::vector> temp(N); + for (size_t row = 0; row < N; row++) { + std::vector v(N); + for (size_t col = 0; col < N; col++) + v[col] = static_cast(A[row][col]); + temp[row] = v; + } + + // start transformations + for (size_t row = 0; row < N; row++) { + for (size_t row2 = row; row2 < N && temp[row][row] == 0; row2++) { + // this to ensure diagonal elements are not 0 + temp[row] = temp[row] + temp[row2]; + inverse[row] = inverse[row] + inverse[row2]; + } + + for (size_t col2 = row; col2 < N && temp[row][row] == 0; col2++) { + // this to further ensure diagonal elements are not 0 + for (size_t row2 = 0; row2 < N; row2++) { + temp[row2][row] = temp[row2][row] + temp[row2][col2]; + inverse[row2][row] = inverse[row2][row] + inverse[row2][col2]; + } + } + + if (temp[row][row] == 0) { + // Probably a low-rank matrix and hence singular + std::cerr << "Low-rank matrix, no inverse!" << std::endl; + return inverse; + } + + // set diagonal to 1 + float divisor = static_cast(temp[row][row]); + temp[row] = temp[row] / divisor; + inverse[row] = inverse[row] / divisor; + // Row transformations + for (size_t row2 = 0; row2 < N; row2++) { + if (row2 == row) + continue; + float factor = temp[row2][row]; + temp[row2] = temp[row2] - factor * temp[row]; + inverse[row2] = inverse[row2] - factor * inverse[row]; + } + } + + return inverse; } /** * matrix transpose **/ template -vector> get_transpose(vector> const &A) -{ - vector> result(A[0].size()); +std::vector> +get_transpose(std::vector> const &A) { + std::vector> result(A[0].size()); - for (size_t row = 0; row < A[0].size(); row++) - { - vector v(A.size()); - for (size_t col = 0; col < A.size(); col++) - v[col] = A[col][row]; + for (size_t row = 0; row < A[0].size(); row++) { + std::vector v(A.size()); + for (size_t col = 0; col < A.size(); col++) + v[col] = A[col][row]; - result[row] = v; - } - return result; + result[row] = v; + } + return result; } template -vector fit_OLS_regressor(vector> const &X, vector const &Y) -{ - vector> X2 = X; //NxF - for (size_t i = 0; i < X2.size(); i++) - X2[i].push_back(1); // add Y-intercept -> Nx(F+1) - vector> Xt = get_transpose(X2); // (F+1)xN - vector> tmp = get_inverse(Xt * X2); // (F+1)x(F+1) - vector> out = tmp * Xt; // (F+1)xN - // cout << endl - // << "Projection matrix: " << X2 * out << endl; - return out * Y; // Fx1,1 -> (F+1)^th element is the independent constant +std::vector fit_OLS_regressor(std::vector> const &X, + std::vector const &Y) { + // NxF + std::vector> X2 = X; + for (size_t i = 0; i < X2.size(); i++) + // add Y-intercept -> Nx(F+1) + X2[i].push_back(1); + // (F+1)xN + std::vector> Xt = get_transpose(X2); + // (F+1)x(F+1) + std::vector> tmp = get_inverse(Xt * X2); + // (F+1)xN + std::vector> out = tmp * Xt; + // cout << endl + // << "Projection matrix: " << X2 * out << endl; + + // Fx1,1 -> (F+1)^th element is the independent constant + return out * Y; } /** @@ -282,59 +281,69 @@ vector fit_OLS_regressor(vector> const &X, vector const &Y) * regression estimates **/ template -vector predict_OLS_regressor(vector> const &X, vector const &beta) -{ - vector result(X.size()); +std::vector predict_OLS_regressor(std::vector> const &X, + std::vector const &beta) { + std::vector result(X.size()); - for (size_t rows = 0; rows < X.size(); rows++) - { - result[rows] = beta[X[0].size()]; // -> start with constant term - for (size_t cols = 0; cols < X[0].size(); cols++) - result[rows] += beta[cols] * X[rows][cols]; - } - return result; // Nx1 + for (size_t rows = 0; rows < X.size(); rows++) { + // -> start with constant term + result[rows] = beta[X[0].size()]; + for (size_t cols = 0; cols < X[0].size(); cols++) + result[rows] += beta[cols] * X[rows][cols]; + } + // Nx1 + return result; } -int main() -{ - size_t N, F; +int main() { + size_t N, F; - cin >> F; // number of features = columns - cin >> N; // number of samples = rows + std::cout << "Enter number of features: "; + // number of features = columns + std::cin >> F; + std::cout << "Enter number of samples: "; + // number of samples = rows + std::cin >> N; - vector> data(N); - vector Y(N); + std::vector> data(N); + std::vector Y(N); - for (size_t rows = 0; rows < N; rows++) - { - vector v(F); - for (size_t cols = 0; cols < F; cols++) - cin >> v[cols]; // get the F features - data[rows] = v; - cin >> Y[rows]; // get the corresponding output - } + std::cout + << "Enter training data. Per sample, provide features ad one output." + << std::endl; - vector beta = fit_OLS_regressor(data, Y); - cout << endl - << endl - << "beta:" << beta << endl; + for (size_t rows = 0; rows < N; rows++) { + std::vector v(F); + std::cout << "Sample# " << rows + 1 << ": "; + for (size_t cols = 0; cols < F; cols++) + // get the F features + std::cin >> v[cols]; + data[rows] = v; + // get the corresponding output + std::cin >> Y[rows]; + } - size_t T; - cin >> T; // number of test sample inputs - vector> data2(T); - // vector Y2(T); + std::vector beta = fit_OLS_regressor(data, Y); + std::cout << std::endl << std::endl << "beta:" << beta << std::endl; - for (size_t rows = 0; rows < T; rows++) - { - vector v(F); - for (size_t cols = 0; cols < F; cols++) - cin >> v[cols]; - data2[rows] = v; - } + size_t T; + std::cout << "Enter number of test samples: "; + // number of test sample inputs + std::cin >> T; + std::vector> data2(T); + // vector Y2(T); - vector out = predict_OLS_regressor(data2, beta); - for (size_t rows = 0; rows < T; rows++) - cout << out[rows] << endl; + for (size_t rows = 0; rows < T; rows++) { + std::cout << "Sample# " << rows + 1 << ": "; + std::vector v(F); + for (size_t cols = 0; cols < F; cols++) + std::cin >> v[cols]; + data2[rows] = v; + } - return 0; + std::vector out = predict_OLS_regressor(data2, beta); + for (size_t rows = 0; rows < T; rows++) + std::cout << out[rows] << std::endl; + + return 0; }