import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sp
## Розмір вибірки.
n = 50
## Значення предиктора.
XV = np.random.uniform(low=-4, high=4, size=n)
XV.sort()
## Матриця плану.
X = np.ones((n,2))
X[:,1] = XV
## Істинні коефіцієнти.
beta = np.array([0, 1.], dtype=np.float64)
## Істинні значення відгуку.
EY = np.dot(X, beta)
## Спостережувані значення відгуку.
Y = EY + np.random.normal(size=n)*np.sqrt(20)
## Отримати оцінки коефіцієнтів.
u,s,vt = np.linalg.svd(X,0)
v = np.transpose(vt)
bhat = np.dot(v, np.dot(np.transpose(u), Y)/s)
## Допасовані значення.
Yhat = np.dot(X, bhat)
## СКП та КСКП.
MSE = ((Y-EY)**2).sum()/(n-X.shape[1])
s = np.sqrt(MSE)
## Ці множники використовуються в побудові інтервалів.
XtX = np.dot(np.transpose(X), X)
V = [np.dot(X[i,:], np.linalg.solve(XtX, X[i,:])) for i in range(n)]
V = np.array(V)
## F-квантиль, що використовується в побудові інтервалу Шефе.
QF = sp.fdtri(X.shape[1], n-X.shape[1], 0.95)
## Нижня та верхня межі смуги Шефе.
D = s*np.sqrt(X.shape[1]*QF*V)
LB,UB = Yhat-D,Yhat+D
## Нижня та верхня межі поточкової смуги.
D = s*np.sqrt(2*V)
LBP,UBP = Yhat-D,Yhat+D
## Створити графік.
plt.clf()
plt.plot(XV, Y, 'o', ms=3, color='grey')
plt.plot(XV, EY, '-', color='black', label = "Істина")
plt.plot(XV, Yhat, '-', color='green', label = "Оцінка")
plt.plot(XV, LB, '-', color='red', label = "Одночасно 95 %-ва ДС")
plt.plot(XV, UB, '-', color='red')
plt.plot(XV, LBP, '-', color='blue', label = "Поточково 95 %-ва ДС")
plt.plot(XV, UBP, '-', color='blue')
plt.legend(frameon=False)
plt.ylim([-20,15])
plt.gca().set_yticks([-20,-10,0,10,20])
plt.xlim([-4,4])
plt.gca().set_xticks([-4,-2,0,2,4])
plt.xlabel("X")
plt.ylabel("Y")
plt.savefig("regression_confidence_band_uk.png")
plt.savefig("regression_confidence_band_uk.svg")