8. 확률 분포

8. 확률 분포 (Probability distributions)

R이 실제 데이터 분석에서 하루도 안 빠지고 쓰이는 가장 큰 이유 중 하나가, 통계 분포를 다루는 내장 함수가 풍부하다는 점이에요. 확률을 계산하거나, 분포 그래프를 그리거나, 분포에서 난수를 뽑는 작업까지 전부 함수 한 줄로 해결됩니다. 이번 장에서는 R이 제공하는 분포 관련 함수들이 어떻게 구성되어 있는지, 그리고 실제 데이터의 분포를 살펴보고 두 집단을 비교하는 검정까지 한 흐름으로 익혀볼게요.

출처: R 공식 매뉴얼

본문

8.1 R은 통계표의 집합이기도 해요

R을 편리하게 쓰는 방법 하나는, R을 포괄적인 통계표 세트로 활용하는 거예요. R에는 다음을 구해주는 함수들이 마련되어 있죠.

누적 분포 함수 P(X ≤ x) — 확률변수가 x 이하일 확률
확률 밀도 함수 (probability density function)
분위수 함수 (quantile function) — q가 주어졌을 때 P(X ≤ x) > q 를 만족하는 가장 작은 x
그리고 그 분포에서 난수(모의) 생성

어떤 분포가 지원되는지 다음 표에서 확인할 수 있어요. 첫 열이 분포 이름, 두 번째 열이 R에서 쓰는 이름(접두어를 붙이는 부분), 세 번째 열이 그 분포에 필요한 추가 인자에요.

Distribution	R name	additional arguments
beta	beta	shape1, shape2, ncp
binomial	binom	size, prob
Cauchy	cauchy	location, scale
chi-squared	chisq	df, ncp
exponential	exp	rate
F	f	df1, df2, ncp
gamma	gamma	shape, scale
geometric	geom	prob
hypergeometric	hyper	m, n, k
log-normal	lnorm	meanlog, sdlog
logistic	logis	location, scale
negative binomial	nbinom	size, prob
normal	norm	mean, sd
Poisson	pois	lambda
signed rank	signrank	n
Student’s t	t	df, ncp
uniform	unif	min, max
Weibull	weibull	shape, scale
Wilcoxon	wilcox	m, n

여기 표의 이름 앞에 문자 하나를 붙여서 역할을 정해요. d 는 밀도(density), p 는 누적 분포 함수(CDF), q 는 분위수 함수(quantile function), r 은 난수 생성(simulation, random deviates)을 뜻해요. 예를 들어 정규분포라면 dnorm, pnorm, qnorm, rnorm처럼 쓰는 식이죠.

첫 번째 인자도 함수마다 정해져 있어요. dxxx 는 x, pxxx 는 q, qxxx 는 p, rxxx 는 n을 받아요. 다만 rhyper, rsignrank, rwilcox만 예외적으로 표본 개수로 nn을 받아요. 또 비중심성 모수(non-centrality parameter)인 ncp는 모든 분포에서 아직 지원되지는 않으니, 자세한 내용은 온라인 도움말을 확인해 보세요.

pxxx 와 qxxx 계열 함수에는 논리값 인자 lower.tail 과 log.p가 있고, dxxx 계열에는 log가 있어요. 이걸 이용하면 예를 들어 누적(또는 "적분된") 위험 함수(hazard function) H(t) = -log(1 - F(t))를 다음처럼 바로 구할 수 있어요.

 - pxxx(t, ..., lower.tail = FALSE, log.p = TRUE)

아니면 dxxx(..., log = TRUE) 로 더 정확한 로그우도(log-likelihood)를 직접 구할 수도 있구요.

덧붙여, 정규분포 표본의 studentized range 분포를 다루는 ptukey와 qtukey, 다항분포(multinomial)를 다루는 dmultinom과 rmultinom 함수도 있어요. 그 외 분포는 기여 패키지(contributed packages)에서 찾을 수 있는데, 특히 SuppDists 패키지가 유명해요.

몇 가지 예시를 볼게요.

> ## 2-tailed p-value for t distribution
> 2*pt(-2.43, df = 13)
> ## upper 1% point for an F(2, 7) distribution
> qf(0.01, 2, 7, lower.tail = FALSE)

난수 생성이 R에서 어떻게 이뤄지는지(RNG 가 어떻게 동작하는지)는 온라인 도움말을 참고하세요.

