from sklearn.datasets import make_regression

X0, y, coef = make_regression(n_samples=100, n_features=1, noise=20, coef=True, random_state=0)

dfX0 = pd.DataFrame(X0, columns=["X1"])
dfX = sm.add_constant(dfX0)
dfy = pd.DataFrame(y, columns=["y"])

model = sm.OLS(dfy, dfX)
result = model.fit()
print(result.params)

"""
Explanation: 확률론적 선형 회귀 모형
OLS(Ordinary Least Square) 방법을 사용하면 데이터에 대한 확률론적인 가정없이도 최적의 가중치를 계산할 수 있다. 그러나 이 경우에는 계산한 가중치가 어느 정도의 신뢰도 또는 안정성을 가지는지 확인할 수 있는 방법이 없다. 이를 확인하고자 하는 시도 중의 하나가 부트스트래핑(bootstrapping) 방법이다.
부트스트래핑
부트스트래핑(bootstrapping)은 회귀 분석에 사용한 데이터가 달라진다면 회귀 분석의 결과는 어느 정도 영향을 받는지를 알기 위한 방법이다. 
데이터가 확률 변수로부터 생성된 표본이거나 혹은 더 큰 모집단 중에서 선택한 표본이라고 가정한다면 회귀 분석의 결과는 분석에 사용한 표본에 의존적임을 알 수 있다. 만약 추가적인 다른 표본을 얻어서 다시 회귀 분석에 사용한다면 회귀 분석 결과 즉, 가중치 벡터의 값은 달라질 것이다.
그러나 현실적으로는 데이터를 추가적으로 얻기가 힘들기 때문에 부트스트래핑 방법에서는 기존의 데이터를 재표본화(re-sampling)하는 방법을 선택한다. 재표본화는 기존의 $D$개의 데이터에서 다시 $D$개의 데이터를 선택하되 중복 선택도 가능하게 한다. (resampling with replacement) 이 경우 이론적으로는 $2^D$개의 새로운 표본 집단을 얻을 수 있다.
직접 부트스트래핑을 실시해 보자. 우선 100개의 가상 데이터를 생성하여 이를 기반으로 회귀 분석을 실시한다.
End of explanation
"""


N = 1000
params_c = np.zeros(N)
params_x1 = np.zeros(N)
for i in range(N):
    idx = np.random.choice(len(dfy), len(dfy), replace=True)
    dfX2 = dfX.ix[idx, :]
    dfy2 = dfy.ix[idx]
    r = sm.OLS(dfy2, dfX2).fit()
    params_c[i] = r.params.const
    params_x1[i] = r.params.X1

"""
Explanation: 다음으로 이 데이터에서 중복을 허락하여 N개의 데이터를 선택한 후 다시 회귀 분석을 한다. 이론적으로 $2^{100}$개의 경우가 있지만 1,000번만 반복해 본다.
End of explanation
"""


ax1 = plt.subplot(121)
sns.distplot(params_c, ax=ax1)
plt.axvline(params_c.mean(), c='r')
plt.title("const parameter")
ax2 = plt.subplot(122)
sns.distplot(params_x1, ax=ax2)
plt.axvline(params_x1.mean(), c='r')
plt.title("x1 parameter")
plt.show()

"""
Explanation: 전체 가중치 집합을 히스토그램으로 나타내면 다음과 같다.
End of explanation
"""


sp.stats.describe(params_c)

sp.stats.describe(params_x1)

"""
Explanation: 평균과 분산은 다음과 같다.
End of explanation
"""


print(result.summary())

"""
Explanation: 가중치 중 상수항의 경우 평균은 -1.6이지만 표분 편차가 $\sqrt{4.81}=2.19$이므로 0일 가능성을 배제할 수 없다.
이 결과를 statsmodels 의 회귀 분석 보고서와 비교하자.
End of explanation
"""


sp.stats.probplot(result.resid, plot=plt)
plt.show()

plt.plot(X0, result.resid, 'o')
plt.axhline(y=0, c='k')
plt.xlabel("X1")
plt.ylabel("Residual")
plt.show()

import statsmodels.stats.api as sms
test = sms.omni_normtest(result.resid)
for x in zip(['Chi^2', 'P-value'], test):
    print("%-12s: %6.3f" % x)

import statsmodels.stats.api as sms
test = sms.jarque_bera(result.resid)
for x in zip(['Jarque-Bera', 'P-value', 'Skew', 'Kurtosis'], test):
    print("%-12s: %6.3f" % x)

