Добавил:
jetu
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз:
Предмет:
Файл:1. Исследование регрессионных моделей
.pyfrom 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()
Соседние файлы в предмете Статистический анализ