Distribution into R

패키지 불러오기

library(knitr) # R의 결과를 표로 재구성하는 데 유용한 패키지
library(kableExtra) # knitr의 표 스타일을 보조하는 패키지
library(ggplot2)
library(patchwork)
library(tidyverse)
library(futurevisions) # 그래픽 컬러 패키지
theme_set(theme_bw()) # 자동으로 ggplot 테마를 bw로 설정

분포(Distribution)

이전 주차의 확률 파트에서 설명한 것과 같이 확률이란 “하나의 사건이 일어날 수 있는 가능성을 수치로 나타낸 것”이라고 할 수 있으며, 일반적으로 사건 A가 나타날 확률을 P(A)라고 표현합니다. 통계학에서 확률이란 일종의 경험적 확률로 한 사건 A가 일어날 확률을 P(A)라고 할 때, \(n\)번의 반복 시행에서 사건 A가 나타난 횟수를 \(r\)이라고 하면, 상대도수 \(\frac{r}{n}\)\(n\)이 커질수록 확률 P(A)에 가까워진다는 것을 알 수 있습니다. 즉, \(\lim r/n = P(A)\)인 경우, 이 P(A)를 사건 A의 통계적 확률이라고 합니다. 이에 기초해서 조건부 확률과 한계확률에 대한 내용을 떠올려주시면 될 것 같습니다.

그렇다면 분포란 무엇일까요? 우리가 관심을 가지고 있는 것은 바로 확률분포(probability distribution)입니다. 확률분포란 "한 번의 시행에서 변수 \(X\)가 취할 수 있는 값에 대하여 \(X\)가 취할 수 있는 값, \(x_1, x_2, x_3, \dots, x_n\)에 대응하는 확률을 \(p_1, p_2, p_3, \dots, p_n\)이라고 할 때, 이 관계를 확률분포라고 합니다. 그럼 확률변수란 무엇이냐? 무작위로 어떠한 값을 가질 수 있는 변수 정도로 이해하실 수 있을 것 같습니다.

확률분포의 종류

통계분석에 앞서서 일반적으로 알아두어야 할 분포 함수(distribution functions)를 살펴보겠습니다. 첫 번째 그림은 확률분포의 종류를 간략하게 보여주는 것입니다. 일종의 자료생성과정(data generating process)가 어떠한 확률분포를 따라 이루어졌을 것이라고 기대하느냐에 따라 우리는 우리가 얻은 표본의 관측치가 모집단 수준에서는 어떠한 형태의 분포를 가지게 될 것이라고 기대할 수 있고, 그 기대에 따라 다른 통계방법을 사용할 수 있습니다. 위의 분포들 중에는 이 스터디에서 다루지 않는 범위에 있는 분포들이 존재하고, 당장은 굳이 아실 필요가 없다고 생각하여 제외한 분포들도 존재합니다(예, 웨이블 분포, 베타 분포 등). 이전 데이터의 종류에서 살펴보았던 연속형과 이산형의 차이에 유의하셔서 이해하시면 좋을 것 같습니다.

기초적인 분포의 개념을 직관적으로 이해하기 위하여 정규분포(normal distributions)와 이항분포(binomial distributions)로부터 분위(quantiles)를 얻어내는 방법을 살펴보겠습니다.

정규분포 (The normal distribution)

일단 정규분포라는게 어떻게 생겨먹은 놈인지 한 번 살펴보도록 하겠습니다.

Normal distribution (Mean = 0, SD = 1)

Normal distribution (Mean = 0, SD = 1)

위의 정규분포를 살펴보면, 일종의 기준점(thresholds)을 볼 수 있습니다. 즉, 전체 관측치의 약 2%가 그 기준점보다 값이 작거나 클것 등이 그런 것이죠. 나중에 가설검정과 유의수준에 대한 논의를 할 때 자세하게 살펴보도록 하겠습니다. 수학의 정석 확률과 통계 파트에서 배웠던 1.96 등과 같은 숫자들을 보게 될 것입니다. 위의 정규분포에서 좌측의 끝에 분포한 경우를 좌측꼬리(lower tail), 오른쪽 끝의 분포를 우측꼬리(upper tail)이라고 부릅니다.