"""
Explanation: 보고서의 std err 항목을 보면 표준 편차의 경우 2.163 이고 마지막의 신뢰 구간(confidence interval)이 -5.920 ~ 2.663 임을 보이고 있다. 부트스트래핑으로 얻은 결과와 유사하다. 이 결과는 다음에 설명할 확률론적 가정에 의해 계산된 값이다.
확률론적 선형 회귀 모형
확률론적 선형 회귀 모형에서는 $y$가 확률 변수로부터 생성된 표본이라고 가정하며 다음과 같은 조건을 만족한다.

선형 정규 분포 가정
종속 변수 $y$는 기댓값 $w^Tx$, 분산 $\sigma^2$ 를 가지는 정규 분포 확률 변수이다.
   $$ p(y \mid x, \theta) = \mathcal{N}(y \mid w^Tx, \sigma^2 ) $$ 

따라서 오차(disturbance)  $ \epsilon = y-w^Tx $ 도 정규 분포 확률 변수이다.
   $$ p(\epsilon \mid \theta) = \mathcal{N}(0, \sigma^2 ) $$ 


외생성(Exogeneity) 가정


오차 $\epsilon$와 독립 변수 $x$는 서로 독립이다. 즉, $ \text{E}[\epsilon] = 0$ 은 $x$의 값에 상관없이 유지된다.
   $$ \text{E}[\epsilon \mid x] = 0$$ 


조건부 독립 가정

$i$번째 샘플의 오차 $\epsilon_i$와 $j$번째 샘플의 오차 $\epsilon_j$는 $x$에 대해 서로 조건부 독립이다.
   $$ \text{Cov}[\epsilon_i, \epsilon_j \mid x] =  0$$ 

샘플의 수가 복수이면 샘플 오차 벡터 $\epsilon$에 대해 다음과 같이 쓸 수 있다. 아래 식에서의 $\epsilon$은 샘플 수와 같은 길이의 벡터이다.
   $$  \text{Cov}[ \epsilon \epsilon^T \mid x ] = \text{Cov}[ \epsilon \epsilon^T ] = \sigma^2 I $$
MLE를 사용한 선형 회귀 분석
앞의 확률론적 선형 회귀 모형과 MLE(Maximum Likelihood Estimation)을 사용하여 가중치 벡터 $w$의 값을 구해보자.
Likelihood는 다음과 같다.
$$
\begin{eqnarray}
p(y_{1:N} \,\big|\, x_{1:N}, \theta)
&=& \prod_{i=1}^N N(y_i \,\big|\, w^T x_i , \sigma^2) \
&=& \prod_{i=1}^N \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left{-\frac{(y_i-w^T x_i)^2}{2\sigma^2} \right}  \
\end{eqnarray}
$$
계산을 용이하기 위해 Log를 취하면 다음과 같다.
$$
\begin{eqnarray}
\text{LL}
&=& \log p(y_{1:N} \,\big|\, x_{1:N}, \theta) \
&=& \log \prod_{i=1}^N \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left{-\frac{(y_i-w^T x_i)^2}{2\sigma^2} \right}  \
&=& -\dfrac{1}{2\sigma^2} \sum_{i=1}^N (y_i-w^T x_i)^2 + \dfrac{1}{2} \sum_{i=1}^N  \log{2\pi}{\sigma^2}  \
\end{eqnarray}
$$
이를 행렬로 표시하면 다음과 같다.
$$
\text{LL}  =  -C_1 (y - Xw)^T(y-Xw) + C_0 = -C_1(w^TX^TXw -2 y^TXw + y^Ty) + C_0
$$
$$
C_1 =  -\dfrac{1}{2\sigma^2}
$$
$$
C_0 =  \dfrac{1}{2} \sum_{i=1}^N  \log{2\pi}{\sigma^2} 
$$
이를 최적화하면 OLS와 동일한 결과를 얻을 수 있다.
$$
\dfrac{\partial}{\partial w} \text{LL} \propto - 2X^TX \hat{w} + 2X^Ty = 0
$$
$$
\hat{w} = (X^TX)^{-1}X^T y
$$
잔차의 분포
위의 확률론적 선형 회귀 모형에 따르면 잔차 $e = y - \hat{w}^Tx$ 도 정규 분포를 따른다. 이는 다음과 같이 증명할 수 있다.
확률론적 선형 회귀 모형의 오차 $\epsilon$와 잔차 $e$는 다음과 같은 관계를 가진다.
$$  \hat{y} = X\hat{w} = X (X^TX)^{-1}X^T y = Hy $$
$$ e = y - \hat{y}= y - Hy = (I - H) y$$
$M = I - H$이라고 정의하면
$$ e = My = M (Xw + \epsilon) $$
최적화 조건에서
$$
X^TX \hat{w} - X^Ty = 0
$$
$$
X^T(X\hat{w} - y) = -X^Te = 0
$$
$$
X^TMy = 0 
$$
이 식은 모든 $y$에 대해 성립하므로
$$
X^TM = 0
$$
$H$가 대칭 행렬이므로 $M = I -H$도 대칭 행렬
$$
MX = 0
$$
$$ e = MXw + M\epsilon = M\epsilon $$
즉, 잔차 $e$는 오차 $\epsilon$의 선형 변환(linear transform)이다. 정규 분포의 선형 변환은 마찬가지로 정규 분포이므로 잔차도 정규 분포를 다른다.
End of explanation
"""


