import numpy as np
import pandas as pd
import math

data = {
    'x1': [32.8, 29.9, 30.12, 31.15, 25.25, 32.8, 29.9, 24.12, 20.15, 33.25], 
    'x2': [0.92, 0.46, 0.73, 0.64, 0.54, 0.74, 0.68, 0.46, 0.36, 0.95],       
    'x3': [20.0, 16.2, 14.2, 10.0, 24.1, 13.2, 18.0, 15.1, 15.2, 75.5]        
}

df = pd.DataFrame(data)
n = len(df)
m = len(df.columns)

print("=== Вхідні дані (Варіант 2) ===")
print(df)
print("-" * 50)

def manual_pearson_corr(x, y):
    x_arr = np.array(x)
    y_arr = np.array(y)
    
    x_mean = np.mean(x_arr)
    y_mean = np.mean(y_arr)
    
    numerator = np.sum((x_arr - x_mean) * (y_arr - y_mean))
    
    sum_sq_diff_x = np.sum((x_arr - x_mean) ** 2)
    sum_sq_diff_y = np.sum((y_arr - y_mean) ** 2)
    denominator = np.sqrt(sum_sq_diff_x * sum_sq_diff_y)
    
    return numerator / denominator

cols = df.columns
R_matrix = np.zeros((m, m))

for i in range(m):
    for j in range(m):
        if i == j:
            R_matrix[i, j] = 1.0
        else:
            r = manual_pearson_corr(df[cols[i]], df[cols[j]])
            R_matrix[i, j] = r

corr_matrix = pd.DataFrame(R_matrix, index=cols, columns=cols)

det_R = np.linalg.det(corr_matrix.values)
chi2_observed = - (n - 1 - (2 * m + 5) / 6) * np.log(det_R)

# df = m(m-1)/2 = 3(2)/2 = 3
# alpha = 0.05
# Таблиця: chi2(0.95, 3) ≈ 7.815
CHI2_CRITICAL_TABLE = 7.815 

print("\n>>> КРОК 1: Перевірка всього масиву (Критерій Пірсона)")
print(f"Кореляційна матриця R:\n{corr_matrix.round(4)}")
print(f"Визначник матриці det(R) = {det_R:.6f}")
print(f"Хі-квадрат спостережуване = {chi2_observed:.4f}")
print(f"Хі-квадрат критичне (табличне, df=3) = {CHI2_CRITICAL_TABLE}")

if chi2_observed > CHI2_CRITICAL_TABLE:
    print("ВИСНОВОК: Мультиколінеарність ІСНУЄ.")
    
    inv_matrix_values = np.linalg.inv(corr_matrix.values)
    inv_corr_matrix = pd.DataFrame(inv_matrix_values, index=cols, columns=cols)
    
    print("\n>>> КРОК 2: Оцінка кожної змінної (F-критерій)")
    print(f"Обернена матриця C:\n{inv_corr_matrix.round(4)}")
    
    # df1 = m - 1 = 2
    # df2 = n - m = 7
    # Таблиця: F(0.95, 2, 7) ≈ 4.737
    F_CRITICAL_TABLE = 4.737
    
    f_results = []
    collinear_vars = []
    
    for col in cols:
        c_kk = inv_corr_matrix.loc[col, col]
        F_obs = (c_kk - 1) * (n - m) / (m - 1)
        
        is_bad = F_obs > F_CRITICAL_TABLE
        if is_bad: collinear_vars.append(col)
        
        f_results.append({
            "Змінна": col,
            "c_kk": round(c_kk, 4),
            "F_спост": round(F_obs, 4),
            "F_крит": F_CRITICAL_TABLE,
            "Результат": "Мультиколінеарна" if is_bad else "Норма"
        })
    
    print(pd.DataFrame(f_results).to_string(index=False))
    
    print("\n>>> КРОК 3: Виявлення пар (t-критерій)")
    
    # df = n - m - 1 = 6
    # alpha = 0.05 (двостороння), тому шукаємо t(0.975, 6)
    # Таблиця: t(0.975, 6) ≈ 2.447
    T_CRITICAL_TABLE = 2.447
    
    pairs_results = []
    conflict_pairs = []
    
    for i in range(m):
        for j in range(i + 1, m):
            var1 = cols[i]
            var2 = cols[j]
            
            c_ij = inv_corr_matrix.loc[var1, var2]
            c_ii = inv_corr_matrix.loc[var1, var1]
            c_jj = inv_corr_matrix.loc[var2, var2]
            
            r_partial = -c_ij / math.sqrt(c_ii * c_jj)
            
            if abs(r_partial) >= 1.0:
                t_obs = float('inf')
            else:
                t_obs = (abs(r_partial) * math.sqrt(n - m - 1)) / math.sqrt(1 - r_partial**2)
            
            is_significant = t_obs > T_CRITICAL_TABLE
            
            if is_significant:
                conflict_pairs.append((var1, var2, r_partial))
            
            pairs_results.append({
                "Пара": f"{var1} - {var2}",
                "r_частковий": round(r_partial, 4),
                "t_спост": round(t_obs, 4),
                "t_крит": T_CRITICAL_TABLE,
                "Зв'язок": "Є мультикол." if is_significant else "Немає"
            })
            
    print(pd.DataFrame(pairs_results).to_string(index=False))
    
    print("\n" + "="*50)
    print("РІШЕННЯ ЩОДО КОРИГУВАННЯ МОДЕЛІ")
    print("="*50)
    
    if conflict_pairs:
        print("Виявлено проблемні пари:")
        variables_to_drop = set()
        for v1, v2, r in conflict_pairs:
            print(f"- {v1} та {v2} (часткова кореляція r = {r:.4f})")
            
            if 'x1' in (v1, v2) and 'x2' in (v1, v2):
                print(f"  -> Залишаємо x1 (Продуктивність праці), прибираємо x2 (Фондомісткість).")
                if v1 == 'x2': variables_to_drop.add(v1)
                else: variables_to_drop.add(v2)
            elif 'x3' in (v1, v2):
                print(f"  -> Прибираємо x3 (Плинність кадрів) як менш стабільний фактор.")
                variables_to_drop.add('x3')
        
        if variables_to_drop:
            print(f"\nОстаточне рішення: виключаємо з моделі змінну(і): {', '.join(variables_to_drop)}")
    else:
        print("Значущої колінеарності між окремими парами за t-критерієм не виявлено. Змінні залишаємо.")
else:
    print("ВИСНОВОК: Мультиколінеарність відсутня. Всі змінні залишаємо в моделі.")