pnorm(70, mean = 50, sd = 10, lower.tail = TRUE)
## [1] 0.9772499
pnorm(70, mean = 50, sd = 10, lower.tail = FALSE)
## [1] 0.02275013
1 - pnorm(70, mean = 50, sd = 10, lower.tail = TRUE)
## [1] 0.02275013

첫번째 코드는 50을 평균으로 하고 10을 표준편차로 하는 정규분포가 있을 때, 그 분포에서 70이라는 숫자는 어디에 위치하는지를 묻는 것입니다.

두 번째도 동일한 의미인데, 두 식의 차이는 lower.tail 옵션을 어떻게 설정하느냐에 달려있습니다.

  • 각 식이 도출한 결과를 보면 이해하시겠지만 근본적으로 두 식은 동일하며, 좌측에서 누적확률을 계산할 것인지 우측으로부터 계산할 것인지의 차이일 뿐입니다.

  • 간단하게 말하면 첫 번째 식은 70이라는 숫자는 이 분포에서 하위 97.7%에 위치한 값이다라고 말하는 것이고 두 번째 식은 상위 2.2%라고 말하는 것입니다.

  • 따라서 두 번째 식은 전체 확률에서 첫 번째 식으로 계산한 확률을 제한 값과 같으므로 세 번째의 형태로 계산할 수도 있습니다.

그렇다면 만약에 양측꼬리 확률(two-tailed probability)에서 최소한 70만큼 ‘극단적인’(extreme) 확률을 얻고 싶을 때는 어떻게 해야할까요? 이 경우는 생각을 좀 달리 해보면 됩니다. 평균을 기점으로 70은 오른쪽 끝쪽에 위치하는 셈입니다. 그만큼 왼쪽 끝에 위치한 값을 상정하고 그 값이 나올 확률을 함께 계산해주어야 합니다.

이 경우 평균 50에서 70은 20만큼 떨어져 있습니다(우측으로, + 방향). 따라서 우리는 좌측으로 20만큼 떨어진 30이 나올 확률을 함께 고려해주어야 하는 것입니다. 이러한 결과를 얻는 데는두 가지 방법이 존재합니다.

## 첫 번째 방법
pnorm(70, mean = 50, sd = 10, lower.tail = FALSE) + 
  pnorm(30, mean = 50, sd = 10, lower.tail = TRUE)
## [1] 0.04550026
## 두 번째 방법
2 * pnorm(70, mean = 50, sd = 10, lower.tail = FALSE)
## [1] 0.04550026

그렇다면 이렇게 구한 분포에서의 누적확률 값을 가지고 분위를 구하여 보겠습니다. 마찬가지로 평균 50에 표준편차가 10인 분포를 상정합니다. 이때 사용할 함수는 qnorm()입니다.

qnorm(0.9772499, mean = 50, sd = 10, lower.tail = TRUE)
## [1] 70.00001

역으로 계산한 것인데, 평균 50에 표준편차 10인 분포에서 좌측부터 앞서 구한 누적확률에 해당하는 값을 구하라는 명령입니다. 앞에서 우리가 입력한 70과 근사한 값을 얻을 수 있습니다. 근소한 차이는 소수점에 의해 발생하는 것으로 이해하시면 됩니다.

다음으로는 주어진 평균 50, 표준편차 10의 분포에서 70이라는 값이 분포에서 차지하는 밀도를 확인해보겠습니다 밀도(density)를 알아보기 위한 함수는 다음과 같습니다.

dnorm(70, mean = 50, sd = 10)
## [1] 0.005399097

그렇다면 이번에는 정규분포에서 관측치를 무작위로 추출(draws)을 해보겠습니다. 이번에도 함수에는 70이라는 값이 사용될 것인데, 여기서 사용되는 70은 특정한 값이 아니라 추출 횟수를 의미합니다.

