1. Исследование регрессионных моделей

from numpy import *
from scipy.stats import norm
from scipy.stats.distributions import chi2
import matplotlib.pyplot as plt


def main():
    alpha1 = 0.2
    h = 2.1
    Y = (4.86, 2.17, 8.92, 4.91, 12.61, 12.80, 1.72, 0.95, 5.96, 10.50, 4.67, 11.39, 6.30, 6.02, 1.98, 8.35, 9.40,
         7.03, 0.01, 7.83, 0.57, 6.71, 9.17, 8.83, 13.45, 4.99, 8.29, 14.21, 7.43, 4.21, 2.21, 8.28, 4.16, 4.51,
         4.64, 4.60, 6.00, 7.38, 6.22, 5.77, 6.69, 8.77, 9.71, 12.76, 2.21, 5.36, 3.40, 7.92, 0.43, 10.82)
    X = (2, 2, 1, 3, 4, 2, 1, 3, 4, 3, 2, 1, 3, 1, 1, 2, 4, 1, 5, 3, 4, 1, 2, 2, 1, 3, 3, 4, 5, 1, 2, 3, 4, 3,
         3, 4, 2, 2, 2, 3, 5, 2, 3, 3, 3, 1, 4, 0, 4, 1)
    n = len(X)

    # Задание 1
    print("Задание 1")
    a1 = sum([x * x for x in X])
    b1 = sum(X)
    s1 = sum([X[i] * Y[i] for i in range(n)])
    a2 = b1
    b2 = n
    s2 = sum(Y)

    Delta = det2x2(a1, b1, a2, b2)
    beta1 = det2x2(s1, b1, s2, b2) / Delta
    beta0 = det2x2(a1, s1, a2, s2) / Delta

    print(f"f(x) = {beta1:.3f}x + {beta0:.3f}")
    x = linspace(min(X) - .5, max(X) + .5, 100)
    plt.plot(x, beta1 * x + beta0)
    plt.scatter(X, Y)
    plt.show()

    # Задание 2
    print("Задание 2")
    sigma = sum([(Y[i] - (beta1 * X[i] + beta0))**2 for i in range(n)]) / (n - 2)
    print(f"Несмещённая оценка дисперсии: {sigma}")

    E = list()
    for i in range(n):
        E.append(Y[i] - (beta1 * X[i] + beta0))
    maxE = max(E)
    minE = min(E)
    k = round((maxE - minE) / h)
    p = (h - (maxE - minE) / k) * k / 2
    plt.hist(E, bins=k, range=(minE - p, maxE + p))
    plt.show()

    x = linspace(minE - p, maxE + p, 100)
    norm_dis = norm(0, sqrt(sigma))
    plt.hist(E, density=True, bins=k, range=(minE - p, maxE + p))
    plt.plot(x, norm_dis.pdf(x), 'C1')
    plt.show()

    intervals = list({"center": minE - p + h * i + h / 2, "counter": 0} for i in range(k))
    for el in E:
        for i in range(k):
            if el < minE - p + h * (i + 1):
                intervals[i]["counter"] += 1
                break

    i_mean = sum([intervals[i]["center"] * intervals[i]["counter"] for i in range(k)]) / k
    i_deviation = sqrt(sum([intervals[i]["center"] * intervals[i]["center"] * intervals[i]["counter"]
                            for i in range(k)]) / k - i_mean * i_mean)
    i_deviation *= n / (n - 1)
    chi2observed = sum([(intervals[i]["counter"] - h * k / i_deviation
                         * norm.pdf((intervals[i]["center"] - i_mean) / i_deviation)) / intervals[i]["counter"]
                        for i in range(k)])
    chi2critical = chi2.ppf(1 - alpha1, df=k - 2)
    print(f"{chi2observed = :.2f} < {chi2critical = :.2f} => Нулевая гипотеза принимается")

    # Задание 5
    print("Задание 5")
    a1 = n
    b1 = sum(X)
    c1 = sum([x * x for x in X])
    s1 = sum(Y)
    a2 = b1
    b2 = c1
    c2 = sum([x * x * x for x in X])
    s2 = sum([X[i] * Y[i] for i in range(n)])
    a3 = b2
    b3 = c2
    c3 = sum([x * x * x * x for x in X])
    s3 = sum([X[i] * X[i] * Y[i] for i in range(n)])

    Delta = det3x3(a1, b1, c1, a2, b2, c2, a3, b3, c3)
    beta0 = det3x3(s1, b1, c1, s2, b2, c2, s3, b3, c3) / Delta
    beta1 = det3x3(a1, s1, c1, a2, s2, c2, a3, s3, c3) / Delta
    beta2 = det3x3(a1, b1, s1, a2, b2, s2, a3, b3, s3) / Delta

    print(f"f(x) = {beta2:.3f}x² + {beta1:.3f}x + {beta0:.3f}")
    x = linspace(min(X) - .5, max(X) + .5, 100)
    plt.plot(x, beta2 * x * x + beta1 * x + beta0)
    plt.scatter(X, Y)
    plt.show()

    # Задание 6
    print("Задание 6")
    sigma = sum([(Y[i] - (beta2 * X[i] * X[i] + beta1 * X[i] + beta0))**2 for i in range(n)]) / (n - 3)
    print(f"Несмещённая оценка дисперсии: {sigma}")

    E = list()
    for i in range(n):
        E.append(Y[i] - (beta2 * X[i] * X[i] + beta1 * X[i] + beta0))
    maxE = max(E)
    minE = min(E)
    k = round((maxE - minE) / h)
    p = (h - (maxE - minE) / k) * k / 2
    plt.hist(E, bins=k, range=(minE - p, maxE + p))
    plt.show()

    x = linspace(minE - p, maxE + p, 100)
    norm_dis = norm(0, sqrt(sigma))
    plt.hist(E, density=True, bins=k, range=(minE - p, maxE + p))
    plt.plot(x, norm_dis.pdf(x), 'C1')
    plt.show()

    intervals = list({"center": minE - p + h * i + h / 2, "counter": 0} for i in range(k))
    for el in E:
        for i in range(k):
            if el < minE - p + h * (i + 1):
                intervals[i]["counter"] += 1
                break
    i_mean = sum([intervals[i]["center"] * intervals[i]["counter"] for i in range(k)]) / k
    i_deviation = sqrt(sum([intervals[i]["center"] * intervals[i]["center"] * intervals[i]["counter"]
                            for i in range(k)]) / k - i_mean * i_mean)
    i_deviation *= n / (n - 1)
    chi2observed = sum([(intervals[i]["counter"] - h * k / i_deviation
                         * norm.pdf((intervals[i]["center"] - i_mean) / i_deviation)) / intervals[i]["counter"]
                        for i in range(k)])
    chi2critical = chi2.ppf(1 - alpha1, df=k - 2)
    print(f"{chi2observed = :.2f} < {chi2critical = :.2f} => Нулевая гипотеза принимается")


def det3x3(a1, b1, c1, a2, b2, c2, a3, b3, c3):
    return a1 * det2x2(b2, c2, b3, c3) - b1 * det2x2(a2, c2, a3, c3) + c1 * det2x2(a2, b2, a3, b3)


def det2x2(a1, b1, a2, b2):
    return a1 * b2 - a2 * b1


if __name__ == '__main__':
    main()