Python Hypothesis


  • 본문의 내용은 Data Science from scratch(밑바닥부터 시작하는 데이터 과학)을 읽고 작성했습니다.
  • 본문의 소스코드의 일부는 Joel Grus의 Github에서 Unlicensed 라이선스로 배포되고 있습니다.

가설과 추론

통계와 확률 이론을 통해 가설을 검정하고 통계적 추론을 할 수 있다.
가설(hypothesis)이란 ‘데이터 과학자는 MATLAB보다 R이나 Python을 선호한다.’ 등과 같은 주장이며, 데이터 통계에 대한 얘기가 될 수도 있다. 통계치들은 다양한 가정에서 특정 분포에 대한 확률변수의 관측치로 이해할 수도 있고, 가정이 얼마나 타당한지 알 수 있게 해주기도 한다.
고전적 가설검정에는 기본적인 가설 귀무가설(H0, null hypothesis)과 비교할 대립가설(H1, alternative hypothesis)로 구성된다. 통계를 통해 H0을 기각할지 말지를 정할 수 있다.

동전 던지기의 상황에 대해서 동전의 앞 뒷면이 나올 확률이 공평하다면 ‘p=0.5’가 귀무가설(H0)이 되며, p가 0.5가 아닌 경우가 대립가설이 된다. 동전을 n번 던져서 나온 횟수 X를 세는 것으로 검정을 진행할 수 있다. 베르누이 분포를 따르면 X가 이항분포를 따르는 확률변수라는 걸 의미한다. 또한 이항분포는 정규분포로 근사할 수 있다.

단측검정

가설을 ‘남자의 키가 클 것이다.’와 같이 세울 수 있다. 이처럼 가설을 세울 때 방향이 정해지는 것을 단측검정이라고 한다. 이미 가설의 방향을 알고 있거나 잠정적으로 결정한 상태에서 검증에 들어간다.

양측검정

반면에 같은 내용이라고 해도 ‘남녀 사이에 차이가 있을 것이다.’처럼 어느 쪽으로 방향을 정하지 않고 가설을 세울 때, 양측검정이라고 한다.

오류

1종 오류

귀무가설H0이 참일 때 H0을 기각하는 false positive(가양성) 오류를 뜻한다. 유의수준(significance) 는 제 1종 오류를 어느정도 허용할 것인가를 의미한다. 유의수준은 보통 5%나 1%로 설정하는 경우가 많다.

2종 오류

2종 오류를 범하지 않을 확률을 구하면 모델의 검정력(power)을 알 수 있다. 제 2종 오류란 H0가 거짓이지만 H0를 기각하지 않는 오류를 의미한다.

유의 수준에 맞게 가설을 기각하고 오류를 줄이면 검정력이 좋아졌다고 볼 수 있다.

p-value

가설을 바라보는 또 다른 관점은 p-value 이다. 이는 확률 값을 기준으로 구간을 선택하지 않고 H0가 참이라는 전제 하에 실제 관측된 값보다 더 극단적이 나올 확률을 구하는 것이다.

Wikipedia - p-value(유의 확률)은 귀무가설이 맞을 경우 대립가설 쪽의 값이 나올 확률을 나타내는 값. 확률 값이라고도 한다. 표본 평균이 귀무가설 값에서 멀수록 작아지게 된다.

p-value가 유의수준보다 작을수록 귀무가설을 지지하는 정도가 약하므로 귀무가설을 기각하게 되며(=대립가설 채택), p-value가 클수록 귀무가설을 지지하는 정도가 커짐으로 귀무가설을 채택하게 된다.

p-value의 모순?

“귀무가설이 ‘p=0.5’이고 유의수준 5%에서 극한값을 매겨서 p-value를 측정하면 p=0.5+-5%에서 멀어질 때마다 극한값이 카운트되니까 p-value 값이 높을 수록 귀무가설에서 멀어지는 게 아닐까?”

(유의수준 = H0을 기각하는 구간)
>> 95% 구간의 신뢰구간 바깥에 양 끝 꼬리 2.5%가 유의수준에 해당하는데, p-value가 5%보다 커야 P값이 5%의 유의수준에서 귀무가설을 기각하지 않을 때 신뢰구간에 들어갈 수 있다.

연속 수정(continuity correction)

