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

import math


DATA_RAW = {
    18: [0.60, 0.64, 0.69, 0.52, 0.65, 0.52, 0.58, 0.63, 0.69, 0.48],
    24: [0.61, 0.72, 0.72, 0.58, 0.46, 0.52, 0.72, 0.59, 0.71, 0.79],
    6:  [0.58, 0.59, 0.66, 0.46, 0.47, 0.56, 0.70, 0.60, 0.60, 0.64],
    12: [0.70, 0.57, 0.67, 0.87, 0.52, 0.66, None, 0.59, 0.65, 0.66],
}

HOURS_ORDER = [18, 24, 6, 12]


def fmt(x, width=6, prec=2):
    if x is None:
        return f"{'—':>{width}}"
    return f"{x:>{width}.{prec}f}"


def print_input_table(data_raw, hours_order):
    max_cols = max(len(v) for v in data_raw.values())
    print("\n--- Вхідні дані ---")
    header = "Години | " + " ".join([f"{j:>6d}" for j in range(1, max_cols + 1)])
    print(header)
    print("-" * len(header))
    for h in hours_order:
        row = data_raw[h] + [None] * (max_cols - len(data_raw[h]))
        print(f"{h:>6d} | " + " ".join(fmt(v) for v in row))


def clean_groups(data_raw, hours_order):
    groups = {}
    stats = {}
    for h in hours_order:
        xs = [v for v in data_raw[h] if v is not None]
        groups[h] = xs
        n_i = len(xs)
        T_i = sum(xs)
        mean_i = T_i / n_i if n_i > 0 else float("nan")
        stats[h] = (n_i, T_i, mean_i)
    return groups, stats


def anova_one_way(groups, stats, hours_order):
    all_x = []
    for h in hours_order:
        all_x.extend(groups[h])

    N = len(all_x)
    T = sum(all_x)
    sum_x2 = sum(x * x for x in all_x)

    C = (T * T) / N

    sum_Ti2_over_ni = 0.0
    for h in hours_order:
        n_i, T_i, _ = stats[h]
        sum_Ti2_over_ni += (T_i * T_i) / n_i

    SS_total = sum_x2 - C
    SS_factor = sum_Ti2_over_ni - C
    SS_error = sum_x2 - sum_Ti2_over_ni

    a = len(hours_order)

    df_total = N - 1
    df_factor = a - 1
    df_error = N - a

    MS_factor = SS_factor / df_factor
    MS_error = SS_error / df_error

    F_data = MS_factor / MS_error

    return {
        "N": N,
        "T": T,
        "sum_x2": sum_x2,
        "C": C,
        "a": a,
        "sum_Ti2_over_ni": sum_Ti2_over_ni,
        "SS_total": SS_total,
        "SS_factor": SS_factor,
        "SS_error": SS_error,
        "df_total": df_total,
        "df_factor": df_factor,
        "df_error": df_error,
        "MS_factor": MS_factor,
        "MS_error": MS_error,
        "F_data": F_data,
    }


def print_table1(data_raw, stats, hours_order):
    max_cols = max(len(v) for v in data_raw.values())
    print("\n--- Таблиця 1. Групи за фактором A + спостереження + суми/середні ---")
    header = "Година | " + " ".join([f"{j:>6d}" for j in range(1, max_cols + 1)]) + " |   T_i  | n_i |  x̄_i"
    print(header)
    print("-" * len(header))

    for h in hours_order:
        row = data_raw[h] + [None] * (max_cols - len(data_raw[h]))
        n_i, T_i, mean_i = stats[h]
        print(
            f"{h:>6d} | "
            + " ".join(fmt(v) for v in row)
            + f" | {T_i:>6.3f} | {n_i:>3d} | {mean_i:>5.4f}"
        )


