Modern Data Science with R_Cpt10

10. Simulation

10.1. Reasoning in reverse

지금까지 데이터는 이해를 돕기 위한 것이라고 배워왔다. 하지만 10챕터에서는 반대로 추측 >>> 데이터의 순서로 유용하게 데이터를 재구성할 것이다. 재구성된 데이터는 해당 시스템에서 생성된 원 데이터에 대한 정보를 준다.


재구성 방법:
1. Conditional inference, 조건부 추론
: 생성한 메커니즘이 실제 시스템 작동 방식을 반영한다면 생성되는 데이터도 실제 데이터와 유사할 것이다.

2. Winnowing out hypotheses, 가설 선별하기
: less desirable한 선택을 제거하여서 남은 것이 유용하도록 만드는 것.
가설을 통해 데이터를 생성한 뒤 실제 데이터와 유사하지 않은 데이터의 가설을 제거한다.

“Making up”으로 표현하는 것은 사기에 가깝다. “Similar”에서 파생된 Simulation을 사용하는 것이 옳다.

10.2. Randomizing functions

: 가능한 가설만큼 가능한 시뮬레이션도 많다. 시뮬레이션에는 많은 기법이 있는데, 먼저 Shuffling은 가장 일반적인 무작위 방법이다. 셔플링은 순서를 무너뜨리고 우연에 의한 순서의 모양만 남기는 것을 말한다. 우리는 시뮬레이션에서 명시적으로 지정하는 것 외의 다른 구조가 들어가는 것을 원하지 않기 때문에 난수를 사용하여 시뮬레이션에서 찾은 구조가 시뮬레이션을 위한 메커니즘 때문이거나 우연적인 것임을 확인할 수 있다. R에서 가장 널리 쓰이는 난수생성기는 runif()다.

runif(5)

[1] 0.4763697 0.4985587 0.2567427 0.4916694 0.1174759

이는 굉장히 단편적인 예시일 뿐이며, 다음과 같이 난수생성기를 통해 다른 난수화 장치를 구성할 수도 있다.

select_one <- function(vec) {
n <- length(vec)
ind <- which.max(runif(n))
vec[ind]
}

select_one(letters)

[1] "b"

select_one(letters)

[1] "u"

데이터 분석에서 널리 사용되는 시뮬레이션은 sample(), resample() 그리고 shuffle() 함수일 것이다. 이 함수들은 데이터프레임이나 벡터를 바꾸지 않고 균일하게 무작위로 샘플링하거나 행을 치환한다. 그 밖의 시뮬레이션 구축에는 중요한 속성을 가진 난수를 생성하는 기능도 중요하다.

1. rnorm(): 정규 분포
2. rexp(): 지수 분포
3. rpois(): 포아송 분포

실제 상황에 필요한 메커니즘에 맞는 함수를 적절히 사용해야 한다.

10.3. Simulating variability

10.3.1. The partially planned rendezvous, 부분적으로 계획된 만남 (예시)

Sally와 Joan이 같이 공부하기 위해 저녁 7~8시 사이에 만나기로 하였다. 두 사람은 10분만 기다렸다가 떠나는 급한 성격의 소유자이다. 도착 시간을 정하거나, 두 사람이 만날 수 있는 확률을 구하시오.

두 사람 모두 1시간 내에 도착할 확률은 동일하므로 균일한 난수를 생성한 뒤 숫자 쌍마다 시차라 10분 이하인 데이터를 확인한다.

n <- 100000

sim_meet <- data.frame(sally <- runif(n, min = 0, max = 60),
                       joan <- runif(n, min = 0, max = 60)) %>%
  mutate(result = ifelse(abs(sally - joan) <= 10,
                         "They meet", "They do not"))

tally(~ result, format = "percent", data = sim_meet)

result
They do not   They meet 
     69.494      30.506 

binom.test(~result, n, success = "They meet", data = sim_meet)

data:  sim_meet$result  [with success = They meet]
number of successes = 30506, number of trials = 1e+05, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.3022069 0.3079244
sample estimates:
probability of success 
               0.30506 

시뮬레이션을 100,000회 시행하였다. 결과에 따르면 두 사람이 만날 확률은 약 30%이며 신뢰구간이 매우 좁으므로 주어진 문제의 답으로 사용하기 충분해보인다.

10.3.2. EX: The jobs report, 일자리 보고서

일자리 보고서가 ‘가우스 분포’, ’표준편차’등의 통계적 단어로 구성되어있다. 상사가 보고서 그 자체를 경제적 지표로 심각하게 받아들이지 않도록 시뮬레이션을 통해 이해시켜 보아라.

jobs_true <- 150
jobs_se <- 65

gen_samp <- function(true_mean, true_sd, num_months = 12, delta = 0, id = 1) {
  samp_year <- rep(true_mean, num_months) +
    rnorm(num_months, mean = delta * (1:num_months), sd = true_sd)
  return(data.frame(jobs_number = samp_year, month = as.factor(1:num_months), id = id))
}

n_sims <- 3
params <- data.frame(sd = c(0, rep(jobs_se, n_sims)),
                     id = c("Truth", paste("Sample", 1:n_sims)))