8.2 데이터의 분포 살펴보기

(일변량) 데이터가 주어졌을 때, 그 분포를 살펴보는 방법은 정말 많아요. 가장 간단한 건 숫자 그 자체를 들여다보는 거예요. summary 와 fivenum 은 조금 다른 두 종류의 요약값을 주고, stem 은 숫자들을 "줄기-잎(stem and leaf)" 그림으로 보여줘요.

> attach(faithful)
> summary(eruptions)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  1.600   2.163   4.000   3.488   4.454   5.100
> fivenum(eruptions)
[1] 1.6000 2.1585 4.0000 4.4585 5.1000
> stem(eruptions)

  The decimal point is 1 digit(s) to the left of the |

  16 | 070355555588
  18 | 000022233333335577777777888822335777888
  20 | 00002223378800035778
  22 | 0002335578023578
  24 | 00228
  26 | 23
  28 | 080
  30 | 7
  32 | 2337
  34 | 250077
  36 | 0000823577
  38 | 2333335582225577
  40 | 0000003357788888002233555577778
  42 | 03335555778800233333555577778
  44 | 02222335557780000000023333357778888
  46 | 0000233357700000023578
  48 | 00000022335800333
  50 | 0370

줄기-잎 그림은 히스토그램과 비슷한데, R에는 히스토그램을 그리는 hist 함수가 따로 있어요.

> hist(eruptions)
## make the bins smaller, make a plot of density
> hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)
> lines(density(eruptions, bw=0.1))
> rug(eruptions) # show the actual data points

더 매끄러운 밀도 그림은 density 로 그릴 수 있고, 위 예시에서 그 density 로 만든 선을 추가했어요. 대역폭 bw 는 시행착오로 고른 값인데, 기본값은 너무 많이 평활화(smoothing)되기 때문이에요(보통 "흥미로운" 밀도에서는 기본값이 과평활화되기 마련이죠). 더 나은 자동 선택 방법도 있고, 이 예시에서는 bw = "SJ" 가 좋은 결과를 줘요.

경험적 누적 분포 함수(empirical CDF)는 ecdf 를 써서 그릴 수 있어요.

> plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE)

이 분포는 어떤 표준 분포와도 거리가 멀어 보이네요. 그렇다면 오른쪽 모드(mode), 즉 3분보다 긴 분출만 따로 보면 어떨까요? 정규분포를 적합해서 그 적합된 CDF를 겹쳐 그려볼게요.

> long <- eruptions[eruptions > 3]
> plot(ecdf(long), do.points=FALSE, verticals=TRUE)
> x <- seq(3, 5.4, 0.01)
> lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3)

이걸 더 꼼꼼히 살펴보려면 분위수-분위수(Q-Q) 그림이 도움이 돼요.

par(pty="s")       # arrange for a square figure region
qqnorm(long); qqline(long)

이 그림은 나름 괜찮은 적합을 보여주지만, 오른쪽 꼬리가 정규분포에서 기대하는 것보다 짧아요. 이번엔 t 분포에서 뽑은 모의 데이터와 비교해 볼게요.

x <- rt(250, df = 5)
qqnorm(x); qqline(x)

이건 (난수 표본이므로 보통) 정규분포에서 기대하는 것보다 더 긴 꼬리를 보여줄 거예요. 생성 분포를 상대로 Q-Q 그림을 그리려면 이렇게 하면 돼요.

qqplot(qt(ppoints(250), df = 5), x, xlab = "Q-Q plot for t dsn")
qqline(x)

마지막으로 정규성 일치 여부를 더 공식적으로 검정하고 싶을 수 있어요. R에는 Shapiro-Wilk 검정이 있어요.

> shapiro.test(long)

         Shapiro-Wilk normality test

data:  long
W = 0.9793, p-value = 0.01052

그리고 Kolmogorov-Smirnov 검정도 있구요.

> ks.test(long, "pnorm", mean = mean(long), sd = sqrt(var(long)))

         One-sample Kolmogorov-Smirnov test

data:  long
D = 0.0661, p-value = 0.4284
alternative hypothesis: two.sided

(여기서 정규분포의 모수를 같은 표본에서 추정했기 때문에 분포 이론이 엄밀하게는 성립하지 않는다는 점에 주의하세요.)

8.3 일표본 검정과 이표본 검정

