from __future__ import division
from collections import Counter
from functools import partial
from linear_algebra import dot, vector_add
from stats import median, standard_deviation, de_mean
from probability import normal_cdf
from gradient_descent import minimize_stochastic
#from simple_linear_regression import total_sum_of_squares
import math, random

def total_sum_of_squares(y):
    """the total squared variation of y_i's from their mean"""
    return sum(v ** 2 for v in de_mean(y))

"""
Explanation: 15. 다중 회귀 분석
End of explanation
"""


def predict(x_i, beta):
    """각 x_i의 첫 번째 항목은 1이라고 가정"""
    return dot(x_i, beta)

"""
Explanation: 14장의 내용을 추가 데이터를 사용해 모델의 성능을 높이기 위해,
더 많은 독립 변수를 사용하는 선형 모델을 시험
14장에서 다뤘던 모델
$y_i=\alpha+\beta x_i+\epsilon_i$
여기에 독립 변수를 추가하면
$\Rightarrow$시간(분)= $\alpha+\beta_1$(친구 수)$+\beta_2$(근무 시간)$+\beta_3$*(박사 학위 취득 여부)+$\epsilon$  
15.1 모델
각 입력값 $x_i$가 숫자 하나가 아니라 $k$개의 숫자인 $x_i1,\cdot\cdot\cdot,x_ik$라고 한다면,
다중 회귀 모델은 다음과 같은 형태를 띈다.
$y_i=\alpha+\beta_1 x_{i1}+\cdot\cdot\cdot+\beta_k x_{ik}+\epsilon_i$  
다중 회귀 분석에서는 보통 파라미터 벡터를 $\beta$라고 부름
여기에 상수항 $\alpha$까지 덧붙이려면, 각 데이터 x_i의 앞부분에 1을 덧붙이면 된다.
beta = [alpha, beta_1,$\cdot\cdot\cdot$, x_ik]  
그리고 각 데이터는 다음과 같이 된다.
x_i = [1, x_i1, $\cdot\cdot\cdot$, x_ik]
이렇게 하면 모델을 다음과 같이 나타낼 수 있다.
End of explanation
"""


def error(x_i, y_i, beta):
    return y_i - predict(x_i, beta)

"""
Explanation: 독립 변수 x는 다음과 같은 벡터들의 열로 표현할 수 있다.
[1,  # 상수항
 49, # 친구의 수
 4,  # 하루 근무 시간
 0]  # 박사 학위 취득 여부

15.2 최소자승법에 대한 몇 가지 추가 가정

$x$의 열은 서로 일차독립 해야 한다.
    이 가정이 성립하지 않는다면 $\beta$를 추정할 수 없다.
$x$의 모든 열은 오류 $\epsilon$과 상관관계가 없어야한다.
    이 가정이 위배되면 잘못된 $\beta$가 추정될 것이다.

독립 변수와 오류 사이에 상관관계가 존재한다면, 최소자승법으로 만들어지는 모델은 편향된 $\beta$를 추정해 준다.
15.3 모델 학습하기
오류 함수
End of explanation
"""


def squared_error(x_i, y_i, beta):
    return error(x_i, y_i, beta) ** 2

"""
Explanation: SGD를 사용하기 위한 오류 제곱 값
End of explanation
"""


def squared_error_gradient(x_i, y_i, beta):
    """i번째 오류 제곱 값의 beta에 대한 기울기"""
    return [-2 * x_ij * error(x_i, y_i, beta)
            for x_ij in x_i]

"""
Explanation: 만약 미적분을 알고 있다면, 오류를 직접 계산할 수도 있다.
End of explanation
"""


def estimate_beta(x, y):
    beta_initial = [random.random() for x_i in x[0]]
    return minimize_stochastic(squared_error, 
                               squared_error_gradient, 
                               x, y, 
                               beta_initial, 
                               0.001)            

"""
Explanation: SGD를 사용해서 최적의 베타를 계산
End of explanation
"""


