Мірошниченко_ГЙМ_ПР№5-6

from __future__ import annotations
import math
from dataclasses import dataclass
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


def generate_data_var2(n: int = 100, seed: int = 42) -> pd.DataFrame:

    rng = np.random.default_rng(seed)

    x = rng.uniform(3500, 11000, size=n)
    eps = rng.normal(0.0, 1.0, size=n)
    y = -1000.0 + 0.5 * x + 300.0 * eps

    df = pd.DataFrame({"x": x, "y": y})
    return df


@dataclass
class CorrelationResult:
    r_xy: float
    r_np: float
    sums: dict


def compute_correlation_and_sums(df: pd.DataFrame) -> CorrelationResult:

    x = df["x"].to_numpy()
    y = df["y"].to_numpy()
    n = len(df)

    Sx = x.sum()
    Sy = y.sum()
    Sx2 = (x ** 2).sum()
    Sy2 = (y ** 2).sum()
    Sxy = (x * y).sum()

    numerator = n * Sxy - Sx * Sy
    denominator = math.sqrt((n * Sx2 - Sx ** 2) * (n * Sy2 - Sy ** 2))
    r_xy = numerator / denominator

    r_np = float(np.corrcoef(x, y)[0, 1])

    Sx3 = (x ** 3).sum()
    Sx4 = (x ** 4).sum()
    Sx5 = (x ** 5).sum()
    Sx6 = (x ** 6).sum()
    Sx2y = ((x ** 2) * y).sum()
    Sx3y = ((x ** 3) * y).sum()

    sums = {
        "n": float(n),
        "Sx": float(Sx),
        "Sx2": float(Sx2),
        "Sx3": float(Sx3),
        "Sx4": float(Sx4),
        "Sx5": float(Sx5),
        "Sx6": float(Sx6),
        "Sy": float(Sy),
        "Sxy": float(Sxy),
        "Sx2y": float(Sx2y),
        "Sx3y": float(Sx3y),
    }

    return CorrelationResult(r_xy=r_xy, r_np=r_np, sums=sums)


def interpret_correlation(r: float) -> str:

    r_abs = abs(r)
    if r_abs < 0.3:
        return "майжу відстутній"
    elif r_abs < 0.5:
        return "слабкий"
    elif r_abs < 0.7:
        return "помірний"
    elif r_abs < 1.0:
        return "сильний"
    else:
        return "функціональний"




@dataclass
class RegressionModel:
    degree: int
    coeffs: np.ndarray
    coeffs_check: np.ndarray
    rss: float
    sigma2: float

    def predict(self, x: np.ndarray) -> np.ndarray:
        powers = np.arange(self.degree + 1)
        return sum(self.coeffs[j] * (x ** powers[j]) for j in range(self.degree + 1))


def solve_regression_models(df: pd.DataFrame, sums: dict) -> list[RegressionModel]:

    x = df["x"].to_numpy()
    y = df["y"].to_numpy()
    n = int(sums["n"])

    Sx = sums["Sx"]
    Sx2 = sums["Sx2"]
    Sx3 = sums["Sx3"]
    Sx4 = sums["Sx4"]
    Sx5 = sums["Sx5"]
    Sx6 = sums["Sx6"]
    Sy = sums["Sy"]
    Sxy = sums["Sxy"]
    Sx2y = sums["Sx2y"]
    Sx3y = sums["Sx3y"]

    models: list[RegressionModel] = []

    A1 = np.array([[n, Sx], [Sx, Sx2]], dtype=float)
    b1 = np.array([Sy, Sxy], dtype=float)
    coeffs1 = np.linalg.solve(A1, b1)

    X1 = np.vstack([np.ones_like(x), x]).T
    coeffs1_check, *_ = np.linalg.lstsq(X1, y, rcond=None)
    resid1 = y - (coeffs1[0] + coeffs1[1] * x)
    rss1 = float(np.sum(resid1 ** 2))
    sigma2_1 = rss1 / (n - 2)

    models.append(RegressionModel(1, coeffs1, coeffs1_check, rss1, sigma2_1))

    A2 = np.array(
        [
            [n, Sx, Sx2],
            [Sx, Sx2, Sx3],
            [Sx2, Sx3, Sx4],
        ],
        dtype=float,
    )
    b2 = np.array([Sy, Sxy, Sx2y], dtype=float)
    coeffs2 = np.linalg.solve(A2, b2)

    X2 = np.vstack([np.ones_like(x), x, x ** 2]).T
    coeffs2_check, *_ = np.linalg.lstsq(X2, y, rcond=None)
    y_hat2 = coeffs2[0] + coeffs2[1] * x + coeffs2[2] * x ** 2
    resid2 = y - y_hat2
    rss2 = float(np.sum(resid2 ** 2))
    sigma2_2 = rss2 / (n - 3)

    models.append(RegressionModel(2, coeffs2, coeffs2_check, rss2, sigma2_2))

    A3 = np.array(
        [
            [n, Sx, Sx2, Sx3],
            [Sx, Sx2, Sx3, Sx4],
            [Sx2, Sx3, Sx4, Sx5],
            [Sx3, Sx4, Sx5, Sx6],
        ],
        dtype=float,
    )
    b3 = np.array([Sy, Sxy, Sx2y, Sx3y], dtype=float)
    coeffs3 = np.linalg.solve(A3, b3)

    X3 = np.vstack([np.ones_like(x), x, x ** 2, x ** 3]).T
    coeffs3_check, *_ = np.linalg.lstsq(X3, y, rcond=None)
    y_hat3 = coeffs3[0] + coeffs3[1] * x + coeffs3[2] * x ** 2 + coeffs3[3] * x ** 3
    resid3 = y - y_hat3
    rss3 = float(np.sum(resid3 ** 2))
    sigma2_3 = rss3 / (n - 4)

    models.append(RegressionModel(3, coeffs3, coeffs3_check, rss3, sigma2_3))

    return models