x <- rnorm(70, mean = 50, sd = 10)
x
##  [1] 72.57308 40.05966 46.37258 44.17327 60.14351 74.55990 63.24215 54.05914
##  [9] 36.88879 63.43313 51.27676 42.93698 48.44602 40.89836 38.45322 63.27353
## [17] 61.26935 55.32106 65.52355 49.31888 57.68101 25.70053 50.03684 35.57660
## [25] 36.04556 39.71587 36.39946 56.48622 35.47841 72.15947 45.02195 55.63756
## [33] 60.69671 35.59950 63.58326 44.27626 55.46370 50.83937 46.59655 44.22666
## [41] 49.97470 38.23799 59.68103 49.14727 53.80333 37.89642 57.98342 78.17659
## [49] 44.52032 69.66339 36.60864 44.03805 53.88238 36.33232 53.21062 45.22242
## [57] 38.41020 54.88670 41.45718 59.78519 62.74681 49.02517 44.10639 48.78397
## [65] 59.74982 65.24785 55.63426 48.71448 51.23788 44.44237

총 70개의 값이 무작위로 추출되어 x라는 벡터에 담겨 있는 것을 확인할 수 있습니다.

이번에는 좀 더 구체적으로 정규분포 사례들을 살펴보겠습니다. 서로 다른 평균(mean)과 표준편차(standard deviation)를 가지는 세 개의 정규분포를 그려보겠습니다.

normal5 <- rnorm(n = 10000, mean = 5, sd = 3)
normal50 <- rnorm(n = 10000, mean = 50, sd = 10)
normal20 <- rnorm(n = 10000, mean = 20, sd = 1)

세 함수 모두 표본의 크기는 1만 개이며, 각자 다른 평균과 표준편차를 따르는 정규분포를 가정하여 무작위로 추출된 값을 담는 벡터로 출력됩니다. 이렇게 값을 갖는 세 개의 벡터를 하나의 데이터로 합쳐보겠습니다. 간단하게 말하면 하나의 표로 합친다는 것과 같습니다.

norm <- bind_rows(tibble(x = normal5, Mean = 5, SD = 3),
                  tibble(x = normal50, Mean = 50, SD = 10),
                  tibble(x = normal20, Mean = 20, SD = 1))
norm <- norm %>% mutate(
  Mean = as.factor(Mean))

그리고 이렇게 만들어진 데이터프레임의 평균 값을 Factor 자료 유형으로 변환해줍니다. 이것이 의미하는 게 뭘까요? 평균5, 평균50, 평균20을 문자열로 인식하게 하여 일종의 “이름”으로 인식하게 만드는 것입니다. 그러면 이제 이렇게 만들어진 데이터로 그래프를 그려보겠습니다.

ggplot(norm, aes(x = x, fill = factor(Mean))) +
  geom_density(adjust = 4, alpha = 1/2) +
  guides(color=guide_legend(title = "Mean, SD")) +
  guides(fill=guide_legend(title = "Mean, SD")) +
  scale_fill_manual(values = futurevisions("mars")) + 
  scale_color_manual(values = futurevisions("mars")) + 
  labs(x = "") + 
  theme(legend.position = "top")
Probability Density Function of Normal Distribution

Probability Density Function of Normal Distribution

ggplot2 패키지는 R에서 그래프를 그리는 데 있어서 유용하게 사용됩니다. 먼저 ggplot2 패키지를 불러오고 나서, ggplot(데이터프레임 이름, aes(x = x축으로 삼을 변수이름))을 설정하면 일단 수학적으로 지정한 데이터의 변수에 대한 기본적인 작업은 진행이 됩니다.