x = [[1,49,4,0],[1,41,9,0],[1,40,8,0],[1,25,6,0],[1,21,1,0],[1,21,0,0],[1,19,3,0],[1,19,0,0],[1,18,9,0],[1,18,8,0],[1,16,4,0],[1,15,3,0],[1,15,0,0],[1,15,2,0],[1,15,7,0],[1,14,0,0],[1,14,1,0],[1,13,1,0],[1,13,7,0],[1,13,4,0],[1,13,2,0],[1,12,5,0],[1,12,0,0],[1,11,9,0],[1,10,9,0],[1,10,1,0],[1,10,1,0],[1,10,7,0],[1,10,9,0],[1,10,1,0],[1,10,6,0],[1,10,6,0],[1,10,8,0],[1,10,10,0],[1,10,6,0],[1,10,0,0],[1,10,5,0],[1,10,3,0],[1,10,4,0],[1,9,9,0],[1,9,9,0],[1,9,0,0],[1,9,0,0],[1,9,6,0],[1,9,10,0],[1,9,8,0],[1,9,5,0],[1,9,2,0],[1,9,9,0],[1,9,10,0],[1,9,7,0],[1,9,2,0],[1,9,0,0],[1,9,4,0],[1,9,6,0],[1,9,4,0],[1,9,7,0],[1,8,3,0],[1,8,2,0],[1,8,4,0],[1,8,9,0],[1,8,2,0],[1,8,3,0],[1,8,5,0],[1,8,8,0],[1,8,0,0],[1,8,9,0],[1,8,10,0],[1,8,5,0],[1,8,5,0],[1,7,5,0],[1,7,5,0],[1,7,0,0],[1,7,2,0],[1,7,8,0],[1,7,10,0],[1,7,5,0],[1,7,3,0],[1,7,3,0],[1,7,6,0],[1,7,7,0],[1,7,7,0],[1,7,9,0],[1,7,3,0],[1,7,8,0],[1,6,4,0],[1,6,6,0],[1,6,4,0],[1,6,9,0],[1,6,0,0],[1,6,1,0],[1,6,4,0],[1,6,1,0],[1,6,0,0],[1,6,7,0],[1,6,0,0],[1,6,8,0],[1,6,4,0],[1,6,2,1],[1,6,1,1],[1,6,3,1],[1,6,6,1],[1,6,4,1],[1,6,4,1],[1,6,1,1],[1,6,3,1],[1,6,4,1],[1,5,1,1],[1,5,9,1],[1,5,4,1],[1,5,6,1],[1,5,4,1],[1,5,4,1],[1,5,10,1],[1,5,5,1],[1,5,2,1],[1,5,4,1],[1,5,4,1],[1,5,9,1],[1,5,3,1],[1,5,10,1],[1,5,2,1],[1,5,2,1],[1,5,9,1],[1,4,8,1],[1,4,6,1],[1,4,0,1],[1,4,10,1],[1,4,5,1],[1,4,10,1],[1,4,9,1],[1,4,1,1],[1,4,4,1],[1,4,4,1],[1,4,0,1],[1,4,3,1],[1,4,1,1],[1,4,3,1],[1,4,2,1],[1,4,4,1],[1,4,4,1],[1,4,8,1],[1,4,2,1],[1,4,4,1],[1,3,2,1],[1,3,6,1],[1,3,4,1],[1,3,7,1],[1,3,4,1],[1,3,1,1],[1,3,10,1],[1,3,3,1],[1,3,4,1],[1,3,7,1],[1,3,5,1],[1,3,6,1],[1,3,1,1],[1,3,6,1],[1,3,10,1],[1,3,2,1],[1,3,4,1],[1,3,2,1],[1,3,1,1],[1,3,5,1],[1,2,4,1],[1,2,2,1],[1,2,8,1],[1,2,3,1],[1,2,1,1],[1,2,9,1],[1,2,10,1],[1,2,9,1],[1,2,4,1],[1,2,5,1],[1,2,0,1],[1,2,9,1],[1,2,9,1],[1,2,0,1],[1,2,1,1],[1,2,1,1],[1,2,4,1],[1,1,0,1],[1,1,2,1],[1,1,2,1],[1,1,5,1],[1,1,3,1],[1,1,10,1],[1,1,6,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],[1,1,4,1],[1,1,9,1],[1,1,9,1],[1,1,4,1],[1,1,2,1],[1,1,9,1],[1,1,0,1],[1,1,8,1],[1,1,6,1],[1,1,1,1],[1,1,1,1],[1,1,5,1]]
num_friends_good = [49,41,40,25,21,21,19,19,18,18,16,15,15,15,15,14,14,13,13,13,13,12,12,11,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,8,8,8,8,8,8,8,8,8,8,8,8,8,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
daily_minutes_good = [68.77,51.25,52.08,38.36,44.54,57.13,51.4,41.42,31.22,34.76,54.01,38.79,47.59,49.1,27.66,41.03,36.73,48.65,28.12,46.62,35.57,32.98,35,26.07,23.77,39.73,40.57,31.65,31.21,36.32,20.45,21.93,26.02,27.34,23.49,46.94,30.5,33.8,24.23,21.4,27.94,32.24,40.57,25.07,19.42,22.39,18.42,46.96,23.72,26.41,26.97,36.76,40.32,35.02,29.47,30.2,31,38.11,38.18,36.31,21.03,30.86,36.07,28.66,29.08,37.28,15.28,24.17,22.31,30.17,25.53,19.85,35.37,44.6,17.23,13.47,26.33,35.02,32.09,24.81,19.33,28.77,24.26,31.98,25.73,24.86,16.28,34.51,15.23,39.72,40.8,26.06,35.76,34.76,16.13,44.04,18.03,19.65,32.62,35.59,39.43,14.18,35.24,40.13,41.82,35.45,36.07,43.67,24.61,20.9,21.9,18.79,27.61,27.21,26.61,29.77,20.59,27.53,13.82,33.2,25,33.1,36.65,18.63,14.87,22.2,36.81,25.53,24.62,26.25,18.21,28.08,19.42,29.79,32.8,35.99,28.32,27.79,35.88,29.06,36.28,14.1,36.63,37.49,26.9,18.58,38.48,24.48,18.95,33.55,14.24,29.04,32.51,25.63,22.22,19,32.73,15.16,13.9,27.2,32.01,29.27,33,13.74,20.42,27.32,18.23,35.35,28.48,9.08,24.62,20.12,35.26,19.92,31.02,16.49,12.16,30.7,31.22,34.65,13.13,27.51,33.2,31.57,14.1,33.42,17.44,10.12,24.42,9.82,23.39,30.93,15.03,21.67,31.09,33.29,22.61,26.89,23.48,8.38,27.81,32.35,23.84]

random.seed(0)
beta = estimate_beta(x, daily_minutes_good)

beta # 분 = 30.63 + 0.972 친구 수 - 1.868 근무 시간 + 0.911 박사 학위 취득 여부

"""
Explanation: 실험 데이터 셋팅 (14장과 동일)
End of explanation
"""


def multiple_r_squared(x, y, beta):
    sum_of_squared_errors = sum(error(x_i, y_i, beta) ** 2
                                for x_i, y_i in zip(x, y))
    return 1.0 - sum_of_squared_errors / total_sum_of_squares(y)

multiple_r_squared(x, daily_minutes_good, beta)

"""
Explanation: 15.4 모델 해석하기
모델의 계수는 해당 항목의 영향력을 나타낸다.
친구 수 증가 =&gt; 대략 1분 증가  
근무 시간 증가 =&gt; 대략 2분 감소  
박사 학위 취득 =&gt; 대략 1분 증가

이러한 해석은 변수간의 관계를 직접적으로 설명해 주지 못한다.
예를 들어, 친구의 수가 다른 사용자들의 근무 시간은 서로 다를 수 있다.
이 모델은 이러한 관계를 잡아내지 못한다.
이러한 문제는 친구 수와 근무 시간을 곱한 새로운 변수로 해결할 수 있다.
변수가 점점 추가되기 시작하면, 각 계수가 유의미한지 살펴봐야 한다.
변수끼리 곱한 값, 변수의 log값, 변수의 제곱 값 등 수 많은 변수를 추가할 수 있기 때문이다.
15.5 적합성(Goodness of fit)
모델의 R 제곱 값을 다시 계산해 보자.
14장에서의 결과
alpha 22.94755241346903
beta 0.903865945605865
r-squared 0.3291078377836305
End of explanation
"""


def bootstrap_sample(data):
    """len(data)개의 항목을 중복을 허용한 무작위 추출"""
    return [random.choice(data) for _ in data]

def bootstrap_statistic(data, stats_fn, num_samples):
    """num_samples개의 bootstrap 샘플에 대해 stats_fn을 적용"""
    return [stats_fn(bootstrap_sample(data)) 
            for _ in range(num_samples)]

"""
Explanation: 회귀 분석 모델에 새로운 변수를 추가하면 R 제곱 값이 어쩔 수 없이 증가한다.
따라서 다중 회귀 분석 모델은 언제나 단순 회귀 분석 모델보다 작은 오류를 갖게 된다.
이러한 이유로 다중 회귀 분석 모델에서는 각 계수의 표준 오차를 살펴 봐야 한다.
계수의 표준 오차는 추정된 $\beta_1$의 계수가 얼마나 확실한지 알려준다.
표준오차 : 모집단 전체를 알 수 없는 상황에서, 모집단이 정상분포(정규분포)를 이루고 있다는 가정 하에, 여러번의 샘플링을 통해 각 표본 집단의 평균들로 이루어진 표준 평균 분포를 얻고, 이 표준 평균 분포의 표준 편차가 표준 오차가 된다.
표준오차는 모평균과 표본평균 사이의 오차를 알려주므로, 모집단의 표준편차가 클수록 표준오차도 커지고, 사례가 많을 수록 작아진다.

오차를 측정하기 위해서는 각 오류 $\epsilon_1$는 독립이며, 평균은 0이고 표준편차는 $\sigma$인 정규분포의 확률변수라는 가정이 필요하다.
표준오차가 클수록 해당 계수는 무의미해 진다.
15.6 여담: bootstrap
알 수 없는 분포에서 생성된 표본 데이터가 주어졌을 때,
이 표본 데이터의 중앙값을 찾으려면?
만약 표본 데이터가 모두 100 근처에 위치하고 있다면, 중앙값 또한 100 근처에 위치할 것이지만,
표본 데이터의 반은 0, 나머지 반은 200 근처에 위치하고 있다면,
추정된 중앙값을 신뢰하기 힘들다.
bootstrap은 중복이 허용된 재추출을 통해 새로운 데이터의 각 항목을 생성한다.
이를 통해 만들어진 데이터로 중앙값을 계산해 볼 수 있다.
End of explanation
"""


# 101개의 데이터가 모두 100에 인접
close_to_100 = [99.5 + random.random() for _ in range(101)]

# 101개의 데이터 중 50개는 0에 인접, 50개는 200에 인접
far_from_100 = ([99.5 + random.random()] +
               [random.random() for _ in range(50)] +
               [200 + random.random() for _ in range(50)])

print(close_to_100)
print(far_from_100)

"""
Explanation: 예를 들어, 다음과 같은 두 가지 데이터를 살펴보자.
End of explanation
"""


print(median(close_to_100))
print(median(far_from_100))

"""
Explanation: 만약 두 데이터의 중앙값을 계산해 보면 둘 다 대략 100에 가까운 것을 확인할 수 있다.
End of explanation
"""


print(bootstrap_statistic(close_to_100, median, 100))

print(bootstrap_statistic(far_from_100, median, 100))

"""
Explanation: 하지만 다음과 같이 bootstrap을 적용해 보면,
close_to_100은 100에 대부분 가깝지만,
far_from_100은 0 또는 200에 가까운 것을 확인할 수 있다.
End of explanation
"""


standard_deviation(bootstrap_statistic(close_to_100, median, 100))

standard_deviation(bootstrap_statistic(far_from_100, median, 100))

"""
Explanation: 첫번째는 표준편차가 0에 가깝지만, 두번째는 표준편차가 100에 가까운 것을 확인할 수 있다.
End of explanation
"""


def estimate_sample_beta(sample):
    x_sample, y_sample = zip(*sample) # magic unzipping trick
    return estimate_beta(x_sample, y_sample)

random.seed(0) # 예시와 동일한 결과를 얻기 위해 설정

bootstrap_betas = bootstrap_statistic(list(zip(x, daily_minutes_good)),
                                     estimate_sample_beta,
                                     100)

bootstrap_betas

"""
Explanation: 데이터가 이렇게 극단적인 경우에는 데이터를 직접 살펴보면 문제를 쉽게 파악할 수 있지만,
대부분의 경우 데이터만 살펴보는 것으로는 부족하다.
15.7 계수의 표준 오차
계수의 표준 오차를 추정할 때도 bootstrap을 적용할 수 있다.
bootstrap을 할 때에는, 하나의 데이터에 속하는 x와 y를 (x_i, y_i) 형태로 묶어줘야 하고,
반환된 데이터를 다시 x_sample, y_sample로 나눠줘야 한다.
End of explanation
"""


bootstrap_standard_errors = [
    standard_deviation([beta[i] for beta in bootstrap_betas])
    for i in range(4)
]

bootstrap_standard_errors

"""
Explanation: 그리고 각 계수의 표준 오차를 추정할 수 있다.
End of explanation
"""


def p_value(beta_hat_j, sigma_hat_j):
    if beta_hat_j > 0:
        return 2 * (1 - normal_cdf(beta_hat_j / sigma_hat_j))
    else:
        return 2 * normal_cdf(beta_hat_j / sigma_hat_j)

print(p_value(30.63, 1.174))
print(p_value(0.972, 0.079))
print(p_value(-1.868, 0.131))
print(p_value(0.911, 0.990))

"""
Explanation: 이제 '과연 $\beta_1$는 0일까?' 같은 가설을 검증해 볼 수 있다.
p-value(유의 확률)은 귀무가설이 맞을 경우 대립가설 쪽의 값이 나올 확률을 나타내는 값. 확률 값이라고도 한다. 표본 평균이 귀무가설 값에서 멀수록 작아지게 된다.
End of explanation
"""


# alpha는 패널티의 강도를 조절하는 하이퍼 파라미터
# 보통 "lamda"라고 표현하지만 파이썬에서는 이미 사용 중인 키워드이다.
def ridge_penalty(beta, alpha):
  return alpha * dot(beta[1:], beta[1:])

def squared_error_ridge(x_i, y_i, beta, alpha):
    """beta를 사용할 때 오류와 패널티의 합을 추정"""
    return error(x_i, y_i, beta) ** 2 + ridge_penalty(beta, alpha)

"""
Explanation: 대부분의 계수들은 0이 아닌 것으로 검증되었으나,
박사 학위 취득 여부에 대한 계수에 의미가 없을 수 있다는 것을 암시
15.8 Regularization

변수가 많아질수록 오버피팅
0이 아닌 계수가 많을수록 모델 해석이 어려움

Regularization은 beta가 커지면 커질수록 해당 모델에게 패널티를 주는 방법이다.
예를 들어, ridge regression의 경우, beta_i를 제곱한 값의 합에 비례하는 패널티를 추가한다.
하지만 상수에 대한 패널티는 주지 않는다.
End of explanation
"""


def ridge_penalty_gradient(beta, alpha):
    """패널티의 기울기"""
    return [0] + [2 * alpha * beta_j for beta_j in beta[1:]]

def squared_error_ridge_gradient(x_i, y_i, beta, alpha):
    """i번 오류 제곱 값과 패널티의 기울기"""
    return vector_add(squared_error_gradient(x_i, y_i, beta),
                      ridge_penalty_gradient(beta, alpha))

def estimate_beta_ridge(x, y, alpha):
    """패널티가 alpha인 리지 회귀를 경사 하강법으로 학습"""
    beta_initial = [random.random() for x_i in x[0]]
    return minimize_stochastic(partial(squared_error_ridge, alpha=alpha), 
                               partial(squared_error_ridge_gradient, 
                                       alpha=alpha), 
                               x, y, 
                               beta_initial, 
                               0.001)

"""
Explanation: 그리고 이전과 동일하게 경사 하강법을 적용할 수 있다.
End of explanation
"""


random.seed(0)
beta_0 = estimate_beta_ridge(x, daily_minutes_good, alpha=0.0)

beta_0

dot(beta_0[1:], beta_0[1:])

multiple_r_squared(x, daily_minutes_good, beta_0)

"""
Explanation: 만약 alpha가 0이라면 패널티는 전혀 없으며, 이전과 동일한 모델이 학습될 것이다.
End of explanation
"""


beta_0_01 = estimate_beta_ridge(x, daily_minutes_good, alpha=0.01)
print(beta_0_01)
print(dot(beta_0_01[1:], beta_0_01[1:]))
print(multiple_r_squared(x, daily_minutes_good, beta_0_01))

beta_0_1 = estimate_beta_ridge(x, daily_minutes_good, alpha=0.1)
print(beta_0_1)
print(dot(beta_0_1[1:], beta_0_1[1:]))
print(multiple_r_squared(x, daily_minutes_good, beta_0_1))

beta_1 = estimate_beta_ridge(x, daily_minutes_good, alpha=1)
print(beta_1)
print(dot(beta_1[1:], beta_1[1:]))
print(multiple_r_squared(x, daily_minutes_good, beta_1))

beta_10 = estimate_beta_ridge(x, daily_minutes_good, alpha=10)
print(beta_10)
print(dot(beta_10[1:], beta_10[1:]))
multiple_r_squared(x, daily_minutes_good, beta_10)

"""
Explanation: 그리고 alpha를 증가시킬수록 적합성은 감소하고, beta의 크기도 감소한다.
End of explanation
"""


def lasso_penalty(beta, alpha):
    return alpha * sum(abs(beta_i) for beta_i in beta[1:])    

"""
Explanation: 패널티가 증가하면 박사 학위 취득 여부 변수는 사라진다.
다른 형태의 패널티를 사용하는 lasso regression도 있다.
End of explanation
"""


