import numpy as np
import pandas as pd

def main():
    Y = np.array([7.60, 3.50, 1.40, 3.90, 6.40, 7.30, 7.60, 8.30, 4.40, 5.60, 7.50, 3.20])
    X1 = np.array([1.40, 0.70, 0.40, 0.80, 1.10, 1.40, 1.50, 0.90, 0.80, 0.70, 1.40, 0.60])
    X2 = np.array([3.80, 1.70, 0.70, 1.80, 3.20, 3.60, 3.90, 4.10, 2.20, 2.90, 3.20, 1.70])
    
    n = len(Y)

    print("--- Вхідні дані ---")
    print(pd.DataFrame({"Y": Y, "X1": X1, "X2": X2}, index=range(1, n+1)))
    print("-" * 40)

    sum_x1 = np.sum(X1)
    sum_x2 = np.sum(X2)
    sum_y  = np.sum(Y)
    
    sum_x1_sq = np.sum(X1**2)
    sum_x2_sq = np.sum(X2**2)
    sum_x1_x2 = np.sum(X1 * X2)
    
    sum_y_x1 = np.sum(Y * X1)
    sum_y_x2 = np.sum(Y * X2)

    A = np.array([
        [n,       sum_x1,    sum_x2],
        [sum_x1,  sum_x1_sq, sum_x1_x2],
        [sum_x2,  sum_x1_x2, sum_x2_sq]
    ])

    B = np.array([sum_y, sum_y_x1, sum_y_x2])

    coeffs = np.linalg.solve(A, B)
    a0, a1, a2 = coeffs

    print(f"\n--- 1. Оцінки параметрів моделі (МНК) ---")
    print(f"Рівняння регресії: Y = {a0:.4f} + {a1:.4f}*X1 + {a2:.4f}*X2")
    print(f"a0 = {a0:.4f}")
    print(f"a1 = {a1:.4f}")
    print(f"a2 = {a2:.4f}")

    def manual_pearson(x, y):
        mean_x = np.mean(x)
        mean_y = np.mean(y)
        numerator = np.sum((x - mean_x) * (y - mean_y))
        denominator = np.sqrt(np.sum((x - mean_x)**2) * np.sum((y - mean_y)**2))
        return numerator / denominator

    r_y_x1 = manual_pearson(Y, X1)
    r_y_x2 = manual_pearson(Y, X2)
    r_x1_x2 = manual_pearson(X1, X2)

    print(f"\n--- 2. Кореляційний аналіз ---")
    print(f"Парні коефіцієнти кореляції:")
    print(f"r(Y, X1)  = {r_y_x1:.4f}")
    print(f"r(Y, X2)  = {r_y_x2:.4f}")
    print(f"r(X1, X2) = {r_x1_x2:.4f}")

    numerator_1 = r_y_x1 - r_y_x2 * r_x1_x2
    denominator_1 = np.sqrt((1 - r_y_x2**2) * (1 - r_x1_x2**2))
    r_y_x1_part = numerator_1 / denominator_1

    numerator_2 = r_y_x2 - r_y_x1 * r_x1_x2
    denominator_2 = np.sqrt((1 - r_y_x1**2) * (1 - r_x1_x2**2))
    r_y_x2_part = numerator_2 / denominator_2

    print(f"\nЧасткові коефіцієнти кореляції:")
    print(f"r(Y, X1 | X2) = {r_y_x1_part:.4f}")
    print(f"r(Y, X2 | X1) = {r_y_x2_part:.4f}")

    R2_numerator = r_y_x1**2 + r_y_x2**2 - 2 * r_y_x1 * r_y_x2 * r_x1_x2
    R2_denominator = 1 - r_x1_x2**2
    R2 = R2_numerator / R2_denominator
    
    R = np.sqrt(R2)

    print(f"\n--- 3. Показники якості моделі ---")
    print(f"Коефіцієнт множинної детермінації (R^2): {R2:.4f}")
    print(f"Коефіцієнт множинної кореляції (R): {R:.4f}")

    print(f"\n--- 4. Аналіз результатів ---")

    def analyze(simple_r, partial_r, var_name, control_var):
        print(f"\nАналіз зв'язку Y <-> {var_name} (при фіксованому {control_var}):")
        print(f"  Парний r: {simple_r:.4f}")
        print(f"  Частковий r: {partial_r:.4f}")
        
        diff = abs(partial_r) - abs(simple_r)
        
        if diff > 0.01:
            print(f"  > Часткова кореляція ЗРОСЛА (по модулю).")
            print(f"  > Висновок: Фактор '{control_var}' приховував (маскував) істинну взаємозалежність.")
        elif diff < -0.01:
            print(f"  > Часткова кореляція ЗМЕНШИЛАСЯ (по модулю).")
            print(f"  > Висновок: Взаємозалежність значною мірою мала місце завдяки фактору '{control_var}'.")
        else:
            print(f"  > Кореляція майже не змінилася.")

    analyze(r_y_x1, r_y_x1_part, "X1 (Площа)", "X2 (Потік)")
    analyze(r_y_x2, r_y_x2_part, "X2 (Потік)", "X1 (Площа)")

if __name__ == "__main__":
    main()