In probability theory, a continuity correction is an adjustment that is made when a discrete distribution is approximated by a continuous distribution.
확률 이론에서 연속 수정은 이항분포를 연속분포로 근사했을 때 일어나는 수정을 뜻한다.
동전의 앞면이 530번 나올 확률을 확인하는 방법으로 529 - 531를 범위로 잡는 것보다 529.5 - 530.5를 범위로 잡는 것이 더 정확하다.

신뢰구간

사건에 대한 분포를 모를 때 관측된 값에 대한 신뢰구간을 사용해 가설을 검정할 수 있다. p의 정확한 값을 알고 있다면 중심극한정리를 이용해 베르누이 확률변수들의 평균은 대략 평균이 p이고 표준편차는 다음과 같은 정규분포로 추정할 수 있다.

math.sqrt(p * (1-q) / 1000)

정확한 p의 값을 모른다면 추정치를 사용할 수 있다.

import math

def normal_cdf(x, mu=0,sigma=1):
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2

def inverse_normal_cdf(p, mu=0, sigma=1, tolerance=0.00001):
    """find approximate inverse using binary search"""

    # if not standard, compute standard and rescale
    if mu != 0 or sigma != 1:
        return mu + sigma * inverse_normal_cdf(p, tolerance=tolerance)

    low_z, low_p = -10.0, 0            # normal_cdf(-10) is (very close to) 0
    hi_z,  hi_p  =  10.0, 1            # normal_cdf(10)  is (very close to) 1
    while hi_z - low_z > tolerance:
        mid_z = (low_z + hi_z) / 2     # consider the midpoint
        mid_p = normal_cdf(mid_z)      # and the cdf's value there
        if mid_p < p:
            # midpoint is still too low, search above it
            low_z, low_p = mid_z, mid_p
        elif mid_p > p:
            # midpoint is still too high, search below it
            hi_z, hi_p = mid_z, mid_p
        else:
            break


#####
#
# probabilities a normal lies in an interval
#
######

# the normal cdf _is_ the probability the variable is below a threshold
normal_probability_below = normal_cdf

# it's above the threshold if it's not below the threshold
def normal_probability_above(lo, mu=0, sigma=1):
    return 1 - normal_cdf(lo, mu, sigma)

# it's between if it's less than hi, but not less than lo
def normal_probability_between(lo, hi, mu=0, sigma=1):
    return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)

# it's outside if it's not between
def normal_probability_outside(lo, hi, mu=0, sigma=1):
    return 1 - normal_probability_between(lo, hi, mu, sigma)

######
#
#  normal bounds
#
######

def normal_upper_bound(probability, mu=0, sigma=1):
    """returns the z for which P(Z <= z) = probability"""
    return inverse_normal_cdf(probability, mu, sigma)

def normal_lower_bound(probability, mu=0, sigma=1):
    """returns the z for which P(Z >= z) = probability"""
    return inverse_normal_cdf(1 - probability, mu, sigma)

def normal_two_sided_bounds(probability, mu=0, sigma=1):
    """returns the symmetric (about the mean) bounds
    that contain the specified probability"""
    tail_probability = (1 - probability) / 2

    # upper bound should have tail_probability above it
    upper_bound = normal_lower_bound(tail_probability, mu, sigma)

    # lower bound should have tail_probability below it
    lower_bound = normal_upper_bound(tail_probability, mu, sigma)

    return lower_bound, upper_bound

def two_sided_p_value(x, mu=0, sigma=1):
    if x >= mu:
        # if x is greater than the mean, the tail is above x
        return 2 * normal_probability_above(x, mu, sigma)
    else:
        # if x is less than the mean, the tail is below x
        return 2 * normal_probability_below(x, mu, sigma)

def count_extreme_values():
    extreme_value_count = 0
    for _ in range(100000):
        num_heads = sum(1 if random.random() < 0.5 else 0    # count # of heads
                        for _ in range(1000))                # in 1000 flips
        if num_heads >= 530 or num_heads <= 470:             # and count how often
            extreme_value_count += 1                         # the # is 'extreme'

    return extreme_value_count / 100000

upper_p_value = normal_probability_above
lower_p_value = normal_probability_below


# p-value Hacking
#
# 동전 1000번 던지기 실험을 1000번 반복해 5%의 유의수준을 두고 가설검정을 한다.
import random

