패키지 불러오기

library(knitr) # R의 결과를 표로 재구성하는 데 유용한 패키지
library(kableExtra) # knitr의 표 스타일을 보조하는 패키지
library(ggplot2)
library(patchwork)
library(tidyverse)

확률 관련 포스팅에서는 몇 가지 간단한 확률을 R을 통해서 증명하고(Some simple probability demonstrations), 정규분포(normal distributions)와 이항분포(binomial distributions)로부터 분위(quantiles)를 얻어내는 방법을 살펴보고자 합니다.

여기서는 일단 간단하게 정규분포랑 이항분포, 그리고 독립사건과 종속사건만 간단하게 살펴보겠습니다.

• 분위(Quantiles)에 대해 어떻게 설명할까 고민을 좀 했는데 좋은 블로그 포스팅을 찾아서 참조용으로 링크하겠습니다.
• 링크는 여기를 보시면 됩니다. Q-Q Plot 설명하면서 Quantile의 정의도 잘 정리해놓으셨습니다.

확률이란?

확률모델은 우리가 어떠한 이론-모델을 가지고 있을 때, 어떠한 데이터를 관측할 확률에 대해 알려줍니다. 동전 던지기를 예로 들어보겠습니다. 동전 던지기를 할 때, 우리는 압면을 관측할 확률에 대해서 대강 알고 있습니다. 누구든지 물어보면 0.5의 확률이라고 말할 것입니다. 동전을 던졌을 때, 얻을 수 있는 가능성이 앞면과 뒷면이라는 두 가지밖에 없고 이 동전에 아무런 문제가 없을 때, 그 기대값이 0.5일 것이라는 일종의 이론, 모델을 가지고 있기 때문입니다. 정확히는 동전을 반복해서 던졌을 때(=관측이 반복될 때) 데이터가 어떻게 나타날 것(앞면 뒷면이 반반)이라는 자료생성과정에 관한 이론을 우리가 가지고 있다는 것입니다. 우리가 어떤 모델, 이론을 가지고 있을 때, \(Y\)라는 결과를 얻을 확률을 아래와 같이 적어보겠습니다.

\[ Pr(Y|M) = Pr(\text{데이터}|\text{모델}) \]

다음으로는 확률의 종류에 대해서 알아보도록 하겠습니다. 예전에 배우셨던 집합 개념을 조금 떠올리실 필요가 있습니다. 먼저, 사건 A와 사건 B가 있다고 하겠습니다. 그리고 이 두 사건이 동시에 일어날 확률을 A와 B에 대한 결합확률(joint probability)이라고 하겠습니다. 집합 기호로 나타내면 \(Pr(A\cap B)\)로 나타낼 수 있으며, \(Pr(A\;\text{and}\;B)\)라고도 합니다.

한계확률은 전체 중에서 특정 사건이 발생할 확률을 말합니다. 예를 들어서, 아래와 같은 상황이 있다고 하겠습니다. A라는 국가와 B라는 국가가 서로 갈등하는 국면에 놓여있는 상황입니다. A와 B 모두가 싸우기로 결정할 확률은 \(Pr(A_{\text{Fight}}\;\cap\;B_{\text{Fight}}) = 0.3\)이라고 할 수 있습니다. 반대로 두 국가 모두가 서로 싸움을 회피하고 상대방의 요구에 순응할 확률은 \(Pr(A_{\text{Comply}}\;\cap\;B_{\text{Comply}}) = 0.4\)입니다. 이 두 경우 모두 결합확률이죠. 그렇다면 A가 싸울 한계확률은 얼마일까요? 전체 중에 \(A_{\text{Fight}}\)일 경우이므로 0.4라고 할 수 있습니다. 마찬가지로 \(Pr(A_{\text{Comply}}) = 0.6\)이며, 같은 논리로 \(Pr(B_{\text{Fight}}, B_{\text{Comply}})\)는 각각 0.5, 0.5 라고 할 수 있습니다.

그럼 위의 표를 바탕으로 조건부확률(conditional distribution)도 한 번 살펴보겠습니다. 조건부확률이란 한계확률 대비 결합확률의 비율이라고 핤 수 있습니다. A와 B라는 사건으로 일반화된 진술로 표현하자면, \(Pr(A|B) = P(A \cap B)/P(B)\)라고 할 수 있으며, 이때 조건부확률은 사건 B 일어났을 경우에 A라는 사건이 나타날 확률을 의미합니다.

하나 더 살펴보도록 하겠습니다. 유권자들에게 A라는 후보에게 투표했는지, 혹은 B라는 후보에게 투표했는지 묻고, 동시에 그들이 남성인지 여성인지에 대해 물었다고 해보겠습니다.

이때, 한계확률과 조건부확률, 그리고 결합확률을 각각 살펴보겠습니다.