df <- params %>%
  group_by(id) %>%
  dplyr::do(gen_samp(true_mean = jobs_true, true_sd = .$sd, id = .$id))

ggplot(data = df, aes(x = month, y = jobs_number)) +
  geom_hline(yintercept = jobs_true, linetype = 2) +
  geom_bar(stat = "identity") +
  facet_wrap(~ id) + ylab("Number of new jobs (in thousands)")

실제 데이터(Truth)는 일자리가 일정하지만, 시뮬레이션으로 만든 샘플은 다르게 해석될 수도 있다. 추론을 통해 결론을 내리기 전에 시스템의 변동성을 이해해야 한다.

10.3.3. EX: Restaurant health and sanitation grades

뉴욕 시의 식당 건강위반 검출을 위해 1년에 한번 불시 검문을 하여 다양한 기준을 통해 등급을 매긴다. 등급을 나누는 현재의 점수 기준이 적절한 지 시뮬레이션을 통해 판단해보아라.(A등급=0에서 13점, B등급=14에서 27점, C등급=28점 이상)

minval <- 7
maxval <- 19

JustScores <- Violations %>%
  filter(score >= minval & score <= maxval) %>%
  select(dba, score) %>%
  unique()

ggplot(data = JustScores, aes(x = score)) +
geom_histogram(binwidth = 0.5) + geom_vline(xintercept = 13, linetype = 2) +
scale_x_continuous(breaks = minval:maxval) +
annotate("text", x = 10.5, y = 10300,
label = "A grade: score of 13 or less")

13점에서 14점으로 식당의 수가 약 5,000개 급감하는 모습이다. 이정도 수치가 단순한 우연으로 발생할 확률이 매우 낮으므로, 등급 측정 시 13점~14점 사이인 경우의 등급컷을 14로 낮춰선 안된다.

10.4. Simulating a complex system

시뮬레이션은 복잡한 시스템을 이애하는 데 유용하다. 다음 예시를 통해 알아보자.

EX: 창구 직원이 1명인 은행이 있다. 첫 고객이 9시 2분에 도착하고 서비스에 5분이 소요된다. 두 번째 고객이 9시 5분에 도착하고 서비스에 3분이 소요된다. 세 번째 고객이 9시 8분에 도착하고 서비스에 2분이 소요된다. 본 시스템을 시뮬레이션으로 표현하되, 고객 수는 포아송 분포, 서비스 처리시간은 지수 분포로 가정하여라.

any_active <- function(df) {

Warning message:
In overlay(...) :
  reverting 'unnamed.chunk.label' to 'unnamed-chunk' for duration of render

  return(max(df$endtime) == Inf)
}

next_customer <- function(df) {
  res <- filter(df, endtime == Inf) %>%
    arrange(arrival)
  return(head(res, 1))
}

update_customer <- function(df, cust_num, end_time) {
  return(mutate(df, endtime = ifelse(custnum == cust_num, end_time, endtime)))
}

시뮬레이션을 실행하는 함수를 정의하였다.

run_sim <- function(n = 1/2, m = 3/2, hours = 6) {
  customers <- rpois(hours * 60, lambda = n)
  arrival <- numeric(sum(customers))
  position <- 1
  for (i in 1:length(customers)) {
    numcust <- customers[i]
    if (numcust != 0) {
      arrival[position:(position + numcust - 1)] <- rep(i, numcust)
      position <- position + numcust
    }
  }
  duration <- rexp(length(arrival), rate = 1/m)
  df <- data.frame(arrival, duration, custnum = 1:length(duration),
                   endtime = Inf, stringsAsFactors = FALSE)
  endtime <- 0
  while (any_active(df)) {
    next_one <- next_customer(df)
    now <- ifelse(next_one$arrival >= endtime, next_one$arrival, endtime)
    endtime <- now + next_one$duration
    df <- update_customer(df, next_one$custnum, endtime)
  }
  df <- mutate(df, totaltime = endtime - arrival)
  return(favstats(~ totaltime, data = df))
}

sim_results <- do(3) * run_sim()
sim_results

결과를 보면 6시간 동안 고객의 수는 n=189~210명 사이였고 고객이 가장 많은 날의 지연시간(mean)이 가장 컸다. 시뮬레이션을 더 많이 수행하면 더 합리적인 근사치를 얻을 수 있을 것이다.

10.5. Key principles of simulation

1. Design: 시뮬레이션을 설계할 문제의 모든 측면을 고려하고 탐색할 가정과 시나리오를 결정하는 것을 말한다. 문제와 관련된 적절한 상황을 잘 반영해야 하며, 간단한 설정에서 복잡성을 추가하는 경우가 많다.
   
   
2. Modularity: 다른 옵션과 매개 변수로 반복하여 실행할 수 있게 시뮬레이션 구성 함수를 짜는 것을 말한다.
   
   
3. Reproducibility and random number seeds: 시뮬레이션은 재현가능하고 대표성이 있어야 함을 말한다. 불필요한 계산을 피하면서 결과가 무작위로 변동하지 않도록 균형을 잡아야 한다. random number seeds는 시뮬레이션 과정에서 동일한 난수 목록을 사용하게 하는 것이다. set.seed()함수를 사용하여 난수 목록을 고정하며, 이후 실행될 때마다 결과가 동일하게 출력된다.