線形回帰は、複数の変数間の定量的関係をモデル化する統計的手法である。目的変数（従属変数）と一つ以上の説明変数（独立変数）の間の線形関係を、最小二乗法を用いて推定する。本稿では、C言語を用いて単回帰と重回帰を実装する方法を示す。

単回帰モデル

単回帰は、一つの説明変数 x と一つの目的変数 y の関係を直線 y = ax + b で近似する。最小二乗法により、残差の二乗和を最小化する係数 a（傾き）と b（切片）を求める。

与えられたデータ点集合 (x_i, y_i) に対して、以下の式で係数を計算する：

• 傾き: a = (n·Σx_iy_i - Σx_i·Σy_i) / (n·Σx_i² - (Σx_i)²)
• 切片: b = (Σx_i²·Σy_i - Σx_i·Σx_iy_i) / (n·Σx_i² - (Σx_i)²)

以下は、温度とアイスクリーム売上を含むCSVファイル（ice_cream.csv）からデータを読み取り、回帰係数を計算するCプログラムである：

#include <stdio.h>

int main() {
    FILE *fp = fopen("ice_cream.csv", "r");
    if (!fp) {
        perror("ファイルを開けません");
        return 1;
    }

    double temp, sales;
    double sum_x = 0.0, sum_y = 0.0, sum_xy = 0.0, sum_xx = 0.0;
    int count = 0;
    char line[1024];

    // ヘッダー行をスキップ
    fgets(line, sizeof(line), fp);

    while (fgets(line, sizeof(line), fp)) {
        if (sscanf(line, "%lf,%lf", &temp, &sales) == 2) {
            sum_x += temp;
            sum_y += sales;
            sum_xy += temp * sales;
            sum_xx += temp * temp;
            count++;
        }
    }

    fclose(fp);

    if (count < 2) {
        printf("データが不足しています\n");
        return 1;
    }

    double slope = (count * sum_xy - sum_x * sum_y) / (count * sum_xx - sum_x * sum_x);
    double intercept = (sum_xx * sum_y - sum_x * sum_xy) / (count * sum_xx - sum_x * sum_x);

    printf("傾き: %.6f\n", slope);
    printf("切片: %.6f\n", intercept);

    return 0;
}

実行結果から得られる回帰式は、売上 = 3.8394 × 温度 - 253.17 となる。

重回帰モデル

重回帰では、複数の説明変数（例：温度、観光客数、晴天日数）と目的変数との線形関係をモデル化する。モデルは以下のように表される：

y = β₀ + β₁x₁ + β₂x₂ + β₃x₃ + ε

ここで、β₀ は切片、β₁〜β₃ は各説明変数の回帰係数、ε は誤差項である。

係数は、設計行列 X と目的変数ベクトル Y を用いて、正規方程式 β = (XᵀX)⁻¹XᵀY を解くことで求められる。この逆行列計算には、ガウス・ジョルダン法を用いた行基本変形を適用する。

以下の実装では、CSVファイル（ice_cream2.csv）からデータを読み込み、設計行列と目的変数行列を構築し、行列演算を介して回帰係数を算出する。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define MAX_SAMPLES 100
#define N_FEATURES 3  // 温度、観光客数、晴天日数

typedef struct {
    double temperature;
    double tourists;
    double sunny_days;
    double sales;
} Sample;

typedef struct {
    double **data;
    int rows;
    int cols;
} Matrix;

Matrix *create_matrix(int r, int c) {
    Matrix *m = malloc(sizeof(Matrix));
    m->rows = r;
    m->cols = c;
    m->data = malloc(r * sizeof(double *));
    for (int i = 0; i < r; i++) {
        m->data[i] = calloc(c, sizeof(double));
    }
    return m;
}

void destroy_matrix(Matrix *m) {
    if (!m) return;
    for (int i = 0; i < m->rows; i++) {
        free(m->data[i]);
    }
    free(m->data);
    free(m);
}

Matrix *transpose(Matrix *m) {
    Matrix *t = create_matrix(m->cols, m->rows);
    for (int i = 0; i < m->rows; i++) {
        for (int j = 0; j < m->cols; j++) {
            t->data[j][i] = m->data[i][j];
        }
    }
    return t;
}

Matrix *multiply(Matrix *A, Matrix *B) {
    if (A->cols != B->rows) {
        fprintf(stderr, "行列の次元が一致しません\n");
        return NULL;
    }
    Matrix *C = create_matrix(A->rows, B->cols);
    for (int i = 0; i < A->rows; i++) {
        for (int j = 0; j < B->cols; j++) {
            for (int k = 0; k < A->cols; k++) {
                C->data[i][j] += A->data[i][k] * B->data[k][j];
            }
        }
    }
    return C;
}

