---

# 0. Загрузка данных
gipers = pd.read_csv(r"Гиперспектр кукурузы.csv")
x = gipers['waveleng'].values
Y = gipers['Spectr'].values

# ==========================================
# 1. Аналитическое вычисление коэффициентов (МНК)
# ==========================================
a1 = ((x - x.mean()) * (Y - Y.mean())).mean() / ((x - x.mean()) ** 2).mean()
a0 = Y.mean() - a1 * x.mean()
print(f"Модель линейной регрессии (аналитически): Y^ = {a0} + {a1} * x")

# График аналитической регрессии
x_space = np.linspace(x.min(), x.max(), 100)
Y_pred_analyt = a0 + a1 * x_space

plt.figure()
plt.scatter(x, Y, label='Обучающая выборка')
plt.plot(x_space, Y_pred_analyt, 'r', label='МНК прогноз')
plt.title("Аналитическая парная линейная регрессия")
plt.xlabel("waveleng")
plt.ylabel("Spectr")
plt.legend()
plt.show()

# ==========================================
# 2. Оценка параметров методом градиентного спуска
# ==========================================
class SimpleRegression(object):
    def __init__(self):
        self.a0 = 0
        self.a1 = 0

    def predict(self, x):
        return self.a0 + self.a1 * x

    def MSE(self, x, Y):
        return ((Y - self.predict(x)) ** 2).mean()

    def MAE(self, x, Y):
        return abs(Y - self.predict(x)).mean()

    def MAPE(self, x, Y):
        return abs((Y - self.predict(x)) / Y).mean()

    def fit(self, x, Y, alpha=1e-9, epsylon=0.01, max_steps=5000):
        steps, errors = [], []
        step = 0
        for _ in range(max_steps):
            # Производные функции потерь по a0 и a1 из методички
            dT_a0 = -2 * sum(Y - self.predict(x))
            dT_a1 = -2 * sum((Y - self.predict(x)) * x)
            
            self.a0 -= alpha * dT_a0
            self.a1 -= alpha * dT_a1
            
            new_error = self.MSE(x, Y)
            step += 1
            steps.append(step)
            errors.append(new_error)
            
            if new_error < epsylon:
                break
        return steps, errors

regr = SimpleRegression()
steps, errors = regr.fit(x, Y, alpha=1e-9, max_steps=5000)
print(f"Модель (градиентный спуск): Y^ = {regr.a0} + {regr.a1} * x")
print("MSE после градиентного спуска:", regr.MSE(x, Y))

# График прогнозов градиентного спуска
Y_pred_gd = regr.predict(x_space)

plt.figure()
plt.scatter(x, Y, label='Обучающая выборка')
plt.plot(x_space, Y_pred_gd, 'g', label='GD прогноз')
plt.title("Парная регрессия (градиентный спуск)")
plt.xlabel("waveleng")
plt.ylabel("Spectr")
plt.legend()
plt.show()

# График кривой обучения (убывание MSE)
plt.figure()
plt.plot(steps, errors)
plt.title("Кривая обучения (изменение MSE)")
plt.xlabel("Шаги")
plt.ylabel("MSE")
plt.show()

# ==========================================
# 3. Сравнение качества моделей и прогнозов
# ==========================================
plt.figure()
plt.scatter(x, Y, color='blue', alpha=0.3, label='Обучающая выборка')
plt.plot(x_space, Y_pred_analyt, 'r', linewidth=2, label='Аналитическое решение')
plt.plot(x_space, Y_pred_gd, 'g--', linewidth=2, label='Градиентный спуск')
plt.title("Сравнение аналитической модели и градиентного спуска")
plt.xlabel("waveleng")
plt.ylabel("Spectr")
plt.legend()
plt.show()

mse_analyt = ((Y - (a0 + a1 * x)) ** 2).mean()
mse_gd = regr.MSE(x, Y)
print(f"MSE аналитической модели: {mse_analyt:.4f}")
print(f"MSE градиентного спуска: {mse_gd:.4f}")