@dataclass
class EmpiricalRegression:
    intervals: pd.DataFrame
    x_mid: np.ndarray
    y_mean: np.ndarray


def build_empirical_regression(df: pd.DataFrame) -> EmpiricalRegression:

    x = df["x"].to_numpy()
    y = df["y"].to_numpy()
    n = len(df)

    m = int(round(1 + 3.322 * math.log10(n)))
    m = max(m, 4)

    x_min, x_max = float(x.min()), float(x.max())
    h = (x_max - x_min) / m

    bounds = [x_min + i * h for i in range(m + 1)]

    rows = []
    x_mid_list = []
    y_mean_list = []

    for i in range(m):
        a_i = bounds[i]
        b_i = bounds[i + 1] if i < m - 1 else bounds[i + 1] + 1e-6

        mask = (x >= a_i) & (x < b_i)
        x_bin = x[mask]
        y_bin = y[mask]

        if len(x_bin) == 0:
            continue

        mid = 0.5 * (a_i + b_i)
        y_mean = float(y_bin.mean())

        rows.append(
            {
                "Інтервал X": f"[{a_i:8.2f}; {b_i:8.2f})",
                "Середина X": mid,
                "N": int(len(x_bin)),
                "Середнє Y (умовне)": y_mean,
            }
        )
        x_mid_list.append(mid)
        y_mean_list.append(y_mean)

    table = pd.DataFrame(rows)
    return EmpiricalRegression(table, np.array(x_mid_list), np.array(y_mean_list))




def print_header(title: str) -> None:
    print("\n" + "=" * 3 + f" {title} " + "=" * 3)


def main():
    print_header("Практична робота №5, 6")
    print("Побудова кореляційного поля і рівняння регресії")
    print("Варіант 2: Дохід (x), Заощадження (y)\n")

    df = generate_data_var2(n=100, seed=42)

    print("Перші 10 спостережень (округлені значення):")
    print(df.round(0).head(10))

    print_header("Коефіцієнт лінійної кореляції")
    corr_res = compute_correlation_and_sums(df)
    print(f"r_xy (за формулою сум) = {corr_res.r_xy: .6f}")
    print(f"r_xy (np.corrcoef)      = {corr_res.r_np: .6f}")
    print(f"Тіснота зв'язку          = {interpret_correlation(corr_res.r_xy)}")

    print("\nПроміжні суми (для нормальних рівнянь):")
    for key in ["n", "Sx", "Sx2", "Sx3", "Sx4", "Sx5", "Sx6", "Sy", "Sxy", "Sx2y", "Sx3y"]:
        print(f"  {key:>4} = {corr_res.sums[key]:.4f}")

    print_header("Коефіцієнти рівнянь регресії (матричний метод)")
    models = solve_regression_models(df, corr_res.sums)

    for model in models:
        a = model.coeffs
        degree = model.degree

        if degree == 1:
            print("Лінійна регресія (1-й порядок):")
            print(f"  y = {a[0]: .6f} + {a[1]: .6f} * x")
        elif degree == 2:
            print("Квадратична регресія (2-й порядок):")
            print(f"  y = {a[0]: .6f} + {a[1]: .6f} * x + {a[2]: .10f} * x^2")
        elif degree == 3:
            print("Кубічна регресія (3-й порядок):")
            print(
                "  y = "
                f"{a[0]: .6f} + {a[1]: .6f} * x + {a[2]: .10f} * x^2 + {a[3]: .10e} * x^3"
            )
        print()

    print_header("Перевірка коефіцієнтів")
    for model in models:
        print(f"Модель ступеня {model.degree}:")
        print(f"  Матричний метод:  {model.coeffs}")
        print(f"  lstsq:            {model.coeffs_check}")
        print()

    print_header("Залишкова дисперсія для моделей")
    print(f"{'Порядок':>8}  {'RSS':>15}  {'σ^2_зал':>15}")
    for model in models:
        print(
            f"{model.degree:8d}  "
            f"{model.rss:15.4f}  "
            f"{model.sigma2:15.4f}"
        )

    best_model = min(models, key=lambda m: m.sigma2)
    if best_model.degree == 1:
        name_best = "лінійна"
    elif best_model.degree == 2:
        name_best = "квадратична"
    else:
        name_best = "кубічна"
    print(f"\nНайкраща модель за критерієм мінімальної σ^2_зал: {name_best} (ступінь {best_model.degree}).")

    print_header("Емпірична регресія (умовне середнє)")
    emp = build_empirical_regression(df)
    print(emp.intervals.to_string(index=False, justify="center", col_space=12))

    x = df["x"].to_numpy()
    y = df["y"].to_numpy()

    x_grid = np.linspace(x.min(), x.max(), 300)

    y1_grid = models[0].predict(x_grid)
    y2_grid = models[1].predict(x_grid)
    y3_grid = models[2].predict(x_grid)

    plt.figure(figsize=(10, 6))

    plt.scatter(x, y, s=20, alpha=0.7, label="Спостереження (x, y)")

    plt.plot(x_grid, y1_grid, label="Лінійна регресія")
    plt.plot(x_grid, y2_grid, linestyle="--", label="Квадратична регресія")
    plt.plot(x_grid, y3_grid, linestyle=":", label="Кубічна регресія")

    plt.plot(
        emp.x_mid,
        emp.y_mean,
        "o-",
        linewidth=2,
        markersize=5,
        label="Емпірична регресія (умовне середнє)",
    )

    plt.xlabel("x – Дохід, грн")
    plt.ylabel("y – Заощадження, грн")
    plt.title("Кореляційне поле з лініями регресії (варіант 2)")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()