그러나 이 단계에서는 아직 가시화(visualization)라고 할 수는 없는데, 주어진 데이터를 컴퓨팅했을 뿐이지 어떻게 가시적인 형태로 보여주라는 명령을 부여하지 않았기 때문입니다. 컴퓨터로 그림을 그려본 사람들은 이해가 쉬운데, ggplot()은 레이어 시스템을 이용해서, 우리가 뭔가 그려진 걸 얻고 싶을 때마다 레이어를 하나씩 추가해서 보여달라고 R에게 요구해야만 합니다.

그리고 이때, 레이어를 추가하는 것은 +로 가능합니다. 하나씩 뜯어서 보면,

  • ggplot(norm, aes(x = x)) + : ggplot2 패키지를 이용하여 norm이라는 데이터 프레임에서 x축에는 x라는 변수를 기준으로 늘어놓아라. 그리고 뒤에 더해지는 레이어 명령을 덮어 씌워라 라는 의미의 코드입니다.

  • geom_density(aes(fill = as.factor(Mean)), adjust = 4, alpha = 1/2) +: 밀도함수의 형태로 그래프를 그리되, Mean이라는 Factor 변수가 같은 것들 끼리 같은 색으로 칠해서 보여주어라 라는 코드입니다.

    • 뒤의 옵션은 세세한 조정이니 굳이 언급하지는 겠습니다.
    • 뒤에 더해지는 옵션들은 레이어를 추가하여 기존의 코드에 덮어 씌워, 그래프를 덧칠하는 역할을 합니다.
  • guides(color=guide_legend(title = "Mean, SD")) + guides(fill=guide_legend(title = "Mean, SD")) +: 우측에 더해지는 레전드의 이름을 어떻게 지으라는 명령입니다.

    • 이 경우에는 색도 있고 그래프를 선명하게 보여주는 선도 있기 때문에 그 각각이 레이어를 이루고 있어 둘 모두에게 이러한 명령어를 적용해야 합니다.
  • scale_color_discrete(labels = c("5, 3", "20, 1", "50, 10")) + scale_fill_discrete(labels = c("5, 3", "20, 1", "50, 10")) +: colorfill을 뒤에 이어지를 라벨로 구분해주라는 코드입니다.

  • ggtitle(Probability Density Function\nNormal Distribution) : 표 전체의 이름을 지정하는 명령어이고, \nR에서는 강제 개행(enter)하는 명령어입니다.

이항분포 (Binomial distribution)

이항분포는 우리에게 n번의 시행에서 k번 성공할 확률을 보여주는 분포입니다. 단순하게 말하자면 0과 1의 값만 갖는다고 가정된 벡터가 있다고 하겠습니다. 이때 벡터의 요소의 총 개수는 100개이고 랜덤으로 0과 1 중 하나가 벡터에 들어간다고 할 수 있습니다. 이때 전체 요소의 수 100개 중 1이 뽑혔을 경우를 계산하면 결국 1을 뽑을 성공 사례의 총합을 100으로 나눔으로써, 전체 대비 성공의 확률을 구하는 것이라 할 수 있습니다. 즉, 정규분포에서와는 다르게 이항분포에서는 평균(mean)이 아니라 비율(proportion)에 초점을 맞추게 됩니다.

pbinom(27, size=100, prob=0.25, lower.tail = TRUE)
## [1] 0.7223805

시도가 성공할 확률을 0.25라고 놓고 시행 횟수를 100번으로 가정한 경우입니다. 즉, 100번 시도했을 때 성공할 확률이 25%인 이항분포를 가지는 데, 실제로는 27번 성공했다고 한다면 누적확률에서 어느 위치에 해당하는지를 묻는 함수입니다. lower.tail()TRUE로 설정되어 있으므로 좌측에서부터 계산한 것입니다. 즉, 27번 성공한 것은 하위 72.2%, 상위 27.8%에 해당한다는 것을 알 수 있습니다.

qbinom(0.7223805, size = 100, prob = 0.25, lower.tail = TRUE)
## [1] 27

