fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. #define ROWS 29
  6. #define COLS 5
  7.  
  8. // Function to multiply two matrices
  9. void matrix_multiply(double a[ROWS][COLS], double b[COLS][ROWS], double result[COLS][COLS]) {
  10.     for (int i = 0; i < COLS; ++i) {
  11.         for (int j = 0; j < COLS; ++j) {
  12.             result[i][j] = 0;
  13.             for (int k = 0; k < ROWS; ++k) {
  14.                 result[i][j] += a[k][i] * b[j][k];
  15.             }
  16.         }
  17.     }
  18. }
  19.  
  20. // Function to transpose a matrix
  21. void matrix_transpose(double src[ROWS][COLS], double dest[COLS][ROWS]) {
  22.     for (int i = 0; i < ROWS; ++i) {
  23.         for (int j = 0; j < COLS; ++j) {
  24.             dest[j][i] = src[i][j];
  25.         }
  26.     }
  27. }
  28.  
  29. // Function to invert a 5x5 matrix (using Gauss-Jordan elimination)
  30. int matrix_inverse(double src[COLS][COLS], double dest[COLS][COLS]) {
  31.     int i, j, k;
  32.     double ratio;
  33.     double aug[COLS][2 * COLS];
  34.  
  35.     // Initialize augmented matrix
  36.     for (i = 0; i < COLS; i++) {
  37.         for (j = 0; j < COLS; j++) {
  38.             aug[i][j] = src[i][j];
  39.             aug[i][j + COLS] = (i == j) ? 1 : 0;
  40.         }
  41.     }
  42.  
  43.     // Apply Gauss-Jordan elimination
  44.     for (i = 0; i < COLS; i++) {
  45.         if (aug[i][i] == 0) {
  46.             // Find a row below the current row with a non-zero element in the current column
  47.             for (k = i + 1; k < COLS; k++) {
  48.                 if (aug[k][i] != 0) {
  49.                     for (j = 0; j < 2 * COLS; j++) {
  50.                         double temp = aug[i][j];
  51.                         aug[i][j] = aug[k][j];
  52.                         aug[k][j] = temp;
  53.                     }
  54.                     break;
  55.                 }
  56.             }
  57.         }
  58.  
  59.         // Make the diagonal contain all ones
  60.         for (j = 0; j < COLS; j++) {
  61.             aug[i][j] /= aug[i][i];
  62.             aug[i][j + COLS] /= aug[i][i];
  63.         }
  64.  
  65.         // Make the other rows contain zeros in the current column
  66.         for (k = 0; k < COLS; k++) {
  67.             if (k != i) {
  68.                 ratio = aug[k][i];
  69.                 for (j = 0; j < 2 * COLS; j++) {
  70.                     aug[k][j] -= ratio * aug[i][j];
  71.                 }
  72.             }
  73.         }
  74.     }
  75.  
  76.     // Copy the inverse matrix
  77.     for (i = 0; i < COLS; i++) {
  78.         for (j = 0; j < COLS; j++) {
  79.             dest[i][j] = aug[i][j + COLS];
  80.         }
  81.     }
  82.  
  83.     return 0;
  84. }
  85.  
  86. // Function to solve the linear system W = inv(X^T * X) * X^T * y
  87. void matrix_vector_multiply(double matrix[COLS][COLS], double vector[COLS], double result[COLS]) {
  88.     for (int i = 0; i < COLS; ++i) {
  89.         result[i] = 0;
  90.         for (int j = 0; j < COLS; ++j) {
  91.             result[i] += matrix[i][j] * vector[j];
  92.         }
  93.     }
  94. }
  95.  
  96. int main() {
  97.     // Provided dataset
  98.     double data[ROWS][3] = {
  99.         {3.85, 3.8, 7.38}, {3.75, 4, 8.51}, {3.7, 4.3, 9.52}, {3.7, 3.7, 7.5}, {3.6, 3.85, 9.33},
  100.         {3.6, 3.8, 8.28}, {3.6, 3.75, 8.75}, {3.8, 3.85, 7.87}, {3.8, 3.65, 7.1}, {3.85, 4, 8},
  101.         {3.9, 4.1, 7.89}, {3.9, 4, 8.15}, {3.7, 4.1, 9.1}, {3.75, 4.2, 8.86}, {3.75, 4.1, 8.9},
  102.         {3.8, 4.1, 8.87}, {3.7, 4.3, 9}, {3.7, 4.1, 8.75}, {3.8, 3.75, 7.95}, {3.8, 3.75, 7.65},
  103.         {3.75, 3.65, 7.27}, {3.7, 3.9, 8}, {3.55, 3.65, 8.5}, {3.6, 4.1, 8.75}, {3.65, 4.25, 9.21},
  104.         {3.7, 3.65, 8.27}, {3.75, 3.75, 7.67}, {3.8, 3.85, 7.93}, {3.7, 4.25, 9.26}
  105.     };
  106.  
  107.     // Initialize matrices
  108.     double X[ROWS][COLS];
  109.     double Xt[COLS][ROWS];
  110.     double XtX[COLS][COLS];
  111.     double XtX_inv[COLS][COLS];
  112.     double y[ROWS];
  113.     double Xty[COLS];
  114.     double W[COLS];
  115.  
  116.     // Populate matrices X and vector y
  117.     for (int i = 0; i < ROWS; ++i) {
  118.         X[i][0] = 1.0; // First column of ones
  119.         X[i][1] = data[i][0]; // x1 column
  120.         X[i][2] = data[i][1]; // x2 column
  121.         X[i][3] = data[i][0] * data[i][0]; // x1^2 column
  122.         X[i][4] = data[i][1] * data[i][1]; // x2^2 column
  123.         y[i] = data[i][2];
  124.     }
  125.  
  126.     // Transpose of X
  127.     matrix_transpose(X, Xt);
  128.  
  129.     // Xt * X
  130.     matrix_multiply(Xt, X, XtX);
  131.  
  132.     // Inverse of Xt * X
  133.     matrix_inverse(XtX, XtX_inv);
  134.  
  135.     // Xt * y
  136.     matrix_vector_multiply(Xt, y, Xty);
  137.  
  138.     // W = inv(XtX) * Xty
  139.     matrix_vector_multiply(XtX_inv, Xty, W);
  140.  
  141.     // Print the result
  142.     printf("Calculated W^T:\n");
  143.     for (int i = 0; i < COLS; ++i) {
  144.         printf("%lf ", W[i]);
  145.     }
  146.     printf("\n");
  147.  
  148.     return 0;
  149. }
  150.  
Success #stdin #stdout 0.01s 5284KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Calculated W^T:
-63962441432531955808425013131559776695487310069760.000000 39723425346549572065586213363881294015894650880.000000 14658701178612767788804270066180028366848.000000 2337792211430634421483470848.000000 42.240000
https://ideone.com/jyggx1
language:
C (gcc 8.3)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!