지금까지는 하나의 표본을 정규분포와 비교했어요. 훨씬 더 흔한 작업은 두 표본의 측면을 비교하는 거예요. 참고로 R에서 여기 쓰일 검정을 포함한 "고전적(classical)" 검정들은 전부 기본적으로 로드되는 stats 패키지 안에 있어요.

얼음의 융해 잠열(latent heat) 데이터를 생각해 볼게요. Rice (1995, p.490)에 나오는, 얼음이 녹을 때의 열량(cal/gm) 데이터예요.

Method A: 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97
          80.05 80.03 80.02 80.00 80.02
Method B: 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97

상자그림(boxplot)은 두 표본을 간단하게 그래프로 비교해 주는 방법이에요.

A <- scan()
79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97
80.05 80.03 80.02 80.00 80.02

B <- scan()
80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97

boxplot(A, B)

이 그림은 첫 번째 집단 결과가 두 번째보다 대체로 높은 경향이 있음을 보여줘요.

두 집단의 평균이 같은지 검정하려면, 짝이 없는(unpaired) t-검정을 쓸 수 있어요.

> t.test(A, B)

         Welch Two Sample t-test

data:  A and B
t = 3.2499, df = 12.027, p-value = 0.00694
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.01385526 0.07018320
sample estimates:
mean of x mean of y
 80.02077  79.97875

정규성을 가정했을 때 유의미한 차이가 있음을 보여주는 결과예요. 기본적으로 R 함수는 두 표본의 분산이 같다고 가정하지 않아요.

두 표본이 정규 모집단에서 나왔다면, F 검정으로 분산이 같은지도 확인할 수 있어요.

> var.test(A, B)

         F test to compare two variances

data:  A and B
F = 0.5837, num df = 12, denom df =  7, p-value = 0.3938
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.1251097 2.1052687
sample estimates:
ratio of variances
         0.5837405

유의미한 차이가 없다는 결과므로, 분산이 같다고 가정하는 고전 t-검정을 쓸 수 있어요.

> t.test(A, B, var.equal=TRUE)

         Two Sample t-test

data:  A and B
t = 3.4722, df = 19, p-value = 0.002551
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.01669058 0.06734788
sample estimates:
mean of x mean of y
 80.02077  79.97875

이 검정들은 모두 두 표본이 정규성을 따른다고 가정해요. 반면 이표본 Wilcoxon(또는 Mann-Whitney) 검정은 귀무가설 하에서 두 표본이 공통의 연속 분포를 따른다는 가정 하나만 필요해요.

> wilcox.test(A, B)

         Wilcoxon rank sum test with continuity correction

data:  A and B
W = 89, p-value = 0.007497
alternative hypothesis: true location shift is not equal to 0

Warning message:
Cannot compute exact p-value with ties in: wilcox.test(A, B)

경고 메시지에 주목하세요. 각 표본에 동률(ties)이 여러 개 있는데, 이는 이 데이터가 (아마 반올림 때문에) 이산 분포에서 왔음을 강하게 시사해요.

두 표본을 그래픽으로 비교하는 방법도 여럿 있어요. 위에서 상자그림 한 쌍을 이미 봤죠. 아래 코드는 두 표본의 경험적 CDF 두 개를 보여줘요.

> plot(ecdf(A), do.points=FALSE, verticals=TRUE, xlim=range(A, B))
> plot(ecdf(B), do.points=FALSE, verticals=TRUE, add=TRUE)

qqplot 은 두 표본의 Q-Q 그림을 그려주고, Kolmogorov-Smirnov 검정은 공통 연속 분포를 가정했을 때 두 ecdf 사이의 최대 수직 거리를 검정해요.

> ks.test(A, B)

         Two-sample Kolmogorov-Smirnov test

data:  A and B
D = 0.5962, p-value = 0.05919
alternative hypothesis: two-sided

Warning message:
cannot compute correct p-values with ties in: ks.test(A, B)

더 알아보기

각 분포 함수의 보다 정확한 인자와 동작은 R 온라인 도움말(?pnorm, ?dnorm 등)을 참고하세요. 특히 ncp(비중심성 모수) 지원 여부는 분포마다 달라요.
난수 생성의 원리와 시드(seed) 관리에 관해서는 ?RNG 도움말을 확인해 보세요.
검정의 전제 조건(정규성, 분산 동일성 가정 등)을 이해하면 t.test, var.test, wilcox.test 같은 검정 결과를 올바르게 해석할 수 있어요.