이렇게 구한 누적확률로 다시금 성공 횟수를 추정해보겠습니다. 정규분포때와 거의 유사합니다. size = 100은 서로 독립적인 사건의 시행, 즉 베르누이 시행(Bernoulli trials)의 횟수를 의미합니다. 그렇다면 100번 시도했을 때 27번 성공할 정확한 확률을 구해보겠습니다.

dbinom(27, size=100, prob=0.25)
## [1] 0.08064075
choose(100, 27)*.25^27*(1-.25)^(100-27)
## [1] 0.08064075

두 가지 방법으로 구할 수 있는데, 첫 번째는 R에 내장된 기본 함수인 dbinom()으로 구하는 것입니다. 두 번째는 choose() 함수를 이용해서 구하는 것입니다. 개인적으로는 dbinom() 함수가 있는데 굳이 choose() 함수 사용법까지 알아야 하나 싶기는 합니다. 그러나 choose() 함수로 보여지는 식이 좀 더 직관적으로 이해하기에는 도움이 됩니다.

rbinom(27, size=100, prob=0.25)
##  [1] 31 21 17 25 30 28 23 27 25 24 30 31 24 27 27 27 22 17 29 21 27 27 10 22 26
## [26] 24 28
sd(rbinom(27, size=100, prob=0.25))
## [1] 3.772747

마지막으로는 0.25의 확률을 가진 100번의 베르누이 시행을 27번을 반복하는 결과를 보여줍니다. 즉, 평균적으로는 100번 시행 중 25번의 성공을 할 것으로 기대되지만, 실제 시행을 27번 해보면 무작위 추출이기 때문에 25를 중심으로 표준편차 4.21의 범위 내에서 여러 값들이 추출되는 것을 확인할 수 있습니다. 무작위로 추출하기 때문에 돌릴 때마다 값은 다르게 나올 것입니다.

다시 한 번, n번의 독립적인 베르누이 사건을 시행할 때, k번 성공할 확률을 구하는 이항분포로 실습해보겠습니다. 구체적으로 모집단의 성공 확률을 p, 모집단의 크기를 n으로 특정하겠습니다.

binom10 <- rbinom(n = 10000, p = .5, size = 10)
binom50 <- rbinom(n = 10000, p = .5, size = 50)
binom100 <- rbinom(n = 10000, p = .5, size = 100)

정규분포와 다르지는 않습니다. 다만 여기서의 size는 전체 동전을 던지는 횟수로 이해하면 되고 10,000은 그렇게 동전을 10번, 50번, 100번 던지는 걸 10,000번 반복한다는 의미라고 할 수 있습니다. 그럼 이제 마찬가지로 하나의 데이터프레임으로 세 번의 시도(binom10, binom50, binom100)의 결과를 합쳐주고 그래프를 그려보겠습니다.

binom <- bind_rows(tibble(k = binom10, Size = 10), 
                   tibble(k = binom50, Size = 50), 
                   tibble(k = binom100, Size = 100)) %>%
  mutate(Size = as.factor(Size))

ggplot(binom, aes(x = k)) +
  geom_density(aes(group = Size, color = Size, fill = Size), 
               adjust = 4, alpha = 1/2) +
  scale_color_manual(values = futurevisions("mars")) + 
  scale_fill_manual(values = futurevisions("mars")) + 
  labs(x = "") + 
  theme(legend.position = "top")
Probability Mass Function of Binomial Distribution

Probability Mass Function of Binomial Distribution

포와송 분포 (Poison distribution)

이 포스팅은 주로 R 코드에 관한 것이기 때문에 분포에 대한 수리통계적인 설명은 가급적 피하도록 하겠습니다. 포와송 분포는 고정된 대규모 모집단(fixed large population)에서 짧은 시간에 걸쳐서 희소한 사건(rare events)의 발생 횟수를 추정하는 데 용이한 분포입니다.

포와송 분포에서 그 희소한 사건의 발생 확률은 시간 단위 별 발생의 평균 횟수로 나타나며 그리스어 람다(Lambda, \(\lambda\))로 표기됩니다. 람다는 평균과 분산을 결정합니다.