def print_anova_results(res, alpha_levels=(0.05, 0.01)):
    df1 = res["df_factor"]
    df2 = res["df_error"]

    fcrit = {alpha: f_isf(alpha, df1, df2) for alpha in alpha_levels}

    print("\n--- Проміжні підсумки ---")
    print(f"Загальна сума T = {res['T']:.4f}")
    print(f"Загальна кількість N = {res['N']}")
    print(f"Сума квадратів усіх x (Σx^2) = {res['sum_x2']:.4f}")
    print(f"Поправковий коефіцієнт C = T^2 / N = {res['C']:.4f}")
    print(f"Σ(T_i^2 / n_i) = {res['sum_Ti2_over_ni']:.4f}")

    print("\n--- Таблиця 2. Підсумкова таблиця дисперсійного аналізу ---")
    print("Джерело варіації        |      SS |  df |      MS |    F_data |  F_0.05 |  F_0.01")
    print("-" * 86)

    ss_a = res["SS_factor"]
    ss_e = res["SS_error"]
    ss_t = res["SS_total"]

    ms_a = res["MS_factor"]
    ms_e = res["MS_error"]
    f_data = res["F_data"]

    f05 = fcrit.get(0.05, float("nan"))
    f01 = fcrit.get(0.01, float("nan"))

    print(f"Між групами (Фактор A)  | {ss_a:>7.4f} | {df1:>3d} | {ms_a:>7.4f} | {f_data:>8.4f} | {f05:>7.4f} | {f01:>7.4f}")
    print(f"Всередині груп (Помилка)| {ss_e:>7.4f} | {df2:>3d} | {ms_e:>7.4f} | {'':>8} | {'':>7} | {'':>7}")
    print(f"Загальна                | {ss_t:>7.4f} | {res['df_total']:>3d} | {'':>7} | {'':>8} | {'':>7} | {'':>7}")

    print("\n--- Критерій Фішера ---")
    print("F_data = MS_фактор / MS_помилка")
    print(f"F_data = {f_data:.6f}")
    print(f"F_crit(α=0.05; df1={df1}, df2={df2}) = {f05:.6f}")
    print(f"F_crit(α=0.01; df1={df1}, df2={df2}) = {f01:.6f}")

    print("\n--- Висновок ---")
    for alpha in alpha_levels:
        fc = fcrit[alpha]
        if f_data < fc:
            print(f"α = {alpha:.2f}: F_data < F_crit ⇒ H0 приймається (вплив фактору A статистично НЕ значущий).")
        else:
            print(f"α = {alpha:.2f}: F_data ≥ F_crit ⇒ H0 відхиляється (вплив фактору A статистично значущий).")


def _betacf(a, b, x, max_iter=200, eps=3e-14):
    am = 1.0
    bm = 1.0
    az = 1.0
    qab = a + b
    qap = a + 1.0
    qam = a - 1.0

    bz = 1.0 - qab * x / qap
    if abs(bz) < 1e-30:
        bz = 1e-30

    aold = 0.0
    for m in range(1, max_iter + 1):
        em = float(m)
        tem = em + em

        d = em * (b - em) * x / ((qam + tem) * (a + tem))
        ap = az + d * am
        bp = bz + d * bm

        d = -(a + em) * (qab + em) * x / ((a + tem) * (qap + tem))
        app = ap + d * az
        bpp = bp + d * bz

        aold = az
        am = ap / bpp
        bm = bp / bpp
        az = app / bpp
        bz = 1.0

        if abs(az - aold) < eps * abs(az):
            return az

    return az


def betai(a, b, x):
    if x <= 0.0:
        return 0.0
    if x >= 1.0:
        return 1.0

    ln_beta = math.lgamma(a) + math.lgamma(b) - math.lgamma(a + b)
    bt = math.exp(a * math.log(x) + b * math.log(1.0 - x) - ln_beta)

    if x < (a + 1.0) / (a + b + 2.0):
        return bt * _betacf(a, b, x) / a
    return 1.0 - bt * _betacf(b, a, 1.0 - x) / b


def f_cdf(x, d1, d2):
    if x <= 0.0:
        return 0.0
    a = d1 / 2.0
    b = d2 / 2.0
    y = (d1 * x) / (d1 * x + d2)
    return betai(a, b, y)


def f_isf(alpha, d1, d2):
    target_cdf = 1.0 - alpha

    lo = 0.0
    hi = 1.0
    while f_cdf(hi, d1, d2) < target_cdf:
        hi *= 2.0
        if hi > 1e9:
            break

    for _ in range(200):
        mid = (lo + hi) / 2.0
        if f_cdf(mid, d1, d2) < target_cdf:
            lo = mid
        else:
            hi = mid
        if hi - lo < 1e-12 * max(1.0, hi):
            break

    return (lo + hi) / 2.0


def main():
    print("=== ПР №9–10. Дисперсійний аналіз даних. ВАРІАНТ 2 ===")
    print("Фактор A: година доби. Об'єкт: Syringa Emodi (вміст каротиноїдів).")

    print_input_table(DATA_RAW, HOURS_ORDER)

    groups, stats = clean_groups(DATA_RAW, HOURS_ORDER)
    print_table1(DATA_RAW, stats, HOURS_ORDER)

    res = anova_one_way(groups, stats, HOURS_ORDER)
    print_anova_results(res, alpha_levels=(0.05, 0.01))

if __name__ == "__main__":
    main()