def run_experiment():
    return [random.random() < 0.5 for _ in range(1000)]

def reject_fairness(experiment):
    num_heads = len([flip for flip in experiment if flip])
    return num_heads < 469 or num_heads > 531

random.seed(0)
experiments = [run_experiment() for _ in range (1000)]
num_rejections = len([experiment for experiment in experiments if reject_fairness(experiment)])

print(num_rejections)
>>46

# 놀랍게도 p-value가 0.05 이하로 나오며, p-value의 관점에서 p=0.5의 가설은 기각된다. (하지만 실험을 반복했을 때 여전히 다른 결과가 나올 수 있다.)

A/B 테스트 해보기

Google Adsense를 관리한다고 가정해보자. 블로그의 포스트 끝 부분에 등록한 광고 A를 N_A명의 사용자가 보았고 n_A명의 사용자가 클릭했다. 또 블로그의 배너에 등록한 광고 B를 N_B명의 사용자가 보았고 n_B명의 사용자가 클릭했다.
각 사용자가 광고를 클릭할 확률을 pA라고 했을 때 n_A/N_A는 평균이 p_A이고 평균이 p_A이고 표준편차가

인 정규분포에 근접하다는 것을 알 수 있고, 확률을 p_B라고 했을 때 n_B/N_B는 평균이 p_B이고 평균이 p_B이고 표준편차가

인 정규분포에 근접하다는 것을 알 수 있다.

def estimated_parameters(N, n):
    p = n / N
    sigma = math.sqrt(p * (1 - p) / N)
    return p, sigma

def a_b_test_statistics(N_A, n_A, N_B, n_B):
    p_A, sigma_A = estimated_parameters(N_A, n_A)
    p_B, sigma_B = estimated_parameters(N_B, n_B)
    return (p_B - p_A) / math.sqrt(sigma_A ** 2 + sigma_B ** 2)

# 두 광고를 모두 1000명이 봤을 때 A는 200명이 클릭을, B는 180명이 클릭을 했다.
z = a_b_test_statistics(1000, 200, 1000, 180)
two_sided_p_value(z)
>>0.254141976542236

값이 크기 때문에 두 분포가 다르다고 결론지을 수 없다.
만약 150명이 B광고를 클릭했다면 통계는 다음처럼 나타난다.

z = a_b_test_statistics(1000, 200, 1000, 150)
two_sided_p_value(z)
>>0.003189699706216853
# 두 광고가 동일하게 효과적일 때 50만큼의 차이가 발생활 확률은 0.003에 불과하다.

베이지안 추론

여태까지의 방법은 통계적 검정을 확률적으로 설명하는 것에 (H_0이 사실이라면 이런 통계치가 발생할 확률은 3%이다.) 중점을 두었다. 또 다른 추론의 방법으로 알려지지 않은 파라미터를 변수로 보는 방법이 있다. 분석가에게 파라미터에 대한 사전분포가 주어지고, 관측된 데이터와 베이즈 정리를 사용해 사후분포를 갱신할 수 있다. 통계적 검정에 대한 확률적으로 결론을 내는 대신에 파라미터에 대해 확률적으로 결론을 낼 수 있다.
동전 던지기, 주사위 굴리기 처럼 알려지지 않은 파라미터가 ‘확률’이라면, 보통 모든 확률 값이 0에서 1 사이로 정의되는 베타분포를 사전분포로 사용한다.

이러한 과정을 통하면 가설에 대해 새로운 관점에서 확률적으로 말할 수 있게 된다. 예컨데 “관측된 데이터와 사전분포로 볼 때 동전의 앞면이 나올 확률이 49% ~ 51%일 경우는 5%밖에 되지 않는다.”가 그것이다. “동전이 공평하다면 이런 편향된 데이터를 관측할 가능성은 5%밖에 되지 않는다”와는 완전히 다른 관점이다.

# Beta 분포
def B(alpha, beta):
    """모든 확률 값의 합이 1이 되도록 해주는 정규화 값"""
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)

def beta_pdf(x, alpha, beta):
    if x < 0 or x > 1:
        return 0
    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)

대수의 법칙 놀이

from collections import Counter
big_number = []
for _ in range(10000):
    big_number.append(random.choice(["True","False"]))

Counter(big_number)
>>Counter({'False': 4989, 'True': 5011})