void read_samples(const char *filename, Sample *samples, int *n) {
    FILE *fp = fopen(filename, "r");
    if (!fp) {
        fprintf(stderr, "ファイル読み込み失敗: %s\n", filename);
        exit(1);
    }

    char buffer[256];
    fgets(buffer, sizeof(buffer), fp); // ヘッダーを読み飛ばす

    *n = 0;
    while (fgets(buffer, sizeof(buffer), fp) && *n < MAX_SAMPLES) {
        if (sscanf(buffer, "%lf,%lf,%lf,%lf",
                   &samples[*n].temperature,
                   &samples[*n].tourists,
                   &samples[*n].sunny_days,
                   &samples[*n].sales) == 4) {
            (*n)++;
        }
    }
    fclose(fp);
}

Matrix *build_design_matrix(Sample *samples, int n, int features) {
    Matrix *X = create_matrix(n, features + 1);
    for (int i = 0; i < n; i++) {
        X->data[i][0] = 1.0;  // 切片用定数項
        X->data[i][1] = samples[i].temperature;
        X->data[i][2] = samples[i].tourists;
        X->data[i][3] = samples[i].sunny_days;
    }
    return X;
}

Matrix *gauss_jordan_solve(Matrix *A, Matrix *B) {
    int n = A->rows;
    int m = B->cols;

    if (n != A->cols || n != B->rows) {
        fprintf(stderr, "行列の次元不一致\n");
        return NULL;
    }

    // 増广行列 [A|B] を作成
    Matrix *aug = create_matrix(n, n + m);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            aug->data[i][j] = A->data[i][j];
        }
        for (int j = 0; j < m; j++) {
            aug->data[i][n + j] = B->data[i][j];
        }
    }

    // ガウス・ジョルダン消去法
    for (int k = 0; k < n; k++) {
        // ピボット選択
        int pivot_row = k;
        for (int i = k + 1; i < n; i++) {
            if (fabs(aug->data[i][k]) > fabs(aug->data[pivot_row][k])) {
                pivot_row = i;
            }
        }

        if (fabs(aug->data[pivot_row][k]) < 1e-12) {
            fprintf(stderr, "特異行列の可能性があります\n");
            destroy_matrix(aug);
            return NULL;
        }

        // 行交換
        if (pivot_row != k) {
            for (int j = 0; j < n + m; j++) {
                double tmp = aug->data[k][j];
                aug->data[k][j] = aug->data[pivot_row][j];
                aug->data[pivot_row][j] = tmp;
            }
        }

        // ピボット行の正規化
        double pivot_val = aug->data[k][k];
        for (int j = k; j < n + m; j++) {
            aug->data[k][j] /= pivot_val;
        }

        // 他の行の掃き出し
        for (int i = 0; i < n; i++) {
            if (i != k) {
                double factor = aug->data[i][k];
                for (int j = k; j < n + m; j++) {
                    aug->data[i][j] -= factor * aug->data[k][j];
                }
            }
        }
    }

    // 解の抽出
    Matrix *solution = create_matrix(n, m);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            solution->data[i][j] = aug->data[i][n + j];
        }
    }

    destroy_matrix(aug);
    return solution;
}

Matrix *compute_regression_coefficients(Matrix *X, Matrix *Y) {
    Matrix *XT = transpose(X);
    Matrix *XTX = multiply(XT, X);
    Matrix *XTY = multiply(XT, Y);
    Matrix *beta = gauss_jordan_solve(XTX, XTY);

    destroy_matrix(XT);
    destroy_matrix(XTX);
    destroy_matrix(XTY);

    return beta;
}

int main() {
    Sample samples[MAX_SAMPLES];
    int n = 0;

    read_samples("ice_cream2.csv", samples, &n);
    if (n < 2) {
        fprintf(stderr, "十分なデータがありません\n");
        return 1;
    }

    Matrix *X = build_design_matrix(samples, n, N_FEATURES);
    Matrix *Y = create_matrix(n, 1);
    for (int i = 0; i < n; i++) {
        Y->data[i][0] = samples[i].sales;
    }

    Matrix *coefficients = compute_regression_coefficients(X, Y);

    if (coefficients) {
        printf("回帰係数:\n");
        printf("切片: %.6f\n", coefficients->data[0][0]);
        printf("温度: %.6f\n", coefficients->data[1][0]);
        printf("観光客数: %.6f\n", coefficients->data[2][0]);
        printf("晴天日数: %.6f\n", coefficients->data[3][0]);
    }

    destroy_matrix(X);
    destroy_matrix(Y);
    destroy_matrix(coefficients);

    return 0;
}

このプログラムは、温度、観光客数、晴天日数の3つの説明変数を用いてアイスクリーム売上を予測する重回帰モデルを構築し、各係数を高精度で計算する。