• 먼저 응답자가 여성일 한계확률은 전체 200명 중에 여성 100명이므로 \(Pr(Female) = 100/200 = 0.5\)라고 할 수 있습니다.

• 그렇다면 이 여성 중에서 A 후보에게 투표했을 조건부확률은 전체 여성 100명 중 A 후보에게 투표한 65명, 즉 \(Pr(\text{A}|Female) = 65/100 = 0.65\)입니다.

• 그렇다면 A후보에 투표했으면서 동시에 여성일 확률은 어떻게 될까요? 일단 이 두 사건이 독립적이라는 가정 하에서 우리는 이 두 사건의 확률을 곱해주면 됩니다. 즉, \(Pr(\text{A}, Female) = Pr(\text{A})\;\;\text{AND}\;\;Pr(Female) = 65/200 = 0.325\)라는 얘기입니다.

조건부확률은 한계확률에 대한 조건부 확률의 비율이라고 이해할 수 있습니다:

• \(Pr(A|B) = Pr(A, \;B)/Pr(B)\)

• \(Pr(\text{A}|Female) = Pr(\text{A}, \;\;Female)/Pr(Female)\)

• \(Pr(\text{A}|Female)/Pr(Female) = 0.325/0.5 = 0.65\)

결합확률, 한계확률, 그리고 조건부확률은 앞으로를 개념들을 이해하는 데 중요하니까 머리에 잘 담아두시고, 이제 다음으로 넘어가겠습니다.

주사위굴리기 게임!

확률을 공부할 때, 지겹도록 등장하는 놈들이 총 셋이 있습니다. 동전, 카드, 그리고 주사위입니다. 인류는 아마도 이 셋을 만듦으로써 스스로를 괴롭히는 통계학을 발전시켰는지도 모르겠습니다… 일단 주사위 굴리기는 직관적으로 확률과 통계를 이해하기 좋은 방식입니다. 먼저 주사위 하나를 한 번 굴리는 것을 시뮬레이팅하는 함수를 코딩해보겠습니다.

die <- as.integer(runif(1, min=1, max=7))
die

## [1] 3

die는 ?runif 라고 입력하여 살펴보면 generates random deviates.라고 기술되어 있는 것을 확인할 수 있습니다. 이어지는 함수를 풀어서 설명하면 다음과 같습니다: > 다음의 결과를 정수의 형태로 저장하라(as.integer) \(\rightarrow\) 무작위로 다음의 범주 내에서 다른 값을 1번 추출하라 \(\rightarrow\) 최소값은 1, 최대값은 6 (1 이상 7미만)을 갖게하라.

그러면 이번에는 두 개의 주사위를 굴려보도록 하겠습니다. 굴리는 횟수는 한 번입니다.

dice <- (as.integer(runif(1, min=1, max=7))) +
  (as.integer(runif(1, min=1, max=7)))
dice

## [1] 7

여기서 + 연산자는 부울리안 논리에 따르면 OR를 의미합니다. 즉 각 주사위를 한 번씩 랜덤으로 굴려 얻는 값을 더한 결과를 dice에 저장하라는 명령입니다. 그럼 이번에는 100번, 1000번, 그리고 100만번을 돌려겠습니다.

# 주사위 두 개를 100번 던져보기
dice100 <-  (as.integer(runif(100, min=1, max=7))) +
  (as.integer(runif(100, min=1, max=7)))

dice100 %>% table() %>% kable(col.names = c("Dice", "Freq"))

# 주사위 두 개를 1,000번 던져보기
dice1000 <-  (as.integer(runif(1000, min=1, max=7))) +
  (as.integer(runif(1000, min=1, max=7)))

dice1000 %>% table() %>% kable(col.names = c("Dice", "Freq"))

# 주사위 두 개를 100만 번 던져보기
dice10000 <-  (as.integer(runif(1e4, min=1, max=7))) + 
  (as.integer(runif(1e4, min=1, max=7)))

dice10000 %>% table() %>% kable(col.names = c("Dice", "Freq"))

이렇게 시뮬레이팅한 세 결과를 히스토그램으로 살펴보겠습니다.

Simulation of Dice Rollings

일단 가시적으로 살펴볼 수 있게 각 코드 이후에 그 결과값의 빈도를 표로 나타내보았습니다. 그리고 그 표를 히스토그램으로 재구성해보았습니다. 역시 N이 늘어날 수록 우리(?)가 사랑하는 그 녀석의 모습이 드러나기 시작합니다.

주사위는 1에서 6까지의 한정된 값을 가지고, 두 개를 합쳐서 굴려봐야 2부터 12까지의 한정된(bounded) 값이긴 하지만 이 주사위 굴리기를 통해서 우리는 지난 번 포스팅에서 살펴보았던 것처럼 정규분포(normal distribution)와 표본 크기(n), 혹은 표집(sampling)의 관계를 간접적으로 다시 한 번 살펴볼 수 있습니다.

동전 던지기

그럼 이번에는 동전을 한 번 던져보겠습니다. 저는 아직까지 앞면 뒷면 이외에 옆면에도 표기를 지닌 동전을 본 적이 없으니, 여기서의 동전도 앞면과 뒷면이라는 두 개의 값만을 가진다고 가정하겠습니다. 삼면이나 사면을 가진 동전을 보신 분들은 부디 댓글로 알려주시길… 나와라 검은 백조야(김웅진 · 김지희 2012, p.53).

아무튼 앞면과 뒷면이 있는 경우에 그 각각이 나올 확률은 0.5, 0.5라고 할 수 있습니다. rbinom()은, randomly [drawn] binomial, 무작위로 이항변수를 추출하라는 함수라고 할 수 있습니다. 백문이 불여일코드.

coin <- rbinom(5, 1, .5)
coin

## [1] 1 1 0 0 0

이어지는 함수를 풀어서 설명하면 다음과 같습니다.: 이항변수로 무작위로 추출하라(rbinom()) \(\rightarrow\) 5번 추출하라 \(\rightarrow\) 최대값은 1 (=최소값은 0) \(\rightarrow\) 추출확률은 0.5. 즉, 1이 나올 확률을 50%로 설정하여 무작위로 추출하라는 것입니다. 그럼 이번에는 100개의 동전을 던져보겠습니다.

# 동전 100개 던지기
coin100 <- rbinom(100, 1, .5)
coin100 %>% table() %>% kable(col.names = c("Coin", "Freq"))

#동전 1000개 던지기
coin1000 <- rbinom(1000, 1, .5)
coin1000 %>%  table() %>% kable(col.names = c("Coin", "Freq"))

#동전 1000개 던지기
coin10000 <- rbinom(10000, 1, .5)
coin10000 %>%  table() %>% kable(col.names = c("Coin", "Freq"))

Simulation of Coin Flippings

이와 같이 이항변수는 나뉘어진(discrete) 값을 가집니다. 히스토그램으로 그리면 0이 나오는 빈도랑 1이 나오는 빈도만 보여주는 것이지요. 이번에는 100개의 동전을 1000번 던지는 경우를 시뮬레이팅해보겠습니다.

coin1Mx <- rep(NA, 1e6)
for (i in 1:1e6) {
  coin1Mx[i] <- sum(rbinom(100, 1, .5))
}

coin1Mx %>% as_tibble() %>% 
  ggplot(aes(value)) + 
  geom_histogram(aes(y = ..density..), 
                 color = "black", fill = "white") + 
  labs(x = "Number of heads") + 
  theme_bw()

100만 개의 셀을 결측치(NA)로 갖는 텅빈 coinMx라는 벡터를 만들어보겠습니다. 그리고 i가 1에서 100만까지 반복되는 loop를 구성합니다. coin1Mx_1부터 coin1Mx_1000000까지 총 100만개의 coin1Mx_n들은 모두 100개의 동전을 던져서 앞면(1)이 나오는 경우의 수를 더한 각각의 값을 가질 것입니다. 따라서 coin1Mx는 100만개의 요소를 지닌 벡터입니다.

이 자료를 활용해서 만약 100개를 던졌을 때, 앞면이 60번 나오는 것이 과연 극단적인 확률일지 아니면 무던한 것일지 확인해보겠습니다. 먼저, 60번 이상 앞면이 나온 경우를 세는 함수를 짜보겠습니다.

coin1Mx[coin1Mx > 60] %>% table() %>% 
  kable(col.names = c("Head > 60", "Freq."))

앞면이 나오는 빈도가 61번 이상부터는 점차 감소하는 것을 확인할 수 있습니다. 그렇다면 이번에는 앞면이 60번 넘게 나올 확률을 구해보겠습니다.

# 앞면이 60번 넘게 나온 시뮬레이션 횟수/전체 시뮬레이션 횟수
length(coin1Mx[coin1Mx > 60])/length(coin1Mx)

## [1] 0.01756

# 앞면이 60번 넘게 나온 시뮬레이션 횟수를 표로 나타내서 총합을 구한 뒤
# 시뮬레이션 전체 횟수로 나누어 줍니다: 위와 동일합니다.
sum(table(coin1Mx[coin1Mx > 60]))/1e6

## [1] 0.01756

length()는 count()랑 같은 개념이라고 할 수 있습니다. 전체 coin1Mx의 수, 즉 100만 건 중에서 앞면이 60번보다 많이 나온 경우의 수가 어느 정도인지를 계산한 것입니다. 즉, 앞면이 60번보다 더 많이 나올 확률은 매우 작다고 할 수 있습니다다. 앞면이 60번 초과하여 나올 확률은 “평균적” or “일반적”인 것은 아니라는 뜻입니다.

독립사건 시뮬레이션

두 개의 독립적인 사건을 시뮬레이션해보겠습니다. 첫 번째 함수는 rand1이라는 자료에 최소값 0, 최대값 10 미만의 값을 가지는 100개의 수를 무작위로 담으라는 명령입니다. 두 번째 함수는 정규분포를 따라 평균이 0이고 표준편차가 2의 값을 갖는 분포에서 100개의 값을 무작위로 추출해 rand2라는 자료에 담으라는 명령입니다.

• 일단 보면 rand2의 경우, 비록 명령으로 평균을 0으로 하라고 설정했지만 무작위 추출 결과 100개의 값의 평균은 0보다 약간 큰 것을 확인할 수 있습니다.

• 이제 두 값을 각각 \(x\)축, \(y\)축으로 설정하여 100개의 값을 좌표에 매칭시킨 그래프로 나타내보겠습니다.

rand1 <- runif(100, min = 0, max = 10)
summary(rand1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2882  2.5464  5.9284  5.2430  7.2565  9.9499

rand2 <- rnorm(100, mean = 0, sd = 2)
summary(rand2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -4.7610 -1.3517  0.1949  0.1156  1.4799  4.5755

rand.dat <- tibble(rand1 = rand1,
                   rand2 = rand2)
rand.dat %>% ggplot(aes(x = rand1, y = rand2)) + 
  geom_point(alpha = 0.3, size = 2, fill = "#752021", shape = 21) + 
  theme_bw()

결과적으로 이 그래프에서 우리는 rand1과 rand2 간에는 어떠한 경향성을 발견하기 힘듭니다. 즉, 독립적이라는 것은 두 변수 간에 어떠한 관계를 상정할 수 없다는 것으로 이해할 수 있습니다.

종속사건 시뮬레이션

이번에는 두 사건이 종속적 관계에 있는, 즉 한 사건에서의 변화가 다른 사건에 영향을 미치는 경우를 시뮬레이션해 보겠습니다.

rand3는 아까 만들어놓은 rand1에다가 0.75를 곱해서 4를 더한 100개의 값에다가 사실 상 rand2와 같은 방식으로 구한 값을 더하여 구한 자료입니다. 즉, 이 값은 rand1의 값에 어떠한 조치를 취하여 얻은 값이므로 rand1에 영향을 받은 결과물이라고 할 수 있습니다. 이렇게 구한 rand3와 rand1 간의 관계를 그래프로 나타내봅시다.

rand.dat <- rand.dat %>%
  mutate(
    rand3 = 4 + 0.75 * rand1 + rnorm(100, mean = 0, sd = 2)
  )

rand.dat %>% ggplot(aes(x = rand1, y = rand3)) + 
  geom_point(alpha = 0.3, size = 2, fill = "#752021",  shape = 21) + 
  theme_bw()

이 그래프의 해석은 나중에 산포도(scatter plot) 및 단순회귀분석(simple regression)을 살펴볼 때 다시 한 번 다룰 것입니다. 일단 여기서는 rand1이 증가할 때, rand3도 증가세를 보이는—둘의 인과적 관계는 확인할 수 없지만 아무튼—양(positive)의 관계를 보이고 있다는 것을 알 수 있습니다.

cor(rand1, rand3)로 두 변수 간의 상관계수를 구해보면¹, 양수의 기울기를 얻게 됩니다.

이번에는 보다 더 불확실성을 가지는 종속사건의 관계를 시뮬레이션해보겠습니다. 여기서 말하는 불확실성은 더 큰 표준편차를 의미합니다.

• 표준편차가 평균에서 각 관측치가 떨어져 있는 거리의 평균이라고 할 때,
• 표준편차 값이 크다는 것은 개별 값들이 평균에서 더 넓게 분포해 있다는 것을 의미합니다.

rand.dat <- rand.dat %>% 
  mutate(
    rand4 = 4 + 0.75 * rand1 + rnorm(100, mean = 0, sd = 5)
  )

rand.dat %>% ggplot(aes(x = rand1, y = rand4)) + 
  geom_point(alpha = 0.3, size = 2, fill = "#752021", shape = 21) + 
  theme_bw()

cor(rand.dat$rand1, rand.dat$rand4)

## [1] 0.385872

## 데이터 프레임을 가지고 상관관계를 구하는 방법: corrr 패키지
## install.packages("corrr")
library(corrr)
rand.dat %>% select(rand1, rand4) %>% correlate(quiet = T) %>% 
  kable(col.names = c("Variables", "rand1", "rand4"), digits = 3)

rand1과 rand4의 산포도는 좀 더 넓게 퍼진 모양새를 보이고 둘의 상관관계는 식에서 상정한 0.75라는 기울기보다 더 낮아지는 것을 확인할 수 있습니다.