lab6.py

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme()  # чтобы графики были красивее

# 1) Исходные данные (вариант 17)
x = np.array([1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0])
y = np.array([2.61, 1.62, 1.17, 0.75, 0.30, 0.75, 1.03, 0.81, 0.57])
n = len(x)

# Крайние и «средние» точки
x0, xn = x[0], x[-1]
y0, yn = y[0], y[-1]

x_a = (x0 + xn)/2
x_g = np.sqrt(x0 * xn)
x_h = 2/(1/x0 + 1/xn)

y_a = (y0 + yn)/2
y_g = np.sqrt(y0 * yn)
y_h = 2/(1/y0 + 1/yn)

# Общая функция для решения СЛАУ нормальных уравнений:
#   [ ΣU^2    ΣU   ] [A] = [ Σ(U·V) ]
#   [ ΣU     N     ] [B]   [ Σ V    ]
def solve_normal(U, V):
    S_UU = np.sum(U*U)
    S_U  = np.sum(U)
    S_V  = np.sum(V)
    S_UV = np.sum(U*V)
    # матрица и правая часть
    M = np.array([[S_UU, S_U],
                  [S_U,  n   ]])
    rhs = np.array([S_UV, S_V])
    A, B = np.linalg.solve(M, rhs)
    return A, B

# Словарь для хранения результатов
models = {}

# --- z1(x) = a*x + b -----------------------------------
# U = x, V = y
A, B = solve_normal(x, y)
# здесь A=a, B=b
models['z1'] = {
    'func': lambda xx,a,b: a*xx + b,
    'a': A, 'b': B,
    'x_mean': x_a, 'y_mean': y_a
}

# --- z2(x) = a * x^b -----------------------------------
# ln y = ln a + b ln x  => U = ln x, V = ln y, A=b, B=ln a
U = np.log(x)
V = np.log(y)
b_lin, ln_a = solve_normal(U, V)
models['z2'] = {
    'func': lambda xx,a,b: a * xx**b,
    'a': np.exp(ln_a), 'b': b_lin,
    'x_mean': x_g, 'y_mean': y_g
}

# --- z3(x) = a * exp(b x) -----------------------------
# ln y = ln a + b x   => U = x, V=ln y
U = x
V = np.log(y)
b_lin, ln_a = solve_normal(U, V)
models['z3'] = {
    'func': lambda xx,a,b: a * np.exp(b*xx),
    'a': np.exp(ln_a), 'b': b_lin,
    'x_mean': x_a, 'y_mean': y_g
}

# --- z4(x) = a ln x + b -------------------------------
# U = ln x, V = y   => A=a, B=b
U = np.log(x)
V = y
A, B = solve_normal(U, V)
models['z4'] = {
    'func': lambda xx,a,b: a*np.log(xx) + b,
    'a': A, 'b': B,
    'x_mean': x_g, 'y_mean': y_a
}

# --- z5(x) = a/x + b -----------------------------------
# U = 1/x, V = y   => A=a, B=b
U = 1/x
V = y
A, B = solve_normal(U, V)
models['z5'] = {
    'func': lambda xx,a,b: a/xx + b,
    'a': A, 'b': B,
    'x_mean': x_h, 'y_mean': y_a
}

# --- z6(x) = 1/(a x + b) ------------------------------
# 1/y = a x + b    => U = x, V = 1/y
U = x
V = 1/y
A, B = solve_normal(U, V)
# теперь 1/y ≈ A*x + B  => a=A, b=B
models['z6'] = {
    'func': lambda xx,a,b: 1/(a*xx + b),
    'a': A, 'b': B,
    'x_mean': x_a, 'y_mean': y_h
}

# --- z7(x) = x/(a x + b) ------------------------------
# 1/y = a + b*(1/x)  => U = 1/x, V=1/y
U = 1/x
V = 1/y
b_lin, a_int = solve_normal(U, V)
# V ≈ b_lin*U + a_int  => a=a_int, b=b_lin
models['z7'] = {
    'func': lambda xx,a,b: xx/(a*xx + b),
    'a': a_int, 'b': b_lin,
    'x_mean': x_h, 'y_mean': y_h
}

# --- z8(x) = a * exp(b/x) -----------------------------
# ln y = ln a + b*(1/x)  => U = 1/x, V = ln y
U = 1/x
V = np.log(y)
b_lin, ln_a = solve_normal(U, V)
models['z8'] = {
    'func': lambda xx,a,b: a * np.exp(b/xx),
    'a': np.exp(ln_a), 'b': b_lin,
    'x_mean': x_h, 'y_mean': y_g
}

# --- z9(x) = 1/(a ln x + b) ---------------------------
# 1/y = a ln x + b  => U = ln x, V = 1/y
U = np.log(x)
V = 1/y
A, B = solve_normal(U, V)
# 1/y ≈ A*ln x + B  => a=A, b=B
models['z9'] = {
    'func': lambda xx,a,b: 1/(a*np.log(xx) + b),
    'a': A, 'b': B,
    'x_mean': x_g, 'y_mean': y_h
}


# 4) Вычисляем локальные ошибки δ_i
errors = {}
for name,m in models.items():
    z_est = m['func'](m['x_mean'], m['a'], m['b'])
    errors[name] = abs(z_est - m['y_mean'])

# 5) Выбираем лучшую модель
best = min(errors, key=errors.get)
a_best, b_best = models[best]['a'], models[best]['b']
print(f"Лучшая модель: {best},  a={a_best:.4f}, b={b_best:.4f}")

# Расчет среднеквадратичного отклонения (сигма)
y_pred_best = models[best]['func'](x, a_best, b_best)
sum_sq_err = np.sum((y - y_pred_best)**2)
sigma = np.sqrt(sum_sq_err / n)
print(f"Среднеквадратичное отклонение (σ): {sigma:.4f}")

# 6) Строим график
x_plot = np.linspace(x.min(), x.max(), 200)
y_plot = models[best]['func'](x_plot, a_best, b_best)

plt.figure(figsize=(8,5))
plt.scatter(x, y, label="Исходные точки", s=50)
plt.plot(x_plot, y_plot, lw=2, label=f"{best}(x)")
plt.xlabel("x")
plt.ylabel("z(x)")
plt.title(f"Аппроксимация моделью {best}")
plt.legend()
plt.tight_layout()
plt.show()