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

МИНОБРНАУКИ РОССИИ

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ

ЭЛЕКТРОТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

«ЛЭТИ» ИМ. В.И. УЛЬЯНОВА (ЛЕНИНА)

Кафедра алгоритмической математики

ОТЧЁТ

по лабораторной работе №1

по дисциплине «Статистический анализ»

Тема: «Исследование регрессионных моделей»

Студент гр. 93—	—
Преподаватель	Чирина А. В.

Санкт-Петербург

2021

Код программы.

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()

Результаты.

Задание 1

f(x) = -0.425x + 7.650

Задание 2

Несмещённая оценка дисперсии: 12.814562236853599

chi2observed = 6.27 < chi2critical = 7.29 => Нулевая гипотеза принимается

Задание 5

f(x) = -0.096x² + 0.082x + 7.131

Задание 6

Несмещённая оценка дисперсии: 13.06022926040737

chi2observed = 6.24 < chi2critical = 7.29 => Нулевая гипотеза принимается