파이썬으로 배우는 통계학
0. 통계적 사고¶
아래 내용은 통계적 사고: 파이썬을 이용한 탐색적 자료 분석 을 읽고 내용을 정리한 것이다. 이 책의 원제는 Thinkstat2(저자 Allen B. Downey) 로 저자는 올린공대 컴퓨터과학과의 교수이다.
1. 탐색적 자료 분석¶
일화적 증거(anecdotal evidence) 는 개인적인 경험을 직접, 혹은 한 다리 건너 듣고 이야기하는 것이다. 일화적 증거의 예시로 “최근에 첫째 아이를 출산한 내 친구 두명은 모두 예정일을 2주 지나서 낳았다. 그러니 내 첫째 아이도 2주가 지나서 나올 것이다."
이런 일화적 증거의 신빙성이 부족한 이유는 다음과 같다:
- 적은 관측치(Small number of observations)
- 선택편의(Selection bias)
- 확증편의(Confirmation bias)
- 부정확함(Inaccuracy)
1.1. 통계적 접근법¶
통계적인 접근법은 다음을 포함한다.
- 자료수집(Data collection): 통상적으로 대규모의 국가적인 조사에서 나온 자료를 사용한다.
- 기술통계(Descriptive statistics): 데이터를 간결하게 요약한 통계량(표본의 몇몇 특징을 수치화한 값이다)을 생성하고 데이터를 시각화 한다.
- 탐색적 데이터분석(Exploratory data analysis): 데이터의 패턴, 차이점, 특이값을 찾는다.
- 추정(Estimation): 사용한 데이터를 토대로 모집단의 특징을 추정한다.
- 가설검증(Hypothesis testing): 데이터셋에서 두 그룹간의 차이가 명백한 것인지 혹은 우연한 사건인지 평가한다.
1.2. 가족성장 국가조사(National Survey of Family Growth)¶
미국질병통제예방센터(Disease Control and Prevention, CDC)는 1973년부터 가족성장 국가조사(National Survey of Family Growth, 이하 NSFG)를 수행하고 있다.
NSFG 조사 목적은“가족생활, 결혼및이혼, 임신, 출산, 피임, 그리고 남녀건강에 대한정보를 수집하는데 있다. 자세한 사항은 공식 웹사이트를 참고하라.
NSFG는 횡단적연구(cross-sectional study)로 특정시점에 집단에 대한 정보를 수집한다. 지금까지 NSFG 조사는 총 7번 수행되었으며 여기에서는 2002년1월에서 2003년 3월까지 수행된 6번째 데이터셋(2002FemPreg.dat)를 사용한다.
1.3. 데이터 가져오기¶
NSFG에서 제공하는 데이터(2002FemPreg.dat)는 Stata 에서 사용되는 파일 형식이다. Stata 파일 형식은 일종의 딕셔너리 형식으로 각 행마다 각 변수의 위치를 식별하는데 사용되는 인덱스, 형식, 변수 정보를 담고 있다.
2002FemPreg.dat파일은 NSFG 홈페이지에서 다운로드 받을 수 있다.
1.4. 데이터프레임¶
원서에서는 2002FemPreg.dat파일을 판다스(pandas) 데이터프레임으로 변환하는 코드를 제공하고 있다. 그러나 여기에서는 편의를 위해 이미 변환된 2002FemPreg.csv파일을 사용한다.
1.5. 변수(Variables)¶
NSFG 데이터셋에는 총 244개의 변수가 있다. 그 중에 중요한 몇가기 변수를 아래에 설명했다.
-
caseid는 응답자 ID로 정수형 숫자이다. -
prglngth는 임신기간을 나타낸다. -
outcome은 출산결과에 대한 정수형 코드값이다. 1은 정상출산을 나타낸다. -
pregordr은 임신에 대한 일련번호다. 음답자의 첫번째 임신은 1이고 두번째 임신은 2이다. -
irthord는 정상적인 출산에 대한 일련번호다. 응답자의 첫번째 아이는 1이며 만약 정상적인 출산이 아닌 경우에는 공백으로 처리한다. -
birthwgt_lb,birthwgt_oz은 출산한 아기의 체중(파운드와 온스) 정보이다. -
agepreg는 임신한 산모의 나이이다. -
finalwgt는 응답자와 연관된 통계적 가중치로 미국인구중 응답자가 대표하는 비중을 나타낸다.
1.6. 변환(Transformation)¶
데이터 정제(data cleaning)은 데이터를가져올때, 오류를점검하고, 특수값을처리하고, 데이터를 다른형식으로변환하고, 계산을 수행하는 것을 말한다.
1.7. 타당성 검사(Validation)¶
데이터를 다른 소프트웨어환경으로 내보내고 또 다른 소프트웨어환경에서 가져올때 오류가 발생하곤 한다. 따라서 데이터의 타당성을 확인하는 단계를 거친다면 나중에 업무시간을 절약하고 오류를 피하는데 도움이 된다.
데이터 타당성을 확보하는 가장 간단한 방법은 기초통계량을 계산해 결과값을 비교하는 것이다. NSFG 데이터의 경우 다음 표와 같은 정보가 이미 알려져 있다.
| value | label | Total |
|---|---|---|
| 1 | LIVE BIRTH | 9148 |
| 2 | INDUCED ABORTION | 1862 |
| 3 | STILLBIRTH | 120 |
| 4 | MISCARRIAGE | 1921 |
| 5 | ECTOPIC PREGNANCY | 190 |
| 6 | CURRENT PREGNANCY | 352 |
우리가 사용할 데이터셋에서도 확인해 보자.
import warnings
warnings.filterwarnings("ignore") # 경고문 끄기
import pandas as pd
df = pd.read_csv("../data/2002FemPreg.csv")
df["outcome"].value_counts().sort_index()
우리가 알고 있던 값과 동일한 것을 확인 할 수 있다.
1.8. 해석(Interpretation)¶
효과적으로 데이터를 다루기 위해 통계적 관점과 문맥적 관점을 동시에 생각해야 한다. 예를 들어 caseid가 10229 인 여성의 경우를 살펴본다.
df[df["caseid"] == 10229].outcome.values
위의 결과 값에서 숫자 4는 유산을 뜻하고 1은 정상적인 출산을 나타낸다. 즉, 위 여성은 총 6번 유산했고 마지막 7번째에 정상적으로 출산을 했다는 것을 알 수 있다.
1.9. 연습문제¶
- NSFG 데이터셋의 신생아 몸무게의 분포를 나타내라.
df.birthwgt_lb.value_counts().sort_index()
위 결과를 통해 대부분의 신생아는 2.75 - 3.17 사이에 모여 있다는 것을 알 수 있다.
- 신생아 몸무게의 평균은 얼마인가?
df.birthwgt_lb.mean()
1.10. 용어설명¶
- 일화적 증거(anecdotal evidence): 과학적인 조사에 의한 것이 아닌 우연하게 수집된 증거.
- 모집단(population): 통계 조사에서 관심을 갖는 집단.
- 종단적연구(cross-sectional study): 특정시점에 모집단에 대한 자료를 수집하는 연구.
- 횡단적연구(longitudinal study): 시간을 두고 모집단을 추적하는연구, 동일한 그룹에서 반복적으로 데이터를 수집한다.
- 레코드(record): 데이터셋에서 하나에 대한 정보
- 표본(sample): 자료 수집하는데 사용된 모집단의 부분집합.
- 대표성(representative): 만약 모집단의 모든 멤버가 표본에 뽑힐 가능성이 동일한다면 대표성이 있다고 한다.
- 오버샘플링(oversampling): 적은 표본 크기로 생기는 오류를 피하기 위해 사용되는 방법
- 원시자료(raw data): 가장 처음에 수집되어 계산및 해석이 전혀 없는 상태의 데이터.
- 재코드(recode): 원시자료에 특정 계산 혹은 다른방법을 통해 수정된 데이터.
- 자료정제(data cleaning): 데이터 타당성 확보, 오류식별, 자료형간의 변환등을 포함하는 과정.
2. 분포(Distribution)¶
2.1. 히스토그램¶
데이터의 분포를 표현할때 가장 일반적으로 사용하는 것은 히스토그램(histogram)이다. 히스토그램은 각 값의 빈도(frequency)를 보여주는 그래프다.
2.3. 간단한 히스토그램 그리기¶
1장에서 살펴본 신생아들의 몸무게에 대한 히스토그램을 그려본다. 판다스의 plot.hist()함수를 사용한다.
df.birthwgt_lb.plot.hist(bins=14)
df.agepreg.plot.hist(bins=100)
산모 나이에 대한 히스토그램의 분포는 꼬리가 좀 더 오른쪽으로 뻗어나갔다. 대부분 의 산모는 20대이지만 30대이후에도 적지않게 존재한다는 의미이다.
2.4.2. 임신기간 히스토그램¶
df[df["outcome"] == 1].prglngth.plot.hist(bins=50)
위의 코드에서 df['outcome']==1는 정상적으로 출산한 경우를 선택하기 위해 사용했다.
임신기간 히스토그램의 단위는 주간(week)으로 가장 흔한값은 39주이다. 그리고 왼쪽 꼬리가 더 길다. 이것은 조산에 의한 것이이다. 그리고 임신기간이 43주가 넘어가면 유도분만을 하기 때문에 오른쪽 꼬리는 짧다.
2.5. 특이값(Outliers)¶
특이값은 일반적인 데이터와 다르게 극단에 존재하는 데이터이다. 특이값은 측정의 오류에 의한 것일 수도 있고 드물게 실제 데이터일수도 있다. 앞에서 배운 히스토그램은 최빈값과 데이터 분포의 모양을 식별하기에는 적합하지만 특이값을 확인하기 어렵다.
특이값을 찾는 가장 좋은 방법은 해당 분야의 전문지식(domain knowledge)이다. 다시말해 데이터 출처와 데이터의 의미를 파악할 수 있어야 특이값을 찾는데 어떤 분석법을 사용할지 알 수 있다.
2.6. 첫째 아이¶
NSFG 데이터셋으로 첫째 아이가 일찍(혹은 늦게) 태어나는 경향이 있는가?를 살펴본다. 먼저 첫째 아이와 나머지 아이들에 대한 임신 기간의 분포를 비교해보자.
import matplotlib.pyplot as plt
%matplotlib inline
live_baby = df[df["outcome"] == 1] # 출산에 성공한 아기
first_baby = live_baby[live_baby["birthord"] == 1] # 첫째 아이
others = live_baby[live_baby["birthord"] != 1] # 나머지 아이
first_baby.prglngth.plot.hist(bins=50, alpha=0.5)
others.prglngth.plot.hist(bins=50, alpha=0.5)
plt.legend(["first_baby", "others"])
히스토그램은 최빈값을 명확하게 보여준다는 점에서 유용하지만 두 그룹간의 분포를 비교하는데는 좋은 선택이 되지 못한다는 것을 알 수 있다.
또한 두 그룹의 높이 차이는 사실 표본크기에 의한 것이다. 위 예제에서“첫째가 아닌아이”보다“첫째아이”숫자가 더 적기 때문이다. 그래서 다음 3장에서 확률질량함수(probability mass functions)를 사용해 이 문제에 대해 다시 살펴본다.
2.7. 평균¶
가장 흔한 요약통계(summary statistics) 는 평균(mean) 으로 분포의 중심경향성을 기술할기 위해 사용된다.
$$\bar x = \frac{1}{n} \sum_{i} x_i$$
“평균(mean)”과“평균(average)”은 때때로 상호 호환적으로 사용된다. 하지만 이 책에서는 다음과 같이 구별한다.
- 표본평균(mean)은 위의 공식으로 계산되는 요약 통계이다.
- 평균(average)은 중심 경향성을 기술하기 위해 사용하는 요약 통계량중 하나다.
2.8. 분산¶
분산은 분포의 변동(variability)과 퍼짐(spread)을 기술하는 요약 통계다. 공식으로 나타낸 분산은 다음과 같다.
$$S^2 = \frac{1}{n} \sum_{i} (x_i - \bar x)^2$$
$(x_i - \bar x)^2$항은“평균으로부터 편차(deviation from the mean)”라고 하며 분산은 평균 편차의 제곱이다. 또한 분산의 제곱근(S)을 표준편차(standard deviation) 이라 부른다.
2.9. 효과크기(Effect size)¶
효과크기(effect size) 는 두 집단간의 차이를 기술하는 요약 통계로 현상이 실제로 모집단에 존재하는 정도을 말한다.
효과크기를 계산하는 방법으로 데이터셋의 집단 간(between groups)의 차이를 집단 내(within groups) 변동성과 비교 하는 코헨(Cohen) d 값을 주로 사용한다. 코헨(Cohen) d 는 다음과 같이 정의된다.
$$d = \frac{\bar x_1 - \bar x_2}{s}$$
$\bar x_1 - \bar x_2$은 집단 평균값이고, $s$는 합동 표준편차(pooled standard deviation)다.
효과크기가 0이라는 것은 비교 집단들 사이의 차이가 없다는 것을 의미하며 귀무 가설이 성립하게 된다. 평균치 비교의 경우 비교 하려는 집단 사이에 평균 차이가 클수록 효과크기는 커지게 된다. 효과크기의 쉬운 예를 들어보면 남녀 성별 비율이 50:50이라는 귀무가설하에서 실제 비율은 53:47이라고 하면 효과크기는 3%가 된다.
즉, 효과크기가 0이라는 의미는 비교하려는 집단사이의 평균이 동일하다는뜻이고, 효과크기의 값이 양수를 갖게 되면 비교집단이 대조집단에 비해 평균치가 크다는 의미이다. 그리고 음수 값을 갖게되면 비교집단의 평균이 대조집단에 비해 작다는것을 의미한다.
2.11. 연습문제¶
- 변수
totalwgt_lb을사용해 첫째 아이가 다른 아이들 보다 더가벼운지더무거운지조사하시오. 집단사이 차이를 정량화하는데 코헨d 값을 계산한다.
import math
def CohenEffectSize(group1, group2):
diff = group1.mean() - group2.mean()
var1 = group1.var()
var2 = group2.var()
n1, n2 = len(group1), len(group2)
pooled_var = (n1 * var1 + n2 * var2) / (n1 + n2)
d = diff / math.sqrt(pooled_var)
return d
cohen_d = CohenEffectSize(first_baby.totalwgt_lb, others.totalwgt_lb)
print(f"코헨d 값은 {cohen_d} 이다.")
위의 결과가 음수이기 때문에 첫째 아이의 몸무게 평균은 나머지 아이들의 평균보다 작다는 것을 알 수 있고, 값이 0에 가깝기 때문에 효과 크기가 아주 작다는 것을 알 수 있다.
코헨은 효과 크기를 임의로 small(0.5미만), medium(0.5이상), large(1이상)로 구분해 직관적으로 이해하기 쉽게 소개하였다.
2.12. 용어사전¶
- 분포(distribution): 표본에 나타나는 값과 개별값의 빈도
- 히스토그램(histogram): 데이터셋의 값들을 빈도로 나타내는 그래프
- 빈도(frequency): 데이터셋에서 값이 나타나는 횟수.
- 최빈값(mode): 데이터셋에서 가장 빈도가 높은값
- 정규분포(normal distribution): 종모양의 분포, 가우스분포로도 알려져있다.
- 균등분포(uniform distribution): 모든값이 동일한 빈도로 분포.
- 꼬리(tail): 분포의 양극단에 존재하는 부분
- 특이값(outlier): 중심에서 많이 떨어진 값.
- 퍼짐(spread): 데이터들이 얼마나 퍼져있는지에 대한 측정값.
- 요약통계(summary statistic): 중심 경향성 혹은 퍼짐 같이 분포의 측면을 정량화하는 통계치.
- 분산(variance): 퍼짐을 정량화하는데 사용되는 요약통계.
- 표준편차(standard deviation): 분산의 제곱근으로 퍼짐을 측정하는데 사용된다.
- 효과크기(effect size): 집단간의 차이처럼 효과크기를 정량화하는 요약통계.
3. 확률 질량 함수¶
3.1. PMF 소개¶
확률 질량 함수(Probability Mass Function; 이하 PMF)은 분포를 표현하는 또 다른 방식이다. PMF는 각각의 값을 확률로 치환하는 방법으로 각각의 확률을 얻기 위해 표본의 크기 n으로 값을 나누는데 이것을 정규화 라고 한다.
여기서는 설명을 위해 저자가 작성한 thinkstats2 모듈을 사용해 PMF를 계산한다. 실무에서는 scipy를 사용하는편이 성능면에서 유리할 것이다.
import thinkstats2
pmf = thinkstats2.Pmf([1, 2, 2, 3, 5])
print(pmf)
import thinkplot
first_pmf = thinkstats2.Pmf(first_baby.prglngth)
other_pmf = thinkstats2.Pmf(others.prglngth)
width = 0.45
thinkplot.PrePlot(2, cols=2)
thinkplot.Hist(first_pmf, align="right", width=width, label="first")
thinkplot.Hist(other_pmf, align="left", width=width, label="others")
thinkplot.Config(xlabel="weeks", ylabel="probability", axis=[27, 46, 0, 0.6])
thinkplot.PrePlot(2)
thinkplot.SubPlot(2)
thinkplot.Pmfs([first_pmf, other_pmf])
thinkplot.Show(xlabel="weeks", axis=[27, 46, 0, 0.6])
3.3. 다른 시각화 방법¶
히스토그램과 PMF은 데이터를 탐색하고 패턴을 식별하는데 유용하다. 위에서 살펴본 첫째아이와 나머지 아이 그룹의 가장 큰차이는 최빈값에 있다. 최빈값 차이를 강조하는 시각화를 해보자.
weeks = range(35, 46)
diffs = []
for week in weeks:
p1 = first_pmf.Prob(week)
p2 = other_pmf.Prob(week)
diff = 100 * (p1 - p2)
diffs.append(diff)
thinkplot.Bar(weeks, diffs)
위 막대 그래프의 y축은 확률이다. 따라서 첫째 아이는 임신 39주차에 덜 태어났고 41-42주차에 좀 더 태어났다는 것을 알 수 있다.
3.4. 학급 크기 역설(Class size paradox)¶
학교에서 학급(Class)를 개설했는데, 학교 측에서 계산한 학급 별 평균 인원보다 학생 측에서 느끼는 학급 별 평균 인원이 더 많은 현상을 말한다.
3.5. 데이터프레임 인덱싱(indexing)¶
자세한 것은 판다스(pandas) 공식 문서를 참고하라.
3.6. 연습문제¶
- 달리기 경주 결과(
Apr25_27thAn_set1.shtml)를 사용해 위에서 배운 학급 크기 역설을 적용한다. 주최자 입장에서의 실제 참가자 속도 와 참가자들의 입장에서 주자들의 속도 를 PMF로 시각화하라.
import numpy as np
def CleanLine(line):
t = line.split()
if len(t) < 6:
return None
place, divtot, div, gun, net, pace = t[0:6]
if "/" not in divtot:
return None
for time in [gun, net, pace]:
if ":" not in time:
return None
return place, divtot, div, gun, net, pace
def ConvertPaceToSpeed(pace):
m, s = [int(x) for x in pace.split(":")]
secs = m * 60 + s
mph = 1 / secs * 60 * 60
return mph
def ReadResults(filename="../data/Apr25_27thAn_set1.shtml"):
results = []
for line in open(filename):
t = CleanLine(line)
if t:
results.append(t)
return results
def GetSpeeds(results, column=5):
speeds = []
for t in results:
pace = t[column]
speed = ConvertPaceToSpeed(pace)
speeds.append(speed)
return speeds
def BinData(data, low, high, n):
data = (np.array(data) - low) / (high - low) * n
data = np.round(data) * (high - low) / n + low
return data
results = ReadResults()
speeds = GetSpeeds(results)
speeds = BinData(speeds, 3, 12, 100)
실제 참가자들의 속도
pmf = thinkstats2.Pmf(speeds, "actual speeds")
thinkplot.Pmf(pmf)
thinkplot.Config(xlabel="Speed (mph)", ylabel="PMF")
경주 참가자들이 느끼는 다른 참가자의 상대 속도
def ObservedPmf(pmf, speed, label=None):
new = pmf.Copy(label=label)
for val in new.Values():
diff = abs(val - speed)
new[val] *= diff
new.Normalize()
return new
biased = ObservedPmf(pmf, 7, label="observed speeds")
thinkplot.Pmf(biased)
thinkplot.Config(xlabel="Speed (mph)", ylabel="PMF")
3.7. 용어사전¶
- 확률질량함수(Probability mass function, PMF): 값을 확률로 치환해 분포로 표현.
- 확률(probability): 표본 크기의 일부로 표현되는 빈도수.
- 정규화(normalization): 확률값을 얻기위해서 표본크기로 빈도수를 나누는과정.
4. 누적 분포 함수¶
4.1. PMF의 한계¶
PMF는 데이터의 수가 적을때는 잘 작동하지만 데이터의 개수가 많이 증가되면 각각의 확률값이 작아져 확률잡음의 효과가 증가하는 한계가 있다. 따라서 이 문제를 해결하기 위해 누적 분포 함수(cumulative distribution function, 이하 CDF)의 개념이 나오게 되었다. CDF를 설명하기 앞서 백분위수(percentile)을 살펴보자.
4.2. 백분위수(Percentiles)¶
백분위수는 데이터를 순서대로 나열했을 때 백분율로 특정 위치의 값을 나타내는 것이다. 일반적으로 가장 작은 것을 0, 가장 큰 것을 100으로 한다.
4.3. CDF¶
누적 분포 함수(CDF)는 데이터를 백분위 순위로 치환하는 것이다.
4.4. CDF 표현하기¶
책의 저자가 작성한 thinkstats2,thinkplot을 사용했다.
first_cdf = thinkstats2.Cdf(first_baby.prglngth, label="prglngth")
thinkplot.Cdf(first_cdf)
thinkplot.Show(xlabel="weeks", ylabel="CDF")
CDF 그래프를 읽는 방법은 백분위수를 찾는 것이다. 위의 그림을 예로 들면 임신기간의 10%는 36주차보다 짧고 90%의 산모들은 41주보다 짧은 임신기간을 갖는다는 것을 알 수 있다.
데이터의 최빈값은 CDF 그래프에서 급격하거나 수직선으로 표현된다. 위 그림에서는 39주차에서 명확하게 표현되어 있다. 또한 30주보다 짧은 데이터는 매우 적기 때문에 평평하게 표현되어 있다.
이렇듯 CDF 그래프는 PMF 그래프보다 더 많은 정보를 보다 명확하게 시각화 할 수 있다는 장점이 있다.
4.5. CDF 비교하기¶
CDF는 서로 다른 데이터의 분포를 비교할때도 사용된다. 앞에서 했던 첫째 아이와 나머지 아이들의 임신기간의 차이를 다시 살펴보자.
first_cdf = thinkstats2.Cdf(first_baby.prglngth, label="first")
others_cdf = thinkstats2.Cdf(others.prglngth, label="others")
thinkplot.Cdf(first_cdf)
thinkplot.Cdf(others_cdf)
thinkplot.Show(xlabel="weeks", ylabel="CDF")
위 그림을 통해 분포 차이를 볼수 있다. 최빈값은 동일하나 첫째 아이가 나머지 아이들에 비해 임신기간이 조금 짧은 것을 볼 수 있다.
4.6. 백분위수 기반 통계량¶
중앙값(median) 또는 중위수는 어떤 주어진 값들을 크기의 순서대로 정렬했을 때 가장 중앙에 위치하는 값을 의미한다. 예를 들어 1, 2, 100의 세 값이 있을 때, 2가 가장 중앙에 있기 때문에 2가 중앙값이다.
사분위수 범위(Interquartile range, IQR)는 자료를 작은 값부터 큰 값까지 순서대로 나열한 후 4등분 하였을 때 각 점에 해당하는 값을 의미한다. 따라서 전체 자료를 100%로 보았을 때, 25% 지점에 해당하는 값은 제1사분위수가 되는 것이다.
- 25% = 제1사분위수 (q1) , the first quartile
- 50% = 제2사분위수(중앙값) (q2), the second quartile
- 75% 제 3사분위수 (q3), the third quartile
여기서 50% 지점에 해당하는 제2사분위수는 자료의 중앙값(중위수, median)과 동일하다.
4.7. 난수(Random numbers)¶
난수란 무작위로 만들어진 수열을 가리킨다. 여기서 무작위란 다음에 나올 수를 절대 예측할 수 없다는 것을 말한다.
weight_cdf = thinkstats2.Cdf(live_baby.totalwgt_lb, label="weight")
sample = np.random.choice(live_baby.totalwgt_lb, 100, replace=True)
ranks = [weight_cdf.PercentileRank(x) for x in sample]
rank_cdf = thinkstats2.Cdf(ranks)
thinkplot.Cdf(rank_cdf)
thinkplot.Show(xlabel="percentile rank", ylabel="CDF")
x = np.random.random(1000)
thinkplot.Pmf(thinkstats2.Pmf(x), linewidth=0.1)
thinkplot.Config(xlabel="Random variate", ylabel="PMF", ylim=(0, 0.05))
위 그래프에서 볼 수 있듯이 데이터들이 균등하게 분포한다는 것을 알 수 있다.
thinkplot.Cdf(thinkstats2.Cdf(x))
thinkplot.Config(xlabel="Random variate", ylabel="CDF")
위 CDF 그래프에서 확인할 수 있듯이, 매우 선형임으로 데이터셋이 균등분포라고 할 수 있다.
4.10. 용어사전¶
- 백분위순위(percentile rank): 분포에서 주어진 값과 동일하거나 적은값의 퍼센티지.
- 백분위수(percentile): 주어진 백분위 순위와 연관된 값.
- 누적분포함수(cumulative distribution function, CDF): 데이터를 누적 확률값으로 치환하는 방법. CDF(x)는 x와 동일하거나 작은 표본비율이다.
- 역CDF(inverse CDF): 누적확률(p)에서 해당 값으로 치환하는 함수.
- 중위값(median): 중심 경향성 측도로 사용되는 백분위의 50번째 데이터.
- 사분위범위(interquartile range): 퍼짐의측도로사용되는75번째와25 번째백분위수간차이.
- 복원(replacement): 표본추출 과정의 속성으로 복원추출(With replacement)은 동일한 값이한번이상 추출될수있다는 의미다. 그리고 비복원추출(Without replacement)은 값이 한번 추출되면 데이터셋에서 제거된다는 의미가 된다.
5. 분포 모델(Modeling distributions)¶
지금까지 배운 분포는 경험적 관측치에 기반한 데이터으로 구성된 경험적 분포(empirical distributions) 이다.
이런 경험적 분포를 수학적인 함수로 단순화하는 작업을 모델(model)이라 하며 모델을 통해 만들어지는 분포를 해석분포(analytic distribution) 라고 한다.
5.1. 지수분포(Exponential distribution)¶
가장 단순한 분포 유형인 지수 분포의 CDF는 다음과 같이 표현한다.
$$CDF(x) = 1 - e^{-\lambda x}$$
여기서 모수 $\lambda$는 분포의 모양을 나타낸다.
thinkplot.PrePlot(3)
for lam in [2.0, 1, 0.5]:
xs, ps = thinkstats2.RenderExpoCdf(lam, 0, 3.0, 50)
label = r"$\lambda=%g$" % lam
thinkplot.Plot(xs, ps, label=label)
thinkplot.Config(title="Exponential CDF", xlabel="x", ylabel="CDF", loc="lower right")
예시로 1997년 12월 18일 호주 브리즈번에서 태어난 아이들의 출생 시간에 대한 데이터셋 babyboom.dat을 사용한다.
def ReadBabyBoom(filename="../data/babyboom.dat"):
var_info = [
("time", 1, 8, int),
("sex", 9, 16, int),
("weight_g", 17, 24, int),
("minutes", 25, 32, int),
]
columns = ["name", "start", "end", "type"]
variables = pd.DataFrame(var_info, columns=columns)
variables.end += 1
dct = thinkstats2.FixedWidthVariables(variables, index_base=1)
df = dct.ReadFixedWidth(filename, skiprows=59)
return df
diffs = ReadBabyBoom().minutes.diff()
cdf = thinkstats2.Cdf(diffs, label="actual")
thinkplot.Cdf(cdf)
thinkplot.Config(xlabel="Time between births (minutes)", ylabel="CDF")
위 그림은 데이터셋의 출생시간간의 간격의 분포에 대해 시각화한 것이다.
보안적인 CDF(complementary CDF, 이하 CCDF) 는 위의 모델에 $Log$를 취하는 것이다. 공식으로 다음과 같이 표현 한다.
$$\log y \approx - \lambda x$$
따라서 로그를 취한 CCDF는 기울기가 $-\lambda$인 직선이 된다. 시각화를 통해 살펴보자.
thinkplot.Cdf(cdf, complement=True)
thinkplot.Config(
xlabel="Time between births (minutes)",
ylabel="CCDF",
yscale="log",
loc="upper right",
)
위 그래프는 명백한 직선이 아니라는 것을 알 수 있다. 따라서 이 데이터셋에 대해 지수분포는 완벽한 모델이 아니라는 것이 역으로 증명된다. 다시말해 아이의 "출생은 균등하게 발생할 것"이라는 가정은 틀렸다고 볼 수 있다.
5.2. 정규분포(normal distribution)¶
다른 말로 가우스 분포(Gaussian distribution) 라고 한다. 정규분포는 다양한 현상에 근사 모델로 사용될 수 있어 자주 사용된다.
정규분포는 모수 평균 $µ$, 표준편차 $σ$로 특성화된다. 그리고 평균이 0이고 표준편차가 1인 정규분포를 표준 정규 분포(standard normal distribution)라고 한다. 수식으로는 다음과 같이 표현된다.
$$ CDF(z) = \frac{1}{\sqrt{2 \pi}} \int^{z}_{-\infty} e^{\frac{-t^2}{2}}dt $$
모수의 범위에 따른 정규 분포 CDF를 시각화하면 다음과 같다
thinkplot.PrePlot(3)
mus = [1.0, 2.0, 3.0]
sigmas = [0.5, 0.4, 0.3]
for mu, sigma in zip(mus, sigmas):
xs, ps = thinkstats2.RenderNormalCdf(mu=mu, sigma=sigma, low=-1.0, high=4.0)
label = r"$\mu=%g$, $\sigma=%g$" % (mu, sigma)
thinkplot.Plot(xs, ps, label=label)
thinkplot.Config(title="Normal CDF", xlabel="x", ylabel="CDF", loc="upper left")
5.3. 정규확률그림(Normal probability plot)¶
n = 1000
thinkplot.PrePlot(3)
mus = [0, 1, 5]
sigmas = [1, 1, 2]
for mu, sigma in zip(mus, sigmas):
sample = np.random.normal(mu, sigma, n)
xs, ys = thinkstats2.NormalProbability(sample)
label = "$\mu=%d$, $\sigma=%d$" % (mu, sigma)
thinkplot.Plot(xs, ys, label=label)
thinkplot.Config(
title="Normal probability plot",
xlabel="standard normal sample",
ylabel="sample values",
)
위의 그림은 정규 분포에 대한 정규확률그림이다. 선들이 근사적인 직선이고 평균에 있는 값보다 벗어난 값을 꼬리에 갖는다.
이제 출생 체중 데이터에 대해 시도해 보자. 모델은 회색선으로 표현하고 실제 데이터는 파란색으로 그린다.
weights = live_baby.totalwgt_lb.dropna()
full_term = df[df.prglngth >= 37]
term_weights = full_term.totalwgt_lb.dropna()
mean, var = thinkstats2.TrimmedMeanVar(weights, p=0.01)
std = np.sqrt(var)
xs = [-4, 4]
fxs, fys = thinkstats2.FitLine(xs, mean, std)
thinkplot.Plot(fxs, fys, linewidth=4, color="0.8")
thinkplot.PrePlot(2)
xs, ys = thinkstats2.NormalProbability(weights)
thinkplot.Plot(xs, ys, label="all live")
xs, ys = thinkstats2.NormalProbability(term_weights)
thinkplot.Plot(xs, ys, label="full term")
thinkplot.Config(
title="Normal probability plot",
xlabel="Standard deviations from mean",
ylabel="Birth weight (lbs)",
)
위 그래프는 전체 정상 출생과 만삭(임신기간이36주이상)에 대한 결과를 보여준다. 두 곡선 모두 평균 근처에서 모형과 일치되고 꼬리 부분에서 차이가 난다. 무거운 아이는 모델로 예측한 것 보다 무겁고 가벼운 아이는 예측보다 더 가볍다.
만삭 데이터만 선택하면 가벼운 아이들의 데이터가 제거되고 그래프의 아래쪽 꼬리가 모델에 좀 더 잘 일치되는 것을 확인 할 수 있다.
5.4. 로그 정규분포(Lognormal distribution)¶
정규분포에 로그를 취하면 로그 정규분포가 된다.
$$CDF_{lognormal}(x) = CDF_{normal}(logx)$$
로그 정규분포의 예시로 미국에서 조사한 성인의 체중 데이터(2008년에 조사한 398,484명의 데이터)를 사용한다.
5.5. 파레토 분포(Pareto distribution)¶
파레토 분포(Pareto distribution)는 경제학자 빌프레도 파레토(Vilfredo Pareto) 이름에서 나왔는데 이것을 사용해 도시와 마을의 크기등을 나타냈다.
$$CDF(x) = 1 - \frac{x}{x_m}^{- \alpha}$$
5.6. 난수 생성하기¶
주어진 분포함수 $p = CDF(x)$로 난수를 생성한다.
import random
def expovariate(lam):
p = random.random()
x = -np.log(1 - p) / lam
return x
t = [expovariate(lam=2) for _ in range(1000)]
cdf = thinkstats2.Cdf(t)
thinkplot.Cdf(cdf, complement=True)
thinkplot.Config(xlabel="Exponential variate", ylabel="CCDF", yscale="log")
5.7. 왜 모델이 필요한가?¶
모델은 일종의 데이터 압축이다. 모델이 데이터셋에 잘 적합 될 때, 적은 모수 집합으로 대량의 데이터를 요약할 수 있는 것이다.
5.8. 연습문제¶
5.9. 용어 사전¶
- 경험적분포(empirical distribution): 표본값의 분포
- 해석분포(analytic distribution): CDF가 해석함수인 분포.
- 도착간격시간(interarrival time): 두사건 사이 경과 시간.
- 보완적 CDF(complementary CDF): 값x에서 x를 넘는 값 비율로 매칭하는 함수, $1 − CDF(x)$.
- 표준정규분포(standard normal distribution): 평균 0과 표준편차1을 갖는 정규분포.
- 정규확률그림(normal probability plot): 표본값과 표준정규분포에서 나온 난수를 대비해 시각화한 그림
6. 확률밀도함수¶
6.1. PDF¶
CDF 미분을 확률밀도함수(probability density function, 이하 PDF) 라고한다.
지수분포의 PDF는 다음과 같다.¶
$$PDF_{expo}(x) = \lambda e^{-\lambda x}$$
정규분포의 PDF는 다음과 같다.¶
$$PDF_{normal}(x) = \frac{1}{\sigma \sqrt{2\pi}} exp [ - \frac{1}{2}(\frac{x-\mu}{\sigma})^2]$$
6.2. 핵밀도추정¶
핵밀도추정(Kernel density estimation, 이하 KDE) 은 표본을 받아 데이터에 적합한 평활 PDF를 찾는 알고리즘이다. KDE는 다음과 같은 유용성이 있다.
- 시각화(Visualization): 데이터셋의 탐색 단계에서는 대체로 KDE가 분포를 가장 잘 시각화한다.
- 보간(Interpolation): KDE를 사용해서 표본에 없는값에대해 밀도를 보간할 수 있다.
- 모의실험(Simulation): 만약 표본 크기가 작다면 KDE를 사용해 표본 분포를 평활할 수 있다.
보간이란 알려진 지점의 값 사이(중간)에 위치한 값을 알려진 값으로부터 추정하는 것을 말한다.
6.3. 분포 프레임워크¶
지금까지 PMF, CDF, PDF를 살펴봤다, 잠시 복습시간을 갖는다. 다음 그림에 각각의 함수가 어떻게 연관 되는지 나타나 있다.
6.8. 왜도¶
왜도(Skewness)는 분포의 형태를 기술한다. 만약 분포가 중심경향성 주변에서 대칭이라면 기울어 지지 않았다. 만약 값들이 오른쪽으로 좀 더 뻗어져있다면“오른쪽으로기울어져(right skewed), 만약 값들이 왼쪽으로 치우쳐있다면,“왼쪽으로 기울어져(left skewed)”이다.
6.9. 연습문제¶
- 소득의 분포는 기울어져 있기로 유명하다. 소득 데이터셋은 http://www.census.gov/hhes/www/cpstables/032013/hhinc/toc.htm 에서 다운로드 받을수 있다. 얼마나 치우쳐 있는지 측정하라.
income_df = pd.read_csv("../data/hinc06(1).csv")
def InterpolateSample(df, log_upper=6.0):
df["log_upper"] = np.log10(df.income)
df["log_lower"] = df.log_upper.shift(1)
df.loc[0, "log_lower"] = 3.0
df.loc[41, "log_upper"] = log_upper
arrays = []
for _, row in df.iterrows():
vals = np.linspace(row.log_lower, row.log_upper, int(row.freq))
arrays.append(vals)
log_sample = np.concatenate(arrays)
return log_sample
log_sample = InterpolateSample(income_df, log_upper=6.0)
log_cdf = thinkstats2.Cdf(log_sample)
thinkplot.Cdf(log_cdf)
thinkplot.Config(xlabel="Household income (log $)", ylabel="CDF")
sample = np.power(10, log_sample)
cdf = thinkstats2.Cdf(sample)
thinkplot.Cdf(cdf)
thinkplot.Config(xlabel="Household income ($)", ylabel="CDF")
def RawMoment(xs, k):
return sum(x**k for x in xs) / len(xs)
def Median(xs):
cdf = thinkstats2.Cdf(xs)
return cdf.Value(0.5)
def Mean(xs):
return RawMoment(xs, 1)
Mean(sample), Median(sample)
평균과 중위값은 각각 74278, 51226임을 알 수 있다.
def CentralMoment(xs, k):
mean = RawMoment(xs, 1)
return sum((x - mean) ** k for x in xs) / len(xs)
def StandardizedMoment(xs, k):
var = CentralMoment(xs, 2)
std = np.sqrt(var)
return CentralMoment(xs, k) / std**k
def Skewness(xs):
return StandardizedMoment(xs, 3)
def PearsonMedianSkewness(xs):
median = Median(xs)
mean = RawMoment(xs, 1)
var = CentralMoment(xs, 2)
std = np.sqrt(var)
gp = 3 * (mean - median) / std
return gp
Skewness(sample), PearsonMedianSkewness(sample)
왜도는 4.9로 데이터셋이 왼쪽으로 치우쳐있다는 것을 알 수 있다. (정규분포의 왜도는 0이다.) 그러나 이 계산는 가장 높은 소득을 백만달러라고 가정한 결과이다. 다시말해 최대값에 대한 정보가 없으면 계산된 왜도는 믿을 수가 없다.
6.10. 용어사전¶
- 확률밀도함수(Probability density function, PDF): 연속CDF 미분으로 값을 확률밀도로 치환하는 함수
- 확률밀도(Probability density): 확률을 만들기 위해 범위값에 대해 적분 할 수 있는양.
- 핵밀도추정(Kernel density estimation, KDE): 표본에기반해서PDF를 추정하는알고리즘.
- 이산화(discretize): 연속 함수 혹은 이산 함수를 가진 분포를 근사함. 평활(smoothing)의반대.
- 원적률(raw moment): 거듭제곱되는 데이터 합계에 기반한 통계량
- 중심적률(central moment): 평균에서 편차 거듭제곱에 기반한 통계량.
- 표준적률(standardized moment): 단위가 없는 적률 비율.
- 왜도(skewness): 분포가 얼마나 비대칭인지 나타내는 측도.
- 표본왜도(sample skewness): 분포 왜도를 정량화 하는데 사용되는 적률 기반 통계량.
- 피어슨중위수왜도계수(Pearson’s median skewness coefficient): 중위수, 평균, 표준편차에 기반한 분포 왜도를 정량화하는데 사용되는 통계량.
- 강건성(robust): 특이값에도 휘둘리지 않는 통계량을 강건하다 표현한다.
7. 변수간 관계¶
변수들에는 관계가 있다. 예를 들면 키와 체중은 관계가 있다; 키가 큰 사람이 체중도 더 나가는 경향이 있다. 물론 키가 작고 뚱뚱한 사람도 있을 수 있지만 키에 대한 정보를 알고 있으면 체중을 추정하는데 더 유리하다.
7.1. 산점도¶
변수간의 관계(relationship)을 확인하는 가장 간단한 방법은 산점도(scatter plot)를 그리는 것이다.
import pandas as pd
df_brfss = pd.read_csv("../data/brfss.csv")
sample = thinkstats2.SampleRows(
df_brfss, 5000
) # SampleRows함수는 데이터에서 일부를 골라낸다.
heights, weights = sample.htm3, sample.wtkg2
thinkplot.Scatter(heights, weights)
thinkplot.Show(xlabel="Height (cm)", ylabel="Weight (kg)", axis=[140, 210, 20, 200])
위 그림은 데이터의 유실을 명백히 보여준다. 실제 키에 대한 데이터는 인치로 측정되었는데, 수치를 반올림해서 센치미터로 변환했기 때문이다. 이런 데이터 유실을 보안하기 위해 임의로 잡음을 추가하는 지터링(jittering) 을 한다.
def Jitter(values, jitter=1.5):
n = len(values)
return np.random.uniform(-jitter, +jitter, n) + values
thinkplot.Scatter(Jitter(heights), Jitter(weights))
thinkplot.Show(xlabel="Height (cm)", ylabel="Weight (kg)", axis=[140, 210, 20, 200])
지터링(jittering)을 통해서 반올림으로 인한 시각적 효과를 줄이고, 변수의 관계 형태를 좀 더 명확해진다. 하지만, 일반적으로 지터링은 시각화 목적으로만 사용해야 하고 분석에 지터링된 데이터 사용은 피해야 한다.
7.2. 변수의 관계를 특징짓기¶
변수의 관계에 대한 더 많은 통찰을 얻기 위해 변수를 구간화(binning)하고 다른 변수 백분위수를 시각화해보자. 넘파이(NumPy)와 판다스를 이용해 데이터를 구간으로 나눈다.
df_brfss = df_brfss.dropna(subset=["htm3", "wtkg2"])
bins = np.arange(135, 210, 5)
indices = np.digitize(df_brfss.htm3, bins)
groups = df_brfss.groupby(indices)
heights = [group.htm3.mean() for i, group in groups]
cdfs = [thinkstats2.Cdf(group.wtkg2) for i, group in groups]
for percent in [75, 50, 25]:
weights = [cdf.Percentile(percent) for cdf in cdfs]
label = "%dth" % percent
thinkplot.Plot(heights, weights, label=label)
위 그림 결과를 보면 140cm에서 200cm 사이 두 변수간의 관계는 대략작인 선형 관계이다. 그리고 이 범위에 99% 이상의 데이터가 몰려있다. 그래서, 극단값들은 신경쓸 필요는 없다.
7.3. 상관(Correlation)¶
상관(correlation) 은 두 변수간 관계를 정량화하는데 사용되는 통계량이다.
7.4. 공분산(Covariance)¶
공분산(Covariance) 은 함께 변화하는 두 변수의 경향성이다. 만약 두 변수 X, Y가 있다면 평균에 대한 편차는 다음과 같다.
7.5. 피어슨 상관(Pearson’s correlation)¶
공분산이 몇몇 연산에서는 유용하지만 해석하기 어렵기 때문에 요약 통계량으로 많이 사용되지 않는다. 반면에 피어슨상관(Pearson's correlation)은 계산하기 쉽고 해석하기도 쉬워 많이 선형 상관을 판별하는데 많이 사용된다.
7.6. 비선형관계(Nonlinear relationships)¶
상관계수가 0이어도 변수간에 관계가 존재하는 것을 비선형관계라고 한다. 다음 그림을 통해
다양한 범위의 상관을 갖는 데이터들을 확인할 수 있다.
위 그림에서 세번째 행이 비선형 관계의 데이터셋이다.
7.7. 스피어만 순위 상관(Spearman’s rank correlation)¶
피어슨 상관은 변수의 관계가 선형이고 정규분포를 따른다면 잘 작동한다. 하지만 데이터에 특이값이 있는 경우에는 스피어만 순위 상관을 대안으로 사용한다.
7.8. 상관과 인과관계(Correlation and causation)¶
만약 변수 A와 B가 상관이 있다면 다음과 같은 세가지 설명이 가능하다.
- A가 B의 발생원인이다.
- B가 A의 발생원인이다.
- 다른 변수가 A와 B의 발생 원인이다.
위와 같은 설명을 인과관계라고 부른다. 인과관계를 증명하기 위해 다음을 수행한다.
- 사건의 순서를 사용한다. 만약 A가 B 보다 먼저 일어났다면 A가 B를 발생시키는 원인이라고 추론할 수 있을 것이다. 그러나 다른 무언가 A와 B를 발생시키는 원인이 된다는 가능성을 배제할 수 없다.
- 확률을 사용한다. 표본을 임의로 두 집단으로 나눠 변수의 평균을 계산하는 무작위 대조군 시험을 통해 인과 관계를 확인할 수 있다.
7.9. 연습문제¶
- NSFG 데이터를 사용해 출생 체중과 산모 연령의 산점도를 그리고, 출생 체중과 산모 연령백분위수를 시각화하라. 피어슨상관과 스피어만 상관을 계산해 두 변수사이의 관계를 설명하라.
live = live_baby.dropna(subset=["agepreg", "totalwgt_lb"])
ages = live.agepreg
weights = live.totalwgt_lb
thinkplot.Scatter(ages, weights)
thinkplot.Config(
xlabel="Age (years)",
ylabel="Birth weight (lbs)",
xlim=[10, 45],
ylim=[0, 15],
legend=False,
)
from scipy.stats import pearsonr, spearmanr
print(f"피어슨 상관 계수: {pearsonr(ages, weights)[0]}")
print(f"스피어만 상관 계수: {spearmanr(ages, weights)[0]}")
def BinnedPercentiles(df):
bins = np.arange(10, 48, 3)
indices = np.digitize(df.agepreg, bins)
groups = df.groupby(indices)
ages = [group.agepreg.mean() for i, group in groups][1:-1]
cdfs = [thinkstats2.Cdf(group.totalwgt_lb) for i, group in groups][1:-1]
thinkplot.PrePlot(3)
for percent in [75, 50, 25]:
weights = [cdf.Percentile(percent) for cdf in cdfs]
label = "%dth" % percent
thinkplot.Plot(ages, weights, label=label)
thinkplot.Config(
xlabel="Mother's age (years)",
ylabel="Birth weight (lbs)",
xlim=[14, 45],
legend=True,
)
BinnedPercentiles(live)
7.10. 용어사전¶
- 산점도(scatter plot): 두 변수 사이 관계를 데이터 각 행마다 점을 찍어 보임으로써 시각화.
- 지터(jitter): 시각화를 위해 데이터에 추가되는 잡음
- 포화(saturation): 다수의 점이 겹쳐져 발생하는 데이터의 손실.
- 상관(correlation): 두 변수 사이 관계를 측정 하는 통계량.
- 표준화(standardize): 변수 집합값을 변환해 평균 0, 분산을 1로 만드는 과정.
- 표준점수(standard score): 통계학적으로 정규분포를 만들고 개개의 경우가 표준편차상에 어떤 위치를 차지하는지를 보여주는 차원없는 수치이다. 표준값, Z값(Z-value), Z 점수(Z score)이라고도 한다.
- 공분산(covariance): 함께 변화하는 두변수 경향성을 나타내는 측도.
- 순위(rank): 인덱스로 정렬된 목록에 요소의 순서.
- 무작위대조군시험: 실험 대상이 무작위 집단으로 나눠 집단간 차이를 살펴보는 실험.
- 처리군집단(treatment group): 대조군시험에서 일종의 개입(intervention)을 받는 집단.
- 대조군집단(control group): 대조군시험에서 처리를 받지 않거나 이미 효과가 알려진 처리를 받는 집단.
- 자연실험(natural experiment): 대상이 집단으로 자연적으로 구분되는 현상을 이용한 실험계획. 여기서 집단은 근사적으로 무작위가 된다.
8. 추정¶
8.1. 통계적 추정¶
추정(estimation) 은 표본을 통해 모집단 특성이 어떠한가에 대해 추측하는 것이다. 그리고 추정에 사용한 통계량(표본평균)을 추정량(estimator) 이라 한다.
어느 추정량이 가장 좋은가는 상황에 따라 달라진다. 예를들어 특이값의 유무와 목적이 무엇인가에 달려있다. 만약 특이값이 없다면 누적평균제곱오차(mean squared error, 이하 MSE) 를 최소화하는 것이 좋은 방법이다.
그러나 주사위 세개를 던져서 합계를 예측하는 문제에서 MSE 값은 10.5로 계산되지만 좋지못한 추정값이다. 왜냐하면 주사위 세개의 합은 결코 10.5 가 되지 못하기 때문이다. 이 경우 10 또는 11이 최선의 값이 되며 이것을 최대우도추정량(maximum likelihood estimator, 이하 MLE) 이라 부른다.
8.2. 모분산 추정¶
어떤 공장에서 일정 크기의 볼트를 만든다고 해보자. 제조 공정에서 볼트의 지름이 일정해야 품질 기준을 만족할 수 있다. 볼트의 지름이 10mm 여야하는데, 평균적으로는 10mm 이지만, 분산이 존재해 불량품이 일정 부분 존재한다. 이런 상황에서 볼트 지름의 분산을 확인하는 방법으로 모분산에 대한 추론을 할 수 있다.
출처: https://3months.tistory.com/480 [Deep Play]
8.3. 표본분포(Sampling distributions)¶
크기 n의 확률 표본(random sample)의 확률 변수(random variable)의 분포(distribution)이다. 또한 확률선택(random selection) 에 의해 발생 되는 추정값의 변동을 표집오차(sampling error)라고한다.
def SimulateSample(mu=90, sigma=7.5, n=9, iters=1000):
xbars = []
for j in range(iters):
xs = np.random.normal(mu, sigma, n)
xbar = np.mean(xs)
xbars.append(xbar)
return xbars
def VertLine(x, y=1):
thinkplot.Plot([x, x], [0, y], color="0.8", linewidth=3)
xbars = SimulateSample()
cdf = thinkstats2.Cdf(xbars)
thinkplot.Cdf(cdf)
ci = cdf.Percentile(5), cdf.Percentile(95)
VertLine(ci[0])
VertLine(ci[1])
thinkplot.Config(xlabel="Sample mean", ylabel="CDF")
표집 분포를 요약하는 두가지 방법이 있다.
- 표준오차(Standard error SE) 는 평균적으로 추정값이 얼마나 차이가 있을지 측정하는 측도다. 각각의 모의실험 마다, 오차($\hat x - \mu$)를 계산해 누적 평균제곱오차에제곱근(root mean squared error, RMSE)을 씌워 계산한다.
- 신뢰구간(confidence interval, CI) 은 표집 분포를 주어진 비율로 포함할 범위가 된다. 예를들어, 신뢰구간 90%는 5번째부터 95번째 백분위수가 된다.
8.4. 선택편향(Sampling bias)¶
'표본 편향'이라고도 한다. 상당히 많은 자료들을 검토하였으나 그 자료를 선택하거나 해석함에 있어 중요한 측면을 간과함으로써 잘못된 결론에 도달하게 만든다.
8.5. 지수분포¶
통계학에서 지수분포(exponential distribution) 는 연속 확률 분포의 일종이다. 사건이 서로 독립적일 때, 일정 시간동안 발생하는 사건의 횟수가 푸아송 분포를 따른다면, 다음 사건이 일어날 때까지 대기 시간은 지수분포를 따른다.
8.6 연습문제¶
- 모수(λ)가 2인 지수 분포에서 표본(n) 10개를 추출한다. 1000번의 실험에서 추정값(L)의표본 분포를 시각화한다. 그리고 추정값의 표준오차와 90% 신뢰구간을 계산하라.
def RMSE(estimates, actual):
e2 = [(estimate - actual) ** 2 for estimate in estimates]
mse = np.mean(e2)
return np.sqrt(mse)
def SimulateSample(lam=2, n=10, iters=1000):
def VertLine(x, y=1):
thinkplot.Plot([x, x], [0, y], color="0.8", linewidth=3)
estimates = []
for _ in range(iters):
xs = np.random.exponential(1.0 / lam, n)
lamhat = 1.0 / np.mean(xs)
estimates.append(lamhat)
stderr = RMSE(estimates, lam)
print("standard error", stderr)
cdf = thinkstats2.Cdf(estimates)
ci = cdf.Percentile(5), cdf.Percentile(95)
print("confidence interval", ci)
VertLine(ci[0])
VertLine(ci[1])
# plot the CDF
thinkplot.Cdf(cdf)
thinkplot.Config(xlabel="estimate", ylabel="CDF", title="Sampling distribution")
return stderr
SimulateSample()
8.7. 용어사전¶
- 추정(estimation): 표본에서 분포를 추정하는 과정.
- 추정량(estimator): 모수를 추정하는데 사용되는 통계량.
- 누적평균제곱오차(mean squared error, MSE): 추정오차측도.
- 제곱근평균제곱오차(root mean squared error, RMSE): MSE 제곱근으로 오차규모를 좀 더일반적으로 표현.
- 최우추정량(maximum likelihood estimator, MLE): 가장 올바를것 같은 추정값을 계산하는추정량.
- (추정량의) 편의(bias of an estimator): 반복되는 실험을 평균낼때, 실제 모수값 보다 높은 혹은 낮은 추정량의 경향성.
- 표집오차(sampling error): 우연에의한변동과한정된표본크기때문에 추정값에생기는오차.
- 표집편의(sampling bias): 모집단을대표하지못하는표집과정때문에 추정값에생기는오차.
- 측정오차(measurement error): 데이터를수집하고기록하는부정확성때문에 추정값에생기는오차.
- 표집분포(sampling distribution): 만약실험을많이반복한다면생기는 통계량분포.
- 표본오차(standard error): 추정값의RMSE로표집오차(하지만다른 오차는제외) 때문에생기는변동성을정량화한다.
- 신뢰구간(confidence interval): 만약실험을많이반복하면추정량예상 범위를 표현하는구간.
9. 가설 검정¶
9.1. 전통적 가설 검정¶
- 첫번째 단계는 검정통계량(test statistic) 을 정해 효과 크기를 정량화 한다. 2. 두번째 단계는 귀무가설(null hypothesis) 을 정의 하는 것으로 명백한 차이가 없다고 가정한다.
- 세번째 단계는 만약 귀무가설이 사실이라면 나타날 확률인 p-value를 계산 하는 것이다.
- 마지막 단계는 결과를 해석 하는 것이다. 만약 p-value가 충분히 낮다면 통계적 유의성(statistically significant) 이 있다고 본다. 다시말해 모집단에서도 동일한 현상이 나타난다고 추론할 수 있다.
9.2. 가설 검정¶
저자는 thinkstats2 모듈에 전통적 가설 검정 구조를 표현하는 HypothesisTest 클래스를 정의해 두었다. 우리는 쓰기만 하면 된다.
9.3. 평균 차이 검정¶
예시로 NSFG 데이터셋에서 첫 번째 아이에 대한 평균 임신기간과 평균 출산 체중이 나머지 아이들과 통계적으로 유의미한 차이가 있는지 살펴보자.
class DiffMeansPermute(thinkstats2.HypothesisTest):
def TestStatistic(self, data):
group1, group2 = data
test_stat = abs(group1.mean() - group2.mean())
return test_stat
def MakeModel(self):
group1, group2 = self.data
self.n, self.m = len(group1), len(group2)
self.pool = np.hstack((group1, group2))
def RunModel(self):
np.random.shuffle(self.pool)
data = self.pool[: self.n], self.pool[self.n :]
return data
data = first_baby.prglngth.values, others.prglngth.values
ht = DiffMeansPermute(data)
pvalue = ht.PValue()
print(f"p-value 값은 {pvalue} 이다.")
결과는 약 0.19로 약 19% 정도 관측된 효과 차이가 예상된다. 다시말해 효과가 통계적으로 유의 하지 않다고 볼 수 있다.
이제 관측 효과 크기를 나타내는 회색선과 검정 통계량 분포를 시각화하자.
ht.PlotCdf()
thinkplot.Config(xlabel="test statistic", ylabel="CDF")
CDF는 0.83에서 관측차이와 교차한다.
9.4. 다른 검정 통계량¶
앞선 검정은 양쪽을 사용하기 때문에양측(two-sided) 검증이라고한다. TestStatistic가차이에 절대값을 취하지 않는다는 것이다. 만약 첫번째 아이가 늦게 출생하는 경향이 있다고 가정한다면 이런 유형의 검증을 단측(one-sided)검증 이라고한다. 왜냐하면, 차이 분포의단측면만 고려하기때문이다.
class DiffStdPermute(DiffMeansPermute):
def TestStatistic(self, data):
group1, group2 = data
test_stat = group1.std() - group2.std()
return test_stat
ht = DiffStdPermute(data)
pvalue = ht.PValue()
print(f"p-value 값은 {pvalue} 이다.")
단측가설, 첫째 아이가 늦게 태어난다는 것이 양측가설보다 좀 더 구체적이다. 그래서 p-value가 더 작게 계산된다. 하지만 그럼에도 그 차이는 통계적으로 유의해보이지 않는다.
9.5. 상관검정¶
NSFG 데이터 셋에서 산모의 나이와 출생 체중 사이 상관관계는 약 0.07 이다. 즉, 나이가 많은 산모 아이가 체중이 더 나가는 것처럼 보인다. 이것은 과연 우연일까? 귀무가설로 산모 연령과 신생아 체중사이에 상관관계는 없다고 가정해 검정해보자.
class CorrelationPermute(thinkstats2.HypothesisTest):
def TestStatistic(self, data):
xs, ys = data
test_stat = abs(thinkstats2.Corr(xs, ys))
return test_stat
def RunModel(self):
xs, ys = self.data
xs = np.random.permutation(xs)
return xs, ys
data = live_baby.agepreg.values, live_baby.totalwgt_lb.values
ht = CorrelationPermute(data)
pvalue = ht.PValue()
print(f"p-value 값은 {pvalue} 이다.")
계산된 p-value는 0 으로 귀무가설을 기각한다. 따라서 산모의 연령과 신생아의 체중에는 관계가 있다고 볼 수 있다.
9.6. 모비율 검정¶
모비율의 검정은 모비율인 p가 이럴것이다. 라고 귀무가설을 세우고 검정 통계량을 계산하는 것이다.
9.7. 카이제곱 검정¶
카이제곱 검정은 모비율 검정에서 일반적으로 사용되는 통계량이다. 수식으로는 다음과 같이 표현한다.
$$ x^2 = \sum_{i} \frac{(O_i - E_i)^2}{E_i} $$
여기서 $O_i$는 관측빈도 $E_i$는 기대 빈도이다.
9.8. 첫번째 아이의 임신기간¶
35 - 43주차에서 첫번째 아이와 나머지 아이의 카이제곱 검정을 수행한다.
class PregLengthTest(thinkstats2.HypothesisTest):
def MakeModel(self):
firsts, others = self.data
self.n = len(firsts)
self.pool = np.hstack((firsts, others))
pmf = thinkstats2.Pmf(self.pool)
self.values = range(35, 44)
self.expected_probs = np.array(pmf.Probs(self.values))
def RunModel(self):
np.random.shuffle(self.pool)
data = self.pool[: self.n], self.pool[self.n :]
return data
def TestStatistic(self, data):
firsts, others = data
stat = self.ChiSquared(firsts) + self.ChiSquared(others)
return stat
def ChiSquared(self, lengths):
hist = thinkstats2.Hist(lengths)
observed = np.array(hist.Freqs(self.values))
expected = self.expected_probs * len(lengths)
stat = sum((observed - expected) ** 2 / expected)
return stat
data = first_baby.prglngth.values, others.prglngth.values
ht = PregLengthTest(data)
p_value = ht.PValue()
print("p-value =", p_value)
print("actual =", ht.actual)
print("ts max =", ht.MaxTestStat())
NSFG 데이터셋의 카이제곱값 101.5 은 그 자체로 많은 것을 의미 하지는 않는다. 하지만 귀무가설 아래에서 1000회 반복했을때 가장 큰 검정 통계량은 25.7 다. 다시 말해 실제 카이 제곱 통계량이 귀무가설 아래에서 가능할 것 같지 않다고 결론 내릴 수 있다. 그래서 외관 효과는 통계적 유의성이 있다고 볼 수 있다. 동시에 이 예제는 카이제곱 검정한계를 보여준다. 바로 두 집단 사이에 차이가 있다는 것은 알 수 있지만 차이가 무엇인지에 대한 구체적인 것은 알 수 없다는 것이다.
9.9 오류¶
일반적으로 전통적 가설 검정에서 p-value가 0.05보다 낮다면 보통 통계적으로 유의미 하다고 판단한다. 그리고 그런 경우에도 다음 두 가지를 고려 한다.
- 만약 효과가 실제로 우연이라면 해당 효과를 잘못 판단할 확률은 얼마나될까? 이것을 이거짓양성률(false positive rate)이라 한다.
- 만약 효과가 진실이라면, 가설검정이 실패할 가능성은 얼마나될까? 이 확률은 거짓 음성률(false negative rate) 이라 한다.
9.10. 검정력(Power)¶
거짓양성률은 상대적으로 계산하기 쉽다. 분계점이 5% 이면, 거짓양성률은 5% 이다. 반면에 거짓음성률은 계산하기 더 어렵다. 왜냐면 실제 효과크기에 의존하고 정상적으로는 실제효과 크기를 모르기 때문이다. 한 가지 선택옵션은 가상효과크기(hypothetical effect size) 조건으로 거짓음성률을 계산하는 것이다.
정확한 양성률(correct positive rate)을 검정력(power)라고 부르며 종종 민감도(sensitivity)라 한다.
9.11. 반복(Replication)¶
현재 데이터에서 통계적 유의성이있는 모든 효과가 새로운 데이터셋에서도 반복되어야 한다.
9.12. 연습문제¶
- 표본크기가 증가하면 가설 검정력이 증가하고 반대로 표본크기가 줄어들면 검정력은 줄어든다. 표본의 크기가 줄어듬에 따라 p-value 가 어떻게 변화하는지 최소한의 표본크기는 얼마인지 알아보라.
def RunTests(live, iters=1000):
n = len(live)
firsts = live[live.birthord == 1]
others = live[live.birthord != 1]
# compare pregnancy lengths
data = firsts.prglngth.values, others.prglngth.values
ht = DiffMeansPermute(data)
p1 = ht.PValue(iters=iters)
data = (firsts.totalwgt_lb.dropna().values, others.totalwgt_lb.dropna().values)
ht = DiffMeansPermute(data)
p2 = ht.PValue(iters=iters)
# test correlation
live2 = live.dropna(subset=["agepreg", "totalwgt_lb"])
data = live2.agepreg.values, live2.totalwgt_lb.values
ht = CorrelationPermute(data)
p3 = ht.PValue(iters=iters)
# compare pregnancy lengths (chi-squared)
data = firsts.prglngth.values, others.prglngth.values
ht = PregLengthTest(data)
p4 = ht.PValue(iters=iters)
print("%d\t%0.2f\t%0.2f\t%0.2f\t%0.2f" % (n, p1, p2, p3, p4))
n = len(live_baby)
for _ in range(7):
sample = thinkstats2.SampleRows(live_baby, n)
RunTests(sample)
n //= 2
예상대로 큰 샘플에서 양성인 데이터를 제거하면 효과 크기가 줄어든다. 그러나 작은 샘플 크기에서도 패턴은 불규칙하게 일어난다.
9.13. 용어사전¶
- 가설검정(hypothesis testing): 외관효과(apparent effect)가 통계적 유의성이 있는지 결정하는 과정.
- 검정통계량(test statistic): 효과크기를 정량화하는데 사용되는 통계량.
- 귀무가설(null hypothesis): 외관효과가 우연에 의한 것이라는 가정에 기반한 시스템모델.
- p-값(p-value): 효과가 우연에 의해서 발생한 확률.
- 통계적유의성(statistically significant): 우연에 의해서 발생할 것 같지않다면, 효과는 통계적 유의성이 있다.
- 순열검정(permutation test): 관측 데이터셋에서 순열을 생성해서 p-값을 계산하는방법.
- 재표본검정(resampling test): 관측 데이터셋에서 복원으로 표본을 생성해서 p-값을 계산하는방법.
- 양측검정(two-sided test): “효과크기가 양으로 혹은 음으로 관측 효과보다 얼마나 큰가?”를 묻는 검정.
- 단측검정(one-sided test): “효과크기가 동일 부호로 관측 효과보다 얼마나 큰가?”를 묻는 검정.
- 카이-제곱검정(chi-squared test): 검정통계량으로 카이-제곱통계량을 사용하는검정.
- 거짓양성(false positive): 사실이 아닐때 효과가 사실이라는 결론.
- 거짓음성(false negative): 우연이 아닐때 효과가 우연 때문이라는 결론.
- 검정력(power): 대립 가설이 사실일때, 이를 사실로서 결정할 확률이다. 검정력이 90% 라고하면, 대립가설이 사실임에도 불구하고 귀무가설을 채택 할 확률은 10%이다.
10. 선형최소제곱(Linear least squares)¶
10.1. 최소제곱 적합(Lease squares fit)¶
변수간의 기울기를 측정하는 가장 흔한 방법이 선형최소제곱 적합(linearleast squares fit) 이다. 선형 적합은 변수 사이의 관계를 하나의 선으로 모델링한다. 따라서 소제곱 적합은 모델과 데이터 사이의 평균제곱오차를 최소화 하는 것이다.
모델과 데이터 사이의 차이를 잔차(residual) 라고 한다. 잔차에 제곱을 하는 이유는 다음과 같다.
- 제곱하게 되면 잔차의 부호를 제거하는 기능이 있다.
- 제곱은 큰 잔차에 더 많은 가중치를 주지만 가장 큰 잔차가 항상 주도적인 경우에는 그렇게 많은 가중치를 주지는 않는다.
- 만약 잔차가 상관관계가 없고 평균과 알수없는 분산을 가진 정규분포라면, 최소제곱적합은 또한 최대우도 추정량이다.
from thinkstats2 import Mean, MeanVar, Var, Cov
def LeastSquares(xs, ys):
meanx, varx = MeanVar(xs)
meany = Mean(ys)
slope = Cov(xs, ys, meanx, meany) / varx
inter = meany - slope * meanx
return inter, slope
def FitLine(xs, inter, slope):
fit_xs = np.sort(xs)
fit_ys = inter + slope * fit_xs
return fit_xs, fit_ys
live = live_baby.dropna(subset=["agepreg", "totalwgt_lb"])
ages = live.agepreg
weights = live.totalwgt_lb
inter, slope = LeastSquares(ages, weights)
fit_xs, fit_ys = FitLine(ages, inter, slope)
inter, slope
thinkplot.Scatter(ages, weights, color="blue", alpha=0.1, s=10)
thinkplot.Plot(fit_xs, fit_ys, color="red", linewidth=3)
thinkplot.Config(
xlabel="Mother's age (years)",
ylabel="Birth weight (lbs)",
axis=[10, 45, 0, 15],
legend=False,
)
선과 함께 출생 체중과 연령을 산점도로 보여준다. 관계가 선형인지, 적합선이관계를 나타내는 좋은 모형인지를 평가 하는데 시각화는 언제나 좋은 생각이다.
10.3. 잔차(Residuals)¶
잔차를 시각화 해보자.
def Residuals(xs, ys, inter, slope):
xs = np.asarray(xs)
ys = np.asarray(ys)
res = ys - (inter + slope * xs)
return res
live["residual"] = Residuals(ages, weights, inter, slope)
bins = np.arange(10, 48, 3)
indices = np.digitize(live.agepreg, bins)
groups = live.groupby(indices)
age_means = [group.agepreg.mean() for _, group in groups][1:-1]
cdfs = [thinkstats2.Cdf(group.residual) for _, group in groups][1:-1]
def PlotPercentiles(age_means, cdfs):
thinkplot.PrePlot(3)
for percent in [75, 50, 25]:
weight_percentiles = [cdf.Percentile(percent) for cdf in cdfs]
label = "%dth" % percent
thinkplot.Plot(age_means, weight_percentiles, label=label)
PlotPercentiles(age_means, cdfs)
thinkplot.Config(xlabel="Mother's age (years)", ylabel="Residual (lbs)", xlim=[10, 45])
이상적으로는 모든 연령집단에서 잔차가 평평한 직선으로 나타나야 한다. 그러나 실제 데이터에는 위의 그림과 같은 모습을 띄게 된다.
10.4. 추정(Estimation)¶
모수 slope 와 inter 는 표본에 기반한 추정값이다
10.5. 적합도(Goodness of fit)¶
선형 모델의 품질은 적합도(goodness of fit) 로 나타내며, 가장 간단하게 측정하는 방법은 잔차의 표준편차를 구하는 것이다.
적합도를 측정하는 또 다른방식은 $R^2$로 표기하고, R-제곱(R-squared) 라고 부르는 결정계수(coefficient of determination)다.
def CoefDetermination(ys, res):
return 1 - Var(res) / Var(ys)
inter, slope = LeastSquares(ages, weights)
res = Residuals(ages, weights, inter, slope)
r2 = CoefDetermination(weights, res)
r2
신생아의 체중과 산모의 나이에 대한 적합도 $R^2$ 값은 0.0047 이다. 그말은 산모 나이가 신생아의 체중을 약 0.5% 확률로 예측한다는 의미가 된다.
10.9. 용어사전¶
- 선형적합(linear fit): 변수사이 관계를 모델링하는 직선.
- 최소제곱적합(least squares fit): 잔차 제곱합을 최소화하는 데이터 모델.
- 잔차(residual): 실제값과 모델값의 차이.
- 적합도(goodness of fit): 모델이 데이터에 얼마나 잘 적합하는지에 대한 척도.
- 결정계수(coefficient of determination): 적합도를 계량화하는 통계량.
- 표집가중치(sampling weight): 표본이 모집단 어느부분을 대표하는지 나타내는 표본 관측과 연관된 값.
11. 회귀¶
회귀(regression)는 데이터를 모형에 적합하는데 사용되는 방법이다. 회귀분석의 목적은 종속변수(dependent variables)라고 불리는 변수집합과 설명변수(explanatory variables)라고불리는 또 다른 집합변수사이관계를 기술하는 것이다.
11.1. StatsModels¶
statsmodels 파이썬 패키지는 몇가지 형태의 회귀분석과 다른 통계 분석 기능을 제공한다.
11.2. 다중회귀(Multiple regression)¶
다중 회귀란 설명변수(독립변수)가 2 개 이상인 회귀모형을 분석대상으로 삼는다.
- 추가적인 독립변수를 도입함으로써 오차항의 값을 줄일 수 있다.
- 단순회귀분석의 단점을 극복 할 수 있다.
import statsmodels.formula.api as smf
formula = "totalwgt_lb ~ agepreg"
model = smf.ols(formula, data=live)
results = model.fit()
results.summary()
inter = results.params["Intercept"]
slope = results.params["agepreg"]
inter, slope
11.3. 비선형관계(Nonlinear relationships)¶
위의 agepreg 변수가 비선형관계일지 모른다는 것을 고려해 변수간의 관계에서 더 많은 정보를 얻기 위해 변수를 추가한다.
11.6. 로지스틱 회귀(Logistic regression)¶
선형회귀는 다른 유형 종속변수를 처리하도록 일반화 될 수 있다. 만약 종속변수가 부울(boolean)이라면, 일반화 모형은 로지스틱회귀(logistic regression)라 불린다. 만약 종속변수가 정수 갯수라면 포아송회귀(Poisson regression)라 한다.
$$ \log o = \beta_0 + \beta_1 x_1 + \beta_2 x_2 $$
로지스틱 회귀(logistic regression)는 영국의 통계학자인 D. R. Cox가 1958년 에 제안한 확률 모델로서 독립 변수의 선형 결합을 이용하여 사건의 발생 가능성을 예측하는데 사용되는 통계 기법이다.
11.7. 모수 추정¶
모수의 값으로 가장 가능성이 높은 하나의 숫자를 찾아내는 작업을 모수 추정(parameter estimation)이라고 한다.
선형 회귀와 달리 로지스틱 회귀는 닫힌 형식 해답(closed form solution)이 없다. 그래서 초기 해답(solution)을 추측해 반복해 해답을 향해 접근한다.
11.9. 정밀도¶
같은 실험을 반복했을 때 측정결과가 거의 유사한 형태의 집합을 이루면 이는 정밀도가 높다고 말한다. 정확성(accuracy)과 혼동하기 쉬우나, 정확성은 실험값이 이론에 상당하는 참값에 근사한 정도를 나타내며, 정밀도가 높더라도 이 결과의 집합들이 실제의 이론값과는 모두 오차가 클 수 있다.
11.11. 용어사전¶
- 회귀(regression): 데이터에 맞는 모델을 추정하기 위한 몇 가지 관련된 과정중 하나.
- 종속변수(dependent variables): 회귀모형에서 예측 하려는 변수. 또한, 내생 변수로도 알려져 있다.
- 설명변수(explanatory variables): 종속변수를 예측 하거나 설명 하는데 사용되는 변수. 또한 독립 변수 혹은 외생 변수로도 알려져 있다.
- 단순회귀(simple regression): 단지 하나의종속변수, 하나의독립변수만을 갖는회귀.
- 다중회귀(multiple regression): 설명 변수 다수를 갖는 회귀, 하지만 종속 변수는 하나다.
- 선형회귀(linear regression): 선형모형에 기반한 회귀.
- 보통최소제곱(ordinary least squares): 잔차 제곱 오차를 최소화함으로써 모수를 추정하는 선형 회귀.
- 거짓관계(spurious relationship): 모형에 포함 되지 않는 통계적 산출물 혹은 요인으로 발생하여 두 변수가 연관된 두 변수 사이 관계.
- 제어변수(control variable): 거짓관계를 제어 목적으로 혹은 제거 하는데 회귀에 포함되는 변수.
- 대리변수(proxy variable): 다른 요인과 관계 때문에 간접적으로 회귀모형에 정보를 기여하는 변수. 그래서 그 요인에 대한 대리(proxy) 역할을 한다
- 범주형변수(categorical variable): 이산형 순서 없는값을 갖는 변수.
- 결합(join): 두 프레임에 행을 매칭 하는 키(key)를 사용해서 두 데이터 프레임에 데이터를 결합.
- 데이터마이닝(data mining): 많은모형을 검정함으로써 변수 사이 관계를 찾아내는 접근법.
- 로지스틱회귀(logistic regression): 종속변수가부울(boolean) 자료 형식일 때 사용되는 회귀.
- 포아송회귀(Poisson regression): 종속변수가 음수가 아닌 정수로 통상적으로 사용되는 회귀 형식.
- 오즈(odds): 확률(p)를 확률과 확률보(complement) 비율로 $\frac{p}{1-p}$로 나타내는 대안적인 방법.
12. 시계열 분석¶
시계열(time series)은 시스템에서 시간에 따라 변화하는 측정값이다. 간단한 예시로 시간에 따라변하는 세계 평균 기온이 있다.
12.1. 데이터 가져와 정제하기¶
예시로 사용할 데이터는 PriceOfWeed 사이트에서 공개하고 있는 대마초의 시세이다. 다운로드한 데이터를 판다스 데이터프레임으로 만든다. 해당 데이터셋에는 다음과 같은 행이 포함되어 있다.
- city(도시): 도시이름을 나타낸 문자열.
- state(주): 알파벳 두 글자로된 미국의 주.
- price(가격): 지불한 가격(달러).
- amount(수량): 구입한 수량(g).
- quality(품질): 대마초의 품질, 고급, 보통, 저급.
- date(날짜): 대마초 구매 날짜
- ppg: 대마초의 그램당 가격(price per gram)
- lat: 위도(Latitude) 정보
- lon: 경도(Longitude) 정보
import pandas as pd
transactions = pd.read_csv("../data/mj-clean.csv", parse_dates=[5])
transactions.head()
위의 대마초 데이터셋의 데이터를 정제하기 위해 아래 함수를 정의한다.
# 하루단위로 발생한 거래의 평균
def GroupByDay(transactions, func=np.mean):
grouped = transactions[["date", "ppg"]].groupby("date")
daily = grouped.aggregate(func)
daily["date"] = daily.index
start = daily.date[0]
one_year = np.timedelta64(1, "Y")
daily["years"] = (daily.date - start) / one_year
return daily
# 하루단위로 평균 품질
def GroupByQualityAndDay(transactions):
groups = transactions.groupby("quality")
dailies = {}
for name, group in groups:
dailies[name] = GroupByDay(group)
return dailies
dailies = GroupByQualityAndDay(transactions)
12.2. 플롯 그리기¶
시간에 따른 대마초의 가격을 품질로 구분해서 시각화한다.
import matplotlib.pyplot as plt
thinkplot.PrePlot(rows=3)
for i, (name, daily) in enumerate(dailies.items()):
thinkplot.SubPlot(i + 1)
title = "Price per gram ($)" if i == 0 else ""
thinkplot.Config(ylim=[0, 20], title=title)
thinkplot.Scatter(daily.ppg, s=10, label=name)
if i == 2:
plt.xticks(rotation=30)
thinkplot.Config()
else:
thinkplot.Config(xticks=[])
품질에 따른 대마초 가격, 위 그래프를 통해 명백하게 알 수 있는 것은 2013년 11월에 데이터가 빈다는 것이다. 이 기간동안의 결측 데이터를 처리하는 방법은 이후에 생각해볼 것이다.
12.3. 선형회귀(Linear regression)¶
종속 변수 y와 한 개 이상의 독립 변수 (또는 설명 변수) X와의 선형 상관 관계를 모델링하는 회귀분석 기법이다.
import statsmodels.formula.api as smf
def RunLinearModel(daily):
model = smf.ols("ppg ~ years", data=daily)
results = model.fit()
return model, results
for name, daily in dailies.items():
model, results = RunLinearModel(daily)
print(name, results.summary())
def PlotFittedValues(model, results, label=""):
years = model.exog[:, 1]
values = model.endog
thinkplot.Scatter(years, values, s=15, label=label)
thinkplot.Plot(years, results.fittedvalues, label="model", color="#ff7f00")
def PlotLinearModel(daily, name):
model, results = RunLinearModel(daily)
PlotFittedValues(model, results, label=name)
thinkplot.Config(
title="Fitted values",
xlabel="Years",
xlim=[-0.1, 3.8],
ylabel="Price per gram ($)",
)
name = "high"
daily = dailies[name]
PlotLinearModel(daily, name)
높은 품질의 대마초 가격 시계열과 선형 최소자승적합 그래프. 그러나 선형회귀는 시계열 분석에 적합하지 않을 수 있는데 그 이유는 다음과 같다.
- 첫째로 장기추세를 선형 혹은 다른 간단한 함수로 예측할 이유가 없다. 일반적으로 가격은수요와 공급에 의해 결정 된다. 모두 시간에 따라 예측불가한 방향으로 변화 한다.
- 둘째 선형회귀모형은 모든 데이터에 동일한 가중치를 둔다. 예측 목적이라면 최근 데이터에 더 많은 가중치를 두어야 한다.
- 마지막으로 선형회귀의 가정중 하나는 잔차는 상관되지 않는 잡음이라는 것이다. 그러나 시계열 데이터에서 연속된 값은 상관되기 때문에 종종 이런 가정은 틀리다.
12.4. 이동 평균(Moving averages)¶
대부분 시계열 분석은 관측된 계열이 다음 세 가지 요소합이라는 가정에 기반한다.
- 추세(Trend): 지속되는 변경 사항을 잡아 내는 평활 함수(smooth function).
- 계절성(Seasonality): 일별, 주별, 월별 주기적 변동 주기.
- 잡음(Noise): 장기 추세 주위 확률 변동.
앞절에서 살펴봤듯 회귀는 계열에서 추세를 뽑아내는 방법이다. 하지만 추세는 찾는 훌륭한 방법 중 하나는 이동평균(moving average)을 구하는 것이다. 이동평균은 데이터를 윈도우(windows) 단위로 나누고 각 윈도우의 값을 평균내 계산한다.
12.5. 결측값(Missing values)¶
시계열 데이터는 종종 일별, 주별, 월별, 년별 주기(cycle)를 나타낸다. 다음 절에서 계절성을 검정하는 방법을 제시한다. 하지만 결측값에는 잘 동작하지 않는다. 그래서 이 문제를 먼저 해결해야한다. 이런 결측값을 채우는 가장 간단하고 흔한 방법이 이동평균이다.
def PlotRollingMean(daily, name):
dates = pd.date_range(daily.index.min(), daily.index.max())
reindexed = daily.reindex(dates)
thinkplot.Scatter(reindexed.ppg, s=15, alpha=0.2, label=name)
roll_mean = reindexed.ppg.rolling(30).mean()
thinkplot.Plot(roll_mean, label="rolling mean", color="#ff7f00")
plt.xticks(rotation=30)
thinkplot.Config(ylabel="price per gram ($)")
PlotRollingMean(daily, name)
def PlotRollingMean(daily, name):
dates = pd.date_range(daily.index.min(), daily.index.max())
reindexed = daily.reindex(dates)
ewma = reindexed.ppg.ewm(span=30).mean()
resid = (reindexed.ppg - ewma).dropna()
fake_data = ewma + thinkstats2.Resample(resid, len(reindexed))
reindexed.ppg.fillna(fake_data, inplace=True)
thinkplot.Scatter(reindexed.ppg, s=15, alpha=0.2, label=name)
thinkplot.Plot(ewma, label="rolling mean", color="#ff7f00")
plt.xticks(rotation=30)
thinkplot.Config(ylabel="price per gram ($)")
PlotRollingMean(daily, name)
위 그래프에 결측 데이터를 채워 시각화한 결과를 확인할 수 있다. 재표본 추출한 잔차가 무작위(random)로 결과는 매번 달라진다. 나중에 결측값에서 생성된 오차를 어떻게 특성화하는지 살펴볼 것이다.
12.6. 계열상관(Serial correlation)¶
만약 월요일에 대마초 가격이 높다면 다음 몇 일동안 가격이 높을 것을 예상된다. 그리고 만약 낮다면 낮게 유지될 것을 예상 할 수 있다. 이와 같이 각값이 다음값과 상관있는 패턴을 계열 상관(serial correlation) 이라 한다.
12.7. 자기상관(Autocorrelation)¶
시간 또는 공간적으로 연속된 일련의 관측치들간에 존재하는 상관관계다.
12.8. 예측¶
시계열분석은 시간에변화하는 시스템을 조사하고 때로는 예측도 할수있다.
def GenerateSimplePrediction(results, years):
n = len(years)
inter = np.ones(n)
d = dict(Intercept=inter, years=years, years2=years**2)
predict_df = pd.DataFrame(d)
predict = results.predict(predict_df)
return predict
def PlotSimplePrediction(results, years):
predict = GenerateSimplePrediction(results, years)
thinkplot.Scatter(daily.years, daily.ppg, alpha=0.2, label=name)
thinkplot.plot(years, predict, color="#ff7f00")
xlim = years[0] - 0.1, years[-1] + 0.1
thinkplot.Config(
title="Predictions",
xlabel="Years",
xlim=xlim,
ylabel="Price per gram ($)",
loc="upper right",
)
name = "high"
daily = dailies[name]
_, results = RunLinearModel(daily)
years = np.linspace(0, 5, 101)
PlotSimplePrediction(results, years)
def SimulateResults(daily, iters=101, func=RunLinearModel):
_, results = func(daily)
fake = daily.copy()
result_seq = []
for _ in range(iters):
fake.ppg = results.fittedvalues + thinkstats2.Resample(results.resid)
_, fake_results = func(fake)
result_seq.append(fake_results)
return result_seq
def GeneratePredictions(result_seq, years, add_resid=False):
n = len(years)
d = dict(Intercept=np.ones(n), years=years, years2=years**2)
predict_df = pd.DataFrame(d)
predict_seq = []
for fake_results in result_seq:
predict = fake_results.predict(predict_df)
if add_resid:
predict += thinkstats2.Resample(fake_results.resid, n)
predict_seq.append(predict)
return predict_seq
def PlotPredictions(daily, years, iters=101, percent=90, func=RunLinearModel):
result_seq = SimulateResults(daily, iters=iters, func=func)
p = (100 - percent) / 2
percents = p, 100 - p
predict_seq = GeneratePredictions(result_seq, years, add_resid=True)
low, high = thinkstats2.PercentileRows(predict_seq, percents)
thinkplot.FillBetween(years, low, high, alpha=0.3, color="gray")
predict_seq = GeneratePredictions(result_seq, years, add_resid=False)
low, high = thinkstats2.PercentileRows(predict_seq, percents)
thinkplot.FillBetween(years, low, high, alpha=0.5, color="gray")
years = np.linspace(0, 5, 101)
thinkplot.Scatter(daily.years, daily.ppg, alpha=0.1, label=name)
PlotPredictions(daily, years)
xlim = years[0] - 0.1, years[-1] + 0.1
thinkplot.Config(
title="Predictions", xlabel="Years", xlim=xlim, ylabel="Price per gram ($)"
)
위 그래프는 선형적합에 기초한 예측, 표집오차, 예측오차에 기인한 변동성을 시각화한다. 짙은 회색 지역은 표집 오차에 대한 90% 신뢰구간을 나타낸다. 즉, 표집 때문에 추정 기울기와 절편에 대한 불확실성이다. 밝은 회색 지역은 예측오차에 대한 90% 신뢰구간을 보여 주는데 표집오차와 확률 변동의 합이다.
12.11. 용어사전¶
- 시계열(time series): 시간(timestamp)과 연관된 데이터셋.
- 윈도우(window): 이동평균을 계산하기 위해 사용되는 시계열에 연속값 서열.
- 이동평균(moving average): 겹쳐지지 않는 일련의 윈도우에 대한 평균을 계산함으로써 시계열에 잠재하는 추세를 추정하는데 사용 되는 여러 통계량 중의 하나.
- 이동평균(rolling mean): 각 윈도우 평균값에 기반한 이동평균.
- 지수가중이동평균(exponentially-weighted moving average, EWMA): 가중 평균에 기반한 이동평균으로 가장 최근값에 가장 높은 가중치를 두고 이전 값에 대해서는 지수적으로 줄어드는 가중치를 둔다.
- 스팬(span): 가중치가 얼마나 빨리 줄어드는지를 결정하는 EWMA 모수.
- 계열상관(serial correlation): 한 시계열과 이동된 혹은 시차 이동한 시계열과 상관.
- 시차(lag): 계열 상관 혹은 자기 상관에서 이동 크기.
- 자기상관(autocorrelation): 임의의 시간차이를 갖는 계열 상관에 대한 좀 더 일반적인 용어.
- 자기상관함수(autocorrelation function): 시차에서 계열상관으로 매핑하는 함수.
- 정상성(stationary): 만약 모수와 잔차분포가 시간에 따라 변화하지 않는다면 모델이 정상성이 있다.
13. 생존분석¶
생존분석(Survival analysis)은 무언가 얼마나 지속하는지를 기술하는 방법이다. 종종 사람생명연구에 사용되지만, 또한 기계나 전자부품의 생존(survial) 혹은 좀 더 일반적으로 사건 전 시간 간격에도 적용된다.
생존 기간을 분석하여 생존함수(survival function)를 추정하는 통계기법으로 치료방법이나 예후 인자 등을 추정하는 데 이용된다. 실험 기간 동안 약물처리, 유전염기서열 변이 와 같이 특정 조건을 만족하는 집단 내 개체가 각 측정 시간에서 살아남아 있는 빈도로 나타낸다.
13.1. 생존곡선(Survival curves)¶
생존 곡선(survival curve; $S(t)$) 는 존속 t를 t보다 더 오래 생존할 확률로 치환하는 함수다. 만약 존속(duration) 분포 즉, 수명(lifetimes)을 알고 있다면 생존곡선을 찾는 것은좀 더 쉽다. CDF의 여분포(complement)가 된다. 여기서는 저자가 만든 survival 모듈을 사용한다.
NSFG 데이터셋에서 결과 코드(outcome code) 1, 3, 4은 각각 정상출산, 사산, 유산을 나타낸다. 여기에서는 분석을 위해서 유발유산(induced abortion), 자궁외임신(ectopic pregnancy), 그리고 임신상태인 경우는 제외한다.
import survival
complete = df.query("outcome in [1, 3, 4]").prglngth
cdf = thinkstats2.Cdf(complete, label="cdf")
def MakeSurvivalFromCdf(cdf, label=""):
ts = cdf.xs
ss = 1 - cdf.ps
return survival.SurvivalFunction(ts, ss, label)
sf = MakeSurvivalFromCdf(cdf, label="survival")
thinkplot.Plot(sf)
thinkplot.Cdf(cdf, alpha=0.2)
plt.legend()
위 그래프에서 임신기간 CDF와 생존함수의 상보적 CDF가 보여진다.
13.2. 위험함수¶
생존 함수에서 위험 함수(hazard function)를 도출할 수 있다. 임신기간에 대해서 위험 함수는 시간 t에서 t까지 계속 되고 나서 t에서 끝나는 임신 비율이다. 좀 더 정확하게는 다음과같다.
$$\lambda(t) = \frac{S(t)-S(t+1)}{S(t)}$$
분자는 PMF(t)로 t에서 끝나는 수명 비율이다. SurvivalFunction는 MakeHazard 함수를 제공 하는데 위험 함수를 계산한다
13.4. 케플란-마이어 추정(Kaplan-Meier estimation)¶
카플란-마이어 생존분석은 생존 분석에서 사용되는 통계 기법이다. 비모수 통계를 이용하여 생존함수를 추정한다. 미국의 통계학자 폴 마이어와 에드워드 카플란에 의해 개발되었다.
생존 분석에서 중심적인 알고리즘 중의 하나인 캐플란-마이어 추정(Kaplan-Meier estimation)을 알아보기 위해 이번 예제는 결혼을 하지 않은 미혼 여성의 관측 정보를 포함해 사용한다.
13.5. 결혼 곡선(Marriage curve)¶
사용한 데이터셋에는 다음의 정보가 포함되어 있다.
- cmbirth: 응답자에 대한 생년월일.
- cmintvw: 응답자를 인터뷰한 날짜.
- cmmarrhx: 응답자가 결혼한 날짜.
- evrmarry: 만약 응답자가 인터뷰한 날짜 이전에 결혼했다면 1, 그렇지않으면 0.
resp6 = pd.read_csv("../data/2002FemResp.csv")
resp6.cmmarrhx.replace([9997, 9998, 9999], np.nan, inplace=True)
resp6["agemarry"] = (resp6.cmmarrhx - resp6.cmbirth) / 12.0
resp6["age"] = (resp6.cmintvw - resp6.cmbirth) / 12.0
complete = resp6[resp6.evrmarry == 1].agemarry.dropna()
ongoing = resp6[resp6.evrmarry == 0].age
from collections import Counter
def EstimateHazardFunction(complete, ongoing, label="", verbose=False):
if np.sum(np.isnan(complete)):
raise ValueError("complete contains NaNs")
if np.sum(np.isnan(ongoing)):
raise ValueError("ongoing contains NaNs")
hist_complete = Counter(complete)
hist_ongoing = Counter(ongoing)
ts = list(hist_complete | hist_ongoing)
ts.sort()
at_risk = len(complete) + len(ongoing)
lams = pd.Series(index=ts, dtype="float64")
for t in ts:
ended = hist_complete[t]
censored = hist_ongoing[t]
lams[t] = ended / at_risk
if verbose:
print(t, at_risk, ended, censored, lams[t])
at_risk -= ended + censored
return survival.HazardFunction(lams, label=label)
hf = EstimateHazardFunction(complete, ongoing) # 위험함수
sf = hf.MakeSurvival() # 생존함수 추정
초혼 연령에 대한 위험 함수 시각화
thinkplot.Plot(hf)
thinkplot.Config(xlabel="Age (years)")
초혼 연령에 대한 생존 함수
thinkplot.Plot(sf)
thinkplot.Config(xlabel="Age (years)", ylabel="Prob unmarried")
13.7. 신뢰구간¶
캐플란-마이어분석(Kaplan-Meier analysis)은생존곡선에대한단하나의추정 값을산출한다. 하지만, 추정값의불확실성을정량화하는것도또한중요하다
초혼연령에대한생존함수와가중재표집에근거한90% 신뢰구간.
def EstimateMarriageSurvival(resp):
complete = resp[resp.evrmarry == 1].agemarry.dropna()
ongoing = resp[resp.evrmarry == 0].age
hf = EstimateHazardFunction(complete, ongoing)
sf = hf.MakeSurvival()
return hf, sf
def ResampleSurvival(resp, iters=101):
_, sf = EstimateMarriageSurvival(resp)
thinkplot.Plot(sf)
low, high = resp.agemarry.min(), resp.agemarry.max()
ts = np.arange(low, high, 1 / 12.0)
ss_seq = []
for _ in range(iters):
sample = thinkstats2.ResampleRowsWeighted(resp)
_, sf = EstimateMarriageSurvival(sample)
ss_seq.append(sf.Probs(ts))
low, high = thinkstats2.PercentileRows(ss_seq, [5, 95])
thinkplot.FillBetween(ts, low, high, color="gray", label="90% CI")
ResampleSurvival(resp6)
thinkplot.Config(
xlabel="Age (years)",
ylabel="Prob unmarried",
xlim=[12, 46],
ylim=[0, 1],
loc="upper right",
)
위 그림은 앞에서 추정한 생존함수와 함께 표집 가중치에 근거한 90% 신뢰구간을 회색선으로 나타낸 것이다. 표집 가중치에 의해 두 곡선사이에 차이가 생긴다는 것을 유념하라.
13.8. 코호트 효과(Cohort effects)¶
생존 분석의 문제 중 하나는 추정 곡선의 다른 부분이 응답자의 다른 집단에 기반한다는것이다. 시점 t에 곡선 부분은 인터뷰 당시에는 적어도 응답자 연령이 적어도 t인 응답자에 기반한다. 그래서 곡선의 가장 왼쪽 부분은 모든 응답자로 부터 데이터가 포함되어있고, 가장 오른쪽에는 가장 나이든 응답자만 포함된다.
데이터의 특성이 시간에 따라 변화 하지 않는다면 문제가 되지않을 것이다. 그러나 출생한 세대에 따라 여성의 혼인 패턴은 다를 것으로 이기 때문에 10년 단위로 응답자 집단을 만들어 효과를 조사해보자.
이렇게 비슷한 사건으로 그룹화한 집단을 코호트(cohorts)라 부르며 집단간 차이를 코호트 효과(cohort effects)라 부른다.
resp5 = pd.read_csv("../data/FemResp_1995.csv")
resp6 = pd.read_csv("../data/FemResp_2002.csv")
resp7 = pd.read_csv("../data/FemResp_2006.csv")
resps = [resp5, resp6, resp7]
def AddLabelsByDecade(groups, **options):
thinkplot.PrePlot(len(groups))
for name, _ in groups:
label = "%d0s" % name
thinkplot.Plot([15], [1], label=label, **options)
def EstimateMarriageSurvivalByDecade(groups, **options):
thinkplot.PrePlot(len(groups))
for _, group in groups:
_, sf = EstimateMarriageSurvival(group)
thinkplot.Plot(sf, **options)
def PlotResampledByDecade(resps, iters=11, predict_flag=False, omit=None):
for i in range(iters):
samples = [thinkstats2.ResampleRowsWeighted(resp) for resp in resps]
sample = pd.concat(samples, ignore_index=True)
groups = sample.groupby("decade")
if omit:
groups = [(name, group) for name, group in groups if name not in omit]
if i == 0:
AddLabelsByDecade(groups, alpha=0.7)
if predict_flag:
PlotPredictionsByDecade(groups, alpha=0.1)
EstimateMarriageSurvivalByDecade(groups, alpha=0.1)
else:
EstimateMarriageSurvivalByDecade(groups, alpha=0.2)
PlotResampledByDecade(resps)
thinkplot.Config(xlabel="Age (years)", ylabel="Prob unmarried", xlim=[13, 45])
plt.legend()
위의 그림은 코호트간의 생존 함수 시각화를 나타낸 것으로 몇가지 패턴이 보인다.
- 50년대 태어난 여성의 초혼 연령이 가장 낮다.
- 이후 세대들의 초혼 연령은 점차 늦어지고 있다.
- 80년대생의 경우 25살이 넘어서 결혼을 하지 않는 비율이 유지되는 것이 보인다.
13.9 외삽법(Extrapolation)¶
앞에서 배운 예시에서 80년대 90년대생의 데이터가 부족해 충분한 곡선의 시각화를 할 수 없었다. 데이터를 빌려(borrowing)와 곡선을 외삽(extrapolate)해보도록 하자.
def PlotPredictionsByDecade(groups, **options):
hfs = []
for _, group in groups:
hf, sf = EstimateMarriageSurvival(group)
hfs.append(hf)
thinkplot.PrePlot(len(hfs))
for i, hf in enumerate(hfs):
if i > 0:
hf.Extend(hfs[i - 1])
sf = hf.MakeSurvival()
thinkplot.Plot(sf, **options)
PlotResampledByDecade(resps, predict_flag=True)
thinkplot.Config(xlabel="Age (years)", ylabel="Prob unmarried", xlim=[13, 45])
plt.legend()
50년대생의 데이터를 빌려와 이후 세대의 초혼 연령을 추정했다. 비혼 여성의 비율이 시대가 흐를수록 늘어남이 보인다.
13.12. 용어사전¶
- 생존분석(survival analysis): 일반적으로 사건이 일어나기까지 시간을 기술하고 예측하는 방법.
- 생존곡선(survival curve): 시점t의 t의 생존 확률로 치환하는 함수.
- 위험함수(hazard function): 시점t의 t시점까지 생존한 사람이 t시점에서 사망한 비율을 치환하는 함수.
- 캐플란-마이어 추정(Kaplan-Meier estimation): 카플란-마이어 생존분석은 관찰 시간에 따라서 사건이 발생한 시점에서의 사건 발생률을 계산하는 방법이다.
- 코호트(cohort): 비슷한 속성으로 그룹화된 집단.
- 코호트 효과(cohort effect): 코호트 사이의 차이.
14. 해석적 방법(Analytic methods)¶
해석적 방법이란 논리적인 추론을 통해 통계 문제를 풀어가는 과정이다. 이번 장에서 해석적 방법 일부를 설명하고 어떻게 동작하는지 알아본다.
14.1. 정규분포¶
인간과 자연 세상에서 일어나는 수많은 일을 설명하는 핵심 개념으로 정규분포는 평균값을 기준으로 좌우 대칭 형태가 나타나며, 좌우 극단으로 갈수록 급격하게 수치가 낮아지는 특징이 있다.
14.2. 표집분포¶
통계의 의사결정을 위한 이론적 분포로 이론적으로 표본의 크기가 n인 표본을 무한히 반복 추출한 후 무한개의 표본들의 평균을 가지고 그린 분포이다.
14.4. 중심극한정리¶
만약 값의 분포가 평균(µ), 표준편차(σ)를 갖는다면 합 분포는 근사적으로 $N(nµ, nσ^2)$이 된다는 것이 중심극한 정리의 의미다.
중심극한정리(Central Limit Theorem, 이하 CLT)는 통계 분석에 가장 유용한 도구중 하나로 몇가지 주의할 점은 다음과 같다.
- 값은 독립적으로 추출 되어야 한다. 만약 값들이 서로 상관 된다면 CLT을 적용 할 수 없다.
- 값들은 동일한 분포에서 나와야 하며 유한한 평균과 분산을 갖는 분포를 가져야한다. 그래서 대부분 파레토 분포는 해당되지 않는다.
- 수렴의 속도는 분포 왜도에 의존 한다. 지수분포에서 나온 값들의 합은 작은 n에 대해서 수렴한다. 로그정규분포에서 나온값들의 합은 더 커다란 크기가 필요 하다.
14.5. 중심극한정리 검정¶
지수 분포값의 합에 대한 분포를 검정한다.
def MakeExpoSamples(beta=2.0, iters=1000):
samples = []
for n in [1, 10, 100]:
sample = [np.sum(np.random.exponential(beta, n)) for _ in range(iters)]
samples.append((n, sample))
return samples
def NormalPlotSamples(samples, plot=1, ylabel=""):
for n, sample in samples:
thinkplot.SubPlot(plot)
thinkstats2.NormalProbabilityPlot(sample)
thinkplot.Config(
title="n=%d" % n,
legend=False,
xticks=[],
yticks=[],
xlabel="random normal variate",
ylabel=ylabel,
)
plot += 1
thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeExpoSamples()
NormalPlotSamples(samples, plot=1, ylabel="sum of expo values")
반면에 파레토분포의 경우는 다음과 같다.
def MakeParetoSamples(alpha=1.0, iters=1000):
samples = []
for n in [1, 10, 100]:
sample = [np.sum(np.random.pareto(alpha, n)) for _ in range(iters)]
samples.append((n, sample))
return samples
thinkplot.PrePlot(num=3, rows=2, cols=3)
samples = MakeParetoSamples()
NormalPlotSamples(samples, ylabel="sum of Pareto values")
파레토 분포는 로그 정규분포보다 더 기울어짐이 심하다. 모수에 따라 많은 파레토 분포는 유한 평균과 분산을 갖지 못한다. 따라서 중심극한정리가 적용 되지 않는다. 심지어 n=100
일때 정규 확률 그림은 직선과 거리가 매우 멀다.
14.7. 상관검정¶
앞서 신생아 체중과 산모 연령사이 상관을 계산하는데 순열 검정(permutationtest)을 사용 해 p-value가 0.001 보다 작아 통계적으로 유의하다는 결론은 내렸다. 이제 해석적 방법을 통해 같은 분석을 수행해보자.
import scipy
def StudentCdf(n):
ts = np.linspace(-3, 3, 101)
ps = scipy.stats.t.cdf(ts, df=n - 2)
rs = ts / np.sqrt(n - 2 + ts**2)
return thinkstats2.Cdf(rs, ps)
def ResampleCorrelations(live):
live2 = live.dropna(subset=["agepreg", "totalwgt_lb"])
data = live2.agepreg.values, live2.totalwgt_lb.values
ht = CorrelationPermute(data)
ht.PValue()
return len(live2), ht.actual, ht.test_cdf
n, r, cdf = ResampleCorrelations(live_baby)
model = StudentCdf(n)
thinkplot.Plot(model.xs, model.ps, color="gray", alpha=0.5, label="Student t")
thinkplot.Cdf(cdf, label="sample")
thinkplot.Config(xlabel="correlation", ylabel="CDF", legend=True, loc="lower right")
위 그림은 재표본추출로 생성한 분포를 따른 분포를 보여준다. 두 분포가 거의 동일(iden tical)하는 것을 알 수 있다. 따라서 실제 분포가 정규 분포는 아니지만, 피어슨상관계수는표본평균과 분산에 기반하며 중심극한정리에 의해 값을 계산할 수 있다.
t = r * np.sqrt((n - 2) / (1 - r**2))
p_value = 1 - scipy.stats.t.cdf(t, df=n - 2)
print(r, p_value)
p-value 결과는 매우작다. 따라서 해석적 방법을 통해서도 신생아의 체중과 산모의 나이에는 통계적 유의성이 있다고 결론내릴 수 있다.
14.8. 카이제곱 검정¶
카이제곱 통계량은 기대값과 정규화된 총 편차는 측정한다. 카이제곱 통계량이 널리 사용되는 이유는 귀무가설 아래에 표집분포가 해석적이기 때문이다. scipy 라이브러리를 사용해 카이제곱 검정을 쉽게 사용할 수 있다.
from scipy.stats import chi2_contingency
data = [[30, 10], [35, 5], [28, 12], [20, 20]]
chi2, pval, dof, expected = chi2_contingency(data)
print(f"카이제곱 검정의 p-value는 {pval} 이다")
14.9. 해석적 방법의 장점¶
- 해석적 방법은 설명과 이해하기 더 쉽다. 귀무 가설을 시뮬레이션하고 검정 통계량을 계산하는것이 좀 더 명확하기 때문이다.
- 좀 더 강건한 모델을 생성한다. 해석적 방법은 더 적은 가정을 요구하고 좀 더 쉽게 모델을 확장할 수 있다.
- 해석적 방법은 디버그(debuggable)할 수 있다. 검증 및 테스트(incremental development and testing)를 할 수 있어 결과에 신뢰를 준다.
그러나 해석적 방법에 장점만 있는 것은 아니다. 해석적 방법은 계산이 오래 걸린다는 단점이 있다.
15. 마치며¶
Thinkstat2은 무료로 공개된 통계학 책으로 영어로 작성되어 있다. 책의 번역은 오래전 이광춘 씨가 인터넷에 공개하였다. 여기에서는 해당 번역본의 어색한 표현을 수정하고 내용을 요약해 정리한 것이다.