포와송 분포의 확률을 이용해서 우리는 단일 시간 단위에서 k라는 희소한 사건을 정확하게 관측할 확률을 기술할 수 있습니다. 나머지는 앞의 정규분포랑 이항분포에서 했던 코드를 기계적으로 반복해서 살펴보는 것과 같습니자. 단, 포와송 분포에서 코드 중에는 이전과 다른 부분이 있으니 그 점만 유의하면 될 것 같습니다.

pois1 <- rpois(n = 10000, lambda = 1)
pois2 <- rpois(n = 10000, lambda = 2)
pois5 <- rpois(n = 10000, lambda = 5)
pois20 <- rpois(n = 10000, lambda = 20)
pois <- tibble(Lambda.1 = pois1, 
               Lambda.2 = pois2,
               Lambda.5 = pois5,
               Lambda.20 = pois20)
pois # 처음 만들어진 pois 데이터

이번에는 만들어진 pois라는 데이터를 tidyrgather() 함수를 이용해서 다른 변수로 재구성해줄 것입니다. 또한 정해진 값에서 변수명으로 변환할 때 문자열을 추출하기 위한 str_extract() 함수는 stringr 패키지에서 로드하도록 하겠습니다.

pois <- pois %>% tidyr::gather("Lambda", "x")
pois # gather 함수를 이용해 wide -> long으로 바뀐 데이터
pois <- pois %>% mutate(
  Lambda = Lambda %>% 
    str_remove("Lambda.") %>%
    parse_factor(., 
                 levels = c("1", "2", "5", "20"), 
                 ordered = T, include_na = F)
)
pois # Lambda라는 문자열을 지우고 숫자만 factor 로 변환

처음 만든 pois 데이터에서 tidyrgather()를 이용해서 Lambda라는 변수 아래에 1, 2, 5, 20의 라벨을 길게(long-shape) 합쳤습니다. 그리고 마지막으로는 Lambda 변수의 문자열 Lambda.을 제외한 라벨 숫자만 남겨놓은 것입니다. 결과적으로 Lambda 변수에는 각 Lambda 값이 얼마였는지를 나타내는 순위를 가진 요인형 변수값만 남았습니다. 이렇게 만들어진 포와송 분포의 값들을 하나의 그래프로 만들어서 시각화해보도록 하겠습니다.

ggplot(pois, aes(x = x)) +
  geom_density(aes(group = Lambda, 
                   color = Lambda, 
                   fill = Lambda), 
               adjust = 4, alpha = 1/2) +
  scale_color_manual(values = futurevisions("mars")) + 
  scale_fill_manual(values = futurevisions("mars")) + 
  labs(x = "") + 
  theme(legend.position = "top")
Probability Mass Function of Poisson Distribution

Probability Mass Function of Poisson Distribution

음이항 분포 (Negative Binomial Distribution)

음이항 분포는 n번째 시도에서 k번째에 성공할 확률을 보여줍니다. 음이항 분포는 다음과 같은 네 조건이 충족될 때 유용합니다.

  1. 모든 시도는 독립적이다.
  2. 각 시도의 결과는 성공 혹은 실패로 구분될 수 있다.
  3. 성공 확률 (p)은 각 시도마다 동일하다.
  4. 마지막 시도는 항상 성공이어야만 한다.

여기서 보면 알 수 있듯이, 처음 세 조건은 이항분포와 동일합니다. 세 개의 음이항 분포 결과를 계산하여 세 개의 벡터로 저장하겠습니다. 그리고 이걸 마찬가지로 하나의 데이터로 합칩니다.

nbinom1 <- rnbinom(n = 10000, p = .3, size = 1)
nbinom5 <- rnbinom(n = 10000, p = .3, size = 5)
nbinom10 <- rnbinom(n = 10000, p = .3, size = 10)
nbinom <- bind_rows(tibble(x = nbinom1, Size = 1), 
                    tibble(x = nbinom5, Size = 5), 
                    tibble(x = nbinom10, Size = 10)) %>%
  mutate(Size = as.factor(Size))

포와송이랑 이항분포에서 사용했던 함수들과 거의 유사합니다. 이렇게 만들어진 그래프는 다음과 같습니다.

ggplot(nbinom, aes(x = x)) +
  geom_density(aes(group = Size, color = Size, fill = Size), 
               adjust = 4, alpha = 1/2) +
  scale_color_manual(values = futurevisions("mars")) + 
  scale_fill_manual(values = futurevisions("mars")) + 
  labs(x = "") + theme(legend.position = "top")
Probability Mass Function of Negative Binomial Distribution

Probability Mass Function of Negative Binomial Distribution

F 분포 (F Distribution)

rf()를 통해서 F 분포에서 무작위로 10,000개의 값을 서로 다른 자유도를 이용하여 추출하여 네 개의 벡터에 저장합니다. 이제는 아시겠지만 rf()r은 random을 의미합니다.

fa <- rf(n = 10000, df1 = 1, df2 = 50)
fb <- rf(n = 10000, df1 = 5, df2 = 100)
fc <- rf(n = 10000, df1 = 50, df2 = 50)
fd <- rf(n = 10000, df1 = 50, df2 = 500)
f <- bind_rows(tibble(x = fa, DF1 = 5, DF2 = 5), 
               tibble(x = fb, DF1 = 5, DF2 = 10), 
               tibble(x = fc, DF1 = 10, DF2 = 5),
               tibble(x = fd, DF1 = 10, DF2 = 10))

이렇게 만들어진 데이터의 결과는 그림과 같습니다. 그리고 이제 x 값이 6보다 작거나 같은 경우로 하위 데이터를 만듭니다. 왜냐하면 굉장히 극단적인 수치들이 희소한 확률로 있기 때문에 그냥 그래프를 그리면 집단 별 차이를 보기가 조금 힘들 수도 있기 때문입니다. 그리고 전체 자유도를 지정해줄 것입니다. 자세한 내용은 나중에 F-distribution을 언급할 때 설명하도록 하겠습니다.

f <- f %>% dplyr::filter(x <= 6) %>%
  mutate(DF = DF1 * 100 + DF2) %>%
  mutate(DF = as.factor(DF))

ggplot(f, aes(x = x)) +
  geom_density(aes(color = DF, fill = DF), 
               adjust = 4, alpha = 1/2) +
  guides(color=guide_legend(title = "DF1, DF2")) +
  guides(fill=guide_legend(title = "DF1, DF2")) +
  scale_color_manual(
    values = futurevisions("mars"),
    labels = c("1, 50", "5, 500", "50, 50", "50, 500")) +
  scale_fill_manual(
    values = futurevisions("mars"),
    labels = c("1, 50", "5, 500", "50, 50", "50, 500"))
Probability Density Function of F Distribution

Probability Density Function of F Distribution

DF1이랑 DF2를 합쳐준다는 것은 F-분포의 특징 때문인데, 간단하게 말하자면 t-분포와 다르게 F-분포는 분자, 분모에 자유도가 하나씩 즉 2개가 필요하기 때문입니다.

오늘 포스팅한 내용 중에서는 처음 두 분포, 정규분포와 이항분포에 대해서는 숙지할 필요가 있고 나머지 분포들은 추후에 다시 다룰 기회가 있을 것이라고 생각합니다. 사실 분포는 엄청 많습니다: 웨이블, 토빗… 뭐 여럿 있기 때문에 다 외우는 것은 어렵고 여기서는 단지 R-code로 구현하는 방법을 러프하게 살펴보았다고 할 수 있을 것입니다.

그리고 ggplot()은 손에 익도록 연습하는 것이 좋습니다. R의 장점 중 하나는 여타의 통계툴들에 비해 그래프 기능이 뛰어나다는 것입니다. ggplot()R이 그러한 명성을 얻게 해 준 공신 중 하나라고 할 수 있습니다.