시뮬레이션과 재표집 - 부트스트랩과 순열 검정

이 장에서 배우는 것

앞 장에서 판매량과 설명 변수의 관계를 모형으로 표현했다. 모형이 내놓은 예측값은 하나의 수지만, 실제 판매량은 날마다 달라진다. 같은 기간에 다른 손님이 방문했다면 평균 판매량도 달라졌을 것이다. 이번에는 관측한 자료를 다시 뽑거나 가상의 수요를 여러 번 만들어, 결과가 얼마나 흔들리는지 살펴본다.

편의점 체인의 판매 자료로 평균 판매량의 불확실성을 계산하고, 두 지점의 차이를 검정하며, 발주량에 따른 재고 부족 확률을 추정한다. 계산을 반복하는 목적은 숫자를 많이 만드는 데 있지 않다. 무엇을 무작위로 바꾸었는지, 무엇을 고정했는지, 그 선택이 어떤 질문에 답하는지 밝히는 데 있다.

• set.seed()로 난수 계산을 재현하고 재현성의 범위를 설명한다.
• sample()과 replicate()로 재표집과 반복 실험을 구성한다.
• 부트스트랩으로 평균 판매량의 신뢰구간을 계산한다.
• 순열 검정으로 두 지점의 판매량 차이를 평가한다.
• 몬테카를로 시뮬레이션으로 재고 부족 확률과 계산 오차를 추정한다.

문제 상황

동네 편의점 체인이 도시락 발주량을 조정하려 한다. 한 지점의 비교 가능한 20일 판매량은 하루 18개에서 22개 사이다. 평균은 20개지만, 이 자료만으로 장기 평균을 정확히 안다고 말하기는 어렵다. 관리자는 평균 판매량의 추정 범위를 알고 싶어 한다.

별도의 소규모 시범 판매에서는 두 지점을 각각 4일 관찰했다. 판매량 차이가 보이지만, 표본이 작다. 두 지점이 같은 판매 분포를 가진다는 가정 아래에서도 이런 차이가 나타날 수 있는지 확인해야 한다. 이 예제에서는 계산 구조가 드러나도록 작은 정수 자료를 사용한다. 실제 업무에서는 상품, 영업시간, 요일, 행사 여부가 비교 가능한지 먼저 검토해야 한다.

발주 문제에는 또 다른 질문이 있다. 평균 수요가 20개라고 해서 매일 20개를 준비하면 충분한 것은 아니다. 수요가 재고를 넘는 날이 얼마나 되는지가 중요하다. 이번 예제에서는 하루 수요가 18개부터 22개까지 같은 확률로 발생한다고 가정하고, 20개를 준비할 때의 부족 확률을 계산한다. 이것은 관측 자료에서 자동으로 확정된 사실이 아니라 시뮬레이션에 넣은 가정이다.

세 질문은 서로 다른 무작위 실험으로 답한다

업무 질문	고정하는 것	바꾸는 것	계산 방법
평균 판매량의 불확실성	관측한 판매량	표본에 포함되는 관측값	부트스트랩
지점 차이의 유의성	합쳐 놓은 판매량	지점 배정	순열 검정
재고 부족 가능성	수요 분포와 재고량	가상의 하루 수요	몬테카를로

완성 프로그램은 이 세 계산을 한 파일에서 실행한다. 숫자 결과는 CSV로 저장하고, 부트스트랩 분포와 재고 부족 확률은 PNG로 저장한다. 콘솔에는 반복 횟수에 따라 달라지는 긴 숫자 목록 대신 확인할 값과 생성한 파일 이름을 출력한다.

난수와 반복 실험을 통제하기

시드는 난수열의 출발점을 정한다

컴퓨터의 난수는 정해진 규칙으로 생성된다. 난수 시드(seed)는 그 계산을 시작할 상태를 지정한다. 같은 난수 생성 방식에서 같은 시드를 설정하고 같은 순서로 함수를 호출하면 같은 결과를 얻는다. 다음 코드는 두 번의 추출이 같은지 확인한다.

set.seed(101)
first <- sample(18:22, size = 8, replace = TRUE)

set.seed(101)
second <- sample(18:22, size = 8, replace = TRUE)

stopifnot(identical(first, second))

시드는 결과를 특정 값으로 고정하는 주문이 아니다. 난수열에서 출발할 위치를 정할 뿐이다. 첫 추출 뒤에 다른 난수 함수를 호출하면 내부 상태가 이동한다. 따라서 코드 중간에 난수 호출을 추가하면 뒤의 결과도 바뀔 수 있다. 완성 코드에서는 부트스트랩과 수요 시뮬레이션에 각각 시드를 지정해 두 계산의 출발점을 분리한다.

재현성을 기록할 때는 시드뿐 아니라 R 버전과 난수 생성 방식도 함께 남기는 편이 좋다. 여기서는 RNGkind()로 난수 생성 방식, 정규 난수 생성 방식, 표본 추출 방식을 명시하고 R.version.string을 파일에 기록한다. 입력 자료와 호출 순서가 같다는 조건도 필요하다. 같은 시드를 썼다는 사실만으로 모든 운영체제와 모든 R 버전에서 파일의 바이트까지 같다고 보장하지는 않는다.

복원 추출과 비복원 추출

sample(x, size, replace)는 x에서 값을 뽑는다. replace = TRUE는 뽑은 항목을 다시 후보에 넣는 복원 추출이다. 같은 관측값이 여러 번 들어갈 수 있고, 어떤 관측값은 한 번도 들어가지 않을 수 있다. 부트스트랩은 이 방식을 사용한다.

replace = FALSE는 한 번 뽑은 항목을 다시 뽑지 않는 비복원 추출이다. 전체 길이만큼 뽑으면 순서를 섞는다. 다만 여기서 말하는 항목은 값의 종류가 아니라 벡터의 위치다. 판매량 20이 여러 위치에 있으면 비복원 추출에서도 값 20이 여러 번 나올 수 있다.

sales <- c(18, 20, 20, 22)
index <- sample.int(length(sales),
                    size = length(sales),
                    replace = TRUE)
resampled <- sales[index]

위처럼 위치를 뽑는 습관은 길이가 하나인 숫자 벡터에서도 안전하다. sample(20, 1)은 값 20만 후보로 삼는 호출이 아니라 1부터 20까지에서 하나를 뽑는 호출이다. 관측 벡터에서 뽑으려면 sample.int(length(x), ...)로 위치를 만든 뒤 x[index]로 접근하면 된다.

replicate는 실험 결과를 모은다

replicate()는 주어진 식을 반복 평가한다. 반복마다 평균 하나를 반환하면 결과는 숫자 벡터가 된다. 완성 코드의 부트스트랩 함수는 표본 추출과 평균 계산을 하나의 반복 단위로 묶는다. 함수가 반환하는 값이 무엇인지 정하면 반복 결과의 구조도 이해하기 쉬워진다.

set.seed(102)
means <- replicate(5, {
  index <- sample.int(length(sales),
                      length(sales),
                      replace = TRUE)
  mean(sales[index])
})

중괄호 안의 마지막 식이 반복 한 번의 반환값이다. 표본 전체를 반환하면 결과가 행렬로 단순화될 수 있다. 목록이 필요하면 simplify = FALSE를 사용한다. 이번에는 평균만 필요하므로 기본 동작을 유지한다.

부트스트랩과 순열 검정의 질문 구분하기

관측 자료를 임시 모집단으로 삼는다

부트스트랩(bootstrap)은 관측 자료에서 원래 표본 크기만큼 복원 추출하고, 관심 있는 통계량을 다시 계산하는 방법이다. 여기서는 20일 판매량에서 20개를 뽑고 평균을 구한다. 이 과정을 5,000번 반복하면 평균 추정량이 흔들리는 모습을 근사할 수 있다.

원자료 평균은 하나지만 재표집 평균은 여러 개다. 그 분포의 아래쪽 2.5% 지점과 위쪽 97.5% 지점을 선택하면 백분위 부트스트랩 신뢰구간을 얻는다. 코드에서는 quantile(boot_means, c(0.025, 0.975), type = 7)로 계산한다. 분위수 계산 방식을 명시해 결과를 해석할 때의 기준도 남긴다.

95% 신뢰구간은 내일 판매량의 95%가 들어가는 범위가 아니다. 평균이라는 모수를 추정하는 절차의 불확실성을 표현한다. 같은 조건에서 표본 수집과 구간 계산을 반복할 때, 적절한 가정 아래 그 구간이 실제 평균을 포함하는 비율이 약 95%가 되도록 기대하는 것이다. 이미 계산한 특정 구간에 참평균이 들어갈 확률을 곧바로 95%라고 해석하지 않는다.

백분위 방법은 간단하지만 모든 상황에서 같은 정확도를 갖지는 않는다. 표본이 매우 작거나 분포가 강하게 치우치거나 추정량에 큰 편향이 있으면 구간의 성질을 별도로 점검해야 한다. 반복 횟수를 늘리는 것은 재표집 계산의 흔들림을 줄일 뿐, 부족한 관측 자료를 늘려 주지는 않는다.

일별 판매량 사이에 강한 의존성이 있을 때도 주의해야 한다. 연속된 행사 기간이나 요일 효과가 있으면 각 날을 독립적으로 뽑는 방식이 원래 자료의 구조를 깨뜨린다. 이번 예제는 비교 가능한 독립적인 날을 관찰했다고 가정한다. 실제 자료에서는 요일별로 나누거나 연속된 날짜 묶음을 뽑는 등 자료 구조에 맞는 재표집 설계가 필요하다.

귀무가설 아래 지점 이름을 바꾼다

순열 검정(permutation test)은 관측값을 다시 생성하는 대신 집단 배정을 바꾼다. 지점 A와 B의 판매량을 합친 뒤 같은 크기의 두 집단으로 재배정한다. 검정 통계량은 A 평균에서 B 평균을 뺀 값으로 정한다. 양측 검정에서는 그 절댓값을 비교한다.

이 예제의 귀무가설은 두 지점 자료가 같은 분포에서 나와 지점 표지를 교환할 수 있다는 것이다. 이를 교환 가능성이라고 한다. 단지 평균이 같다는 조건만으로 서로 다른 분산이나 다른 관측 구조를 무시해도 된다는 뜻은 아니다. 날짜가 짝지어진 자료라면 짝을 유지하는 재배정이 필요하고, 관찰 자료에서 지점의 운영 조건이 다르면 검정 결과만으로 지점 자체의 인과 효과를 주장할 수 없다.

지점별 표본이 4개이므로 합쳐진 8개 중 A에 들어갈 4개 위치를 고르는 방법은 70개다. 이번에는 무작위로 일부만 뽑지 않고 combn()으로 70개를 모두 계산한다. 집단 내부의 순서는 평균 차이에 영향을 주지 않으므로 서로 다른 집단 배정만 열거하면 충분하다.

자료는 A가 18, 19, 20, 21이고 B가 22, 23, 24, 25다. 관측 평균 차이는 −4다. 70개 배정 중 절댓값이 4 이상인 것은 가장 작은 네 값이 A가 되는 경우와 가장 큰 네 값이 A가 되는 경우, 두 개다. 따라서 정확한 양측 p값은 2/70이다. p값은 귀무가설이 참일 확률이 아니라, 귀무가설에 따른 재배정에서 관측 통계량만큼 또는 더 극단적인 결과가 나오는 비율이다.

배정이 너무 많으면 무작위 순열로 근사할 수 있다. 그때는 관측 배정을 포함하는 방식으로 (extreme + 1) / (B + 1)을 쓰는 구성을 사용할 수 있다. 완성 코드는 전수 열거이므로 mean()으로 비율을 계산하며 이 보정을 추가하지 않는다.

몬테카를로로 재고 부족 확률 추정하기

몬테카를로(Monte Carlo) 방법은 가정한 확률 모형에서 실험을 반복해 관심 있는 양을 근사한다. 이번에는 18개부터 22개까지의 정수 수요가 각각 같은 확률로 발생한다고 가정한다. 재고가 20개라면 수요 21개와 22개인 날에 부족하다.

부족 사건은 demand > stock이다. 수요와 재고가 모두 20이면 모든 수요를 충족했으므로 부족으로 세지 않는다. 다만 하루가 끝날 때 재고가 없는 사건을 분석하려면 다른 정의가 필요하다. 업무 용어를 비교 연산자로 바꾸는 순간, 분석하려는 사건이 결정된다.

가상의 수요를 100,000번 생성하고 부족 여부를 논리 벡터로 만든다. R에서 논리값의 평균은 참의 비율이므로 mean(demand > stock)이 부족 확률 추정값이다. 이 모형에서는 가능한 다섯 수요 중 두 가지가 부족하므로 이론값은 0.4다. 계산 결과가 이 값 근처인지 검사할 수 있다.

독립적으로 생성한 부족 여부는 0 또는 1을 갖는다. 추정 확률을 p̂, 반복 횟수를 M이라고 하면 몬테카를로 표준오차는 대략 √{p̂(1−p̂)/M}이다. 이 값은 반복 계산 때문에 생기는 오차를 나타낸다. 수요 분포를 잘못 가정한 데서 생기는 오차나, 실제 판매 자료가 적어서 생기는 불확실성은 포함하지 않는다.

표준오차를 절반으로 줄이려면 반복 횟수를 대략 네 배로 늘려야 한다. 이미 충분히 작은 계산 오차를 더 줄이기 전에 수요 모형이 현실에 맞는지 살펴보는 편이 유용하다. 특히 판매량은 재고가 충분할 때만 수요를 온전히 보여 준다. 품절된 날의 판매량을 그대로 수요로 쓰면 부족 위험을 작게 평가할 수 있다.

함수의 동작과 인수는 사실 확인이 필요할 때 공식 문서에서 확인할 수 있다. 난수 상태와 생성 방식, 표본 추출, 분위수 계산 문서는 이 장의 계산에서 사용한 기능을 확인하는 근거다.

완성 코드

다음 내용을 main.R로 저장한다. 실행 위치에 파일을 쓸 권한이 있어야 한다. 외부 자료나 패키지는 필요하지 않다. 숫자 결과는 충분한 자릿수로 CSV에 저장하고, 실행 환경은 별도 텍스트 파일에 남긴다. 그래프에는 한글 글꼴 설치 여부에 따른 차이를 피하기 위해 영문 축 이름을 사용한다.

options(warn = 2)

RNGkind(kind = "Mersenne-Twister",
        normal.kind = "Inversion",
        sample.kind = "Rejection")

check <- function(ok, message) {
  if (!isTRUE(ok)) {
    stop(message, call. = FALSE)
  }
  invisible(TRUE)
}

bootstrap_mean <- function(x, times) {
  check(is.numeric(x) && length(x) >= 2L,
        "x must contain at least two numeric values")
  check(all(is.finite(x)), "x must contain finite values")
  check(length(times) == 1L && is.finite(times) &&
          times >= 2 && times == floor(times),
        "times must be an integer of at least two")

  replicate(times, {
    index <- sample.int(length(x), length(x), replace = TRUE)
    mean(x[index])
  })
}

sales <- rep(18:22, times = 4L)
boot_times <- 5000L

set.seed(20261001)
boot_means <- bootstrap_mean(sales, boot_times)

set.seed(20261001)
boot_again <- bootstrap_mean(sales, boot_times)
same_bootstrap <- identical(boot_means, boot_again)
check(same_bootstrap, "bootstrap reproduction failed")

boot_ci <- unname(
  quantile(boot_means, probs = c(0.025, 0.975), type = 7)
)
check(length(boot_means) == boot_times,
      "incorrect bootstrap length")
check(all(boot_means >= min(sales) &
          boot_means <= max(sales)),
      "bootstrap mean outside data range")
check(boot_ci[1L] <= boot_ci[2L], "invalid interval order")

branch_a <- c(18, 19, 20, 21)
branch_b <- c(22, 23, 24, 25)
pooled <- c(branch_a, branch_b)
observed_diff <- mean(branch_a) - mean(branch_b)

assignments <- combn(seq_along(pooled), length(branch_a))
perm_diffs <- apply(assignments, 2L, function(index) {
  mean(pooled[index]) - mean(pooled[-index])
})
tolerance <- 1e-12
perm_p <- mean(
  abs(perm_diffs) >= abs(observed_diff) - tolerance
)
stopifnot(ncol(assignments) == 70L)
check(abs(perm_p - 2 / 70) < tolerance,
      "exact permutation check failed")

stock <- 20L
mc_times <- 100000L
demand_values <- 18:22

set.seed(20261002)
demand <- sample(demand_values, mc_times, replace = TRUE)
shortage <- demand > stock
shortage_prob <- mean(shortage)
mc_se <- sqrt(shortage_prob * (1 - shortage_prob) / mc_times)
theoretical_prob <- mean(demand_values > stock)

check(length(shortage) == mc_times, "incorrect simulation length")
check(is.finite(mc_se) && mc_se > 0,
      "invalid Monte Carlo standard error")

stock_grid <- 18:22
stock_probs <- vapply(stock_grid, function(s) {
  mean(demand > s)
}, numeric(1))
check(all(diff(stock_probs) <= 0),
      "shortage probability must not increase with stock")
check(stock_probs[length(stock_probs)] == 0,
      "maximum stock must cover all simulated demands")

results <- data.frame(
  metric = c("observed_mean", "bootstrap_ci_lower",
             "bootstrap_ci_upper", "observed_difference",
             "exact_permutation_p", "shortage_probability",
             "monte_carlo_se", "theoretical_shortage_probability"),
  value = c(mean(sales), boot_ci, observed_diff, perm_p,
            shortage_prob, mc_se, theoretical_prob)
)
stock_results <- data.frame(
  stock = stock_grid,
  shortage_probability = stock_probs
)

options(digits = 17)
write.csv(results, "simulation_results.csv", row.names = FALSE)
write.csv(stock_results, "stock_risk.csv", row.names = FALSE)
writeLines(
  c(R.version.string,
    paste("RNG:", paste(RNGkind(), collapse = ", ")),
    "Bootstrap seed: 20261001",
    "Demand seed: 20261002",
    paste("Bootstrap repetitions:", boot_times),
    paste("Monte Carlo repetitions:", mc_times)),
  "run_info.txt"
)

png("bootstrap.png", width = 960, height = 640)
hist(boot_means, breaks = seq(18, 22, by = 0.1),
     col = "#e1edf7", border = "#315675",
     main = "Bootstrap means",
     xlab = "Mean daily sales", ylab = "Frequency")
abline(v = mean(sales), col = "#173047", lwd = 2)
abline(v = boot_ci, col = "#c0392b", lwd = 2, lty = 2)
legend("topright",
       legend = c("Observed mean", "95% percentile interval"),
       col = c("#173047", "#c0392b"),
       lty = c(1, 2), lwd = 2, bty = "n")
invisible(dev.off())

png("shortage.png", width = 960, height = 640)
plot(stock_grid, stock_probs, type = "b", pch = 19,
     ylim = c(0, 1), xaxt = "n",
     col = "#315675", lwd = 2,
     main = "Stock and shortage probability",
     xlab = "Stock", ylab = "Estimated probability")
axis(1, at = stock_grid)
points(stock, shortage_prob, col = "#c0392b", pch = 19, cex = 1.5)
invisible(dev.off())

cat(sprintf("Observed mean: %.1f\n", mean(sales)))
cat(sprintf("Observed difference (A - B): %.1f\n", observed_diff))
cat(sprintf("Exact permutation p-value: %.4f\n", perm_p))
cat(sprintf("Theoretical shortage probability: %.1f\n",
            theoretical_prob))
cat("Bootstrap reproducible:", same_bootstrap, "\n")
cat("Checks: passed\n")
cat("Saved: simulation_results.csv\n")
cat("Saved: stock_risk.csv\n")
cat("Saved: run_info.txt\n")
cat("Saved: bootstrap.png\n")
cat("Saved: shortage.png\n")

줄별 해설

options(warn = 2)는 경고가 발생하면 오류로 바꾼다. 이 프로그램에서 경고를 숨기지 않고 실행을 멈추게 하는 확인 장치다. RNGkind()는 난수 계산 방식을 명시한다. R은 이 스크립트를 읽어 실행하므로 별도의 컴파일 명령은 필요하지 않다.

check()는 검사 조건이 하나의 참인지 확인한다. 조건이 거짓이거나 결측값이면 메시지와 함께 멈춘다. bootstrap_mean()의 첫 검사들은 빈 자료, 결측값, 무한대, 잘못된 반복 횟수가 계산에 들어가지 않도록 한다. 결측값을 조용히 제거하지 않고 자료를 먼저 확인하게 만든 선택이다.

sample.int()는 자료의 위치를 원래 표본 크기만큼 복원 추출한다. mean(x[index])는 그 표본의 평균 하나를 반환한다. replicate()가 이를 모아 길이 5,000의 숫자 벡터를 만든다. 원자료 sales는 다섯 판매량을 네 번 반복한 20개 관측값이며 평균은 20이다.

같은 시드를 다시 설정한 뒤 boot_again을 계산하는 부분은 재현성 검사다. 두 벡터의 모든 값과 구조가 같은지 identical()로 확인한다. 반복 결과를 두 번 만드는 것은 여기서는 학습과 검사를 위한 비용이다. 실제 분석에서는 검사를 별도 실행 흐름으로 옮길 수 있다.

quantile()은 두 분위수에 이름을 붙여 반환한다. unname()으로 이름을 제거하면 이후 결과 표를 만들 때 값의 위치가 명확해진다. 평균이 원자료 최솟값과 최댓값 사이인지 확인하는 검사는 시드에 의존하는 특정 숫자 대신 계산의 성질을 검증한다.

combn()의 각 열은 A에 배정할 네 위치다. apply(..., 2L, ...)는 열마다 평균 차이를 계산한다. 양의 위치 벡터는 A를 선택하고 pooled[-index]는 나머지를 B로 선택한다. 값이 중복되는 자료에서도 위치를 기준으로 배정을 세어야 한다.

순열 p값 계산은 관측 차이와 같은 크기의 결과도 포함한다. 작은 tolerance는 소수 계산의 미세한 차이로 경계값이 빠지는 것을 줄이기 위한 것이다. 허용 오차는 통계량의 규모에 맞춰 정해야 한다. 이 예제의 평균 차이는 정확히 표현 가능한 값이므로 이 보정이 결과를 바꾸지는 않는다.

수요 생성에는 부트스트랩과 다른 시드를 사용한다. shortage의 평균은 부족 확률이고, mc_se는 그 추정의 몬테카를로 표준오차다. theoretical_prob는 같은 확률의 다섯 수요를 직접 세어 얻은 비교 기준이다. 추정값이 이론값과 정확히 같아야 한다는 검사는 하지 않는다. 유한한 무작위 실험에는 오차가 있기 때문이다.

vapply()는 여러 재고량에 대해 숫자 하나씩을 받도록 결과 형태를 지정한다. 모든 재고량에 같은 수요 벡터를 사용하므로 재고가 늘어날 때 추정 부족 확률이 증가하지 않는다. 이 성질과 최대 수요를 모두 충족하는 재고에서 부족 확률이 0이라는 성질을 검사한다.

CSV의 value 열에는 신뢰구간 경계와 추정 확률을 그대로 저장한다. options(digits = 17)은 숫자를 파일에 쓸 때 사용할 표시 정밀도를 설정한다. PNG를 열고 그림을 그린 뒤 dev.off()로 장치를 닫아 파일을 완성한다. invisible()은 장치를 닫을 때의 반환값을 출력하지 않게 한다.

실행 결과

터미널에서 파일을 저장한 디렉터리로 이동한 뒤 실행한다. 다음 콘솔 출력은 코드의 출력 형식과 일치한다. 부트스트랩 구간과 몬테카를로 추정값의 상세 숫자는 simulation_results.csv에 들어 있다.

Rscript main.R
Observed mean: 20.0
Observed difference (A - B): -4.0
Exact permutation p-value: 0.0286
Theoretical shortage probability: 0.4
Bootstrap reproducible: TRUE
Checks: passed
Saved: simulation_results.csv
Saved: stock_risk.csv
Saved: run_info.txt
Saved: bootstrap.png
Saved: shortage.png

bootstrap.png에서는 재표집 평균이 20 근처에 모이며 붉은 점선 두 개가 백분위 신뢰구간의 경계를 표시한다. 이 그림의 가로축은 개별 날짜의 판매량이 아니라 재표집 평균이다. shortage.png에서는 재고가 늘어날수록 부족 확률이 감소한다. 재고 22개에서는 가정한 수요를 모두 충족하므로 확률이 0이다.

순열 p값 0.0286은 이 예제의 교환 가능성 가정 아래 관측 차이가 드문 편이라는 근거다. 실제 관리 판단에는 평균 차이의 크기와 운영 조건도 함께 필요하다. 재고 부족 추정값이 0.4 근처라고 해도 현실의 부족 확률을 확인한 것은 아니다. 먼저 같은 확률의 다섯 수요라는 모형을 받아들였다는 조건이 붙는다.

실무에서 자주 틀리는 것

반복 안에서 같은 시드를 다시 설정한다

반복마다 같은 시드를 설정하면 같은 표본을 계속 뽑는다. 평균 벡터가 길어도 서로 다른 실험을 수행한 것이 아니다.

# 틀린 코드
bad <- replicate(5000, {
  set.seed(101)
  mean(sample(sales, length(sales), replace = TRUE))
})

# 고친 코드
set.seed(101)
good <- bootstrap_mean(sales, 5000L)

시드는 반복을 시작하기 전에 설정한다. 재현성을 검사할 때만 전체 반복의 출발점으로 다시 돌아간다.

부트스트랩에서 복원 추출을 빠뜨린다

원래 길이만큼 비복원 추출하면 값의 순서만 바뀐다. 평균은 매번 같아져 불확실성이 사라진 것처럼 보인다.

# 틀린 코드
bad <- replicate(5000, {
  mean(sample(sales, length(sales), replace = FALSE))
})

# 고친 코드
set.seed(102)
good <- replicate(5000, {
  index <- sample.int(length(sales), length(sales),
                      replace = TRUE)
  mean(sales[index])
})

추출 방법은 반복 횟수보다 먼저 점검해야 한다. 잘못된 실험을 많이 반복해도 질문에 맞는 결과가 되지 않는다.

양측 검정에서 한쪽 꼬리만 센다

관측 차이가 음수일 때 perm_diffs >= observed_diff만 계산하면 반대 방향의 큰 차이를 적절히 평가하지 못한다.

# 틀린 코드
bad_p <- mean(perm_diffs >= observed_diff)

# 고친 코드
good_p <- mean(
  abs(perm_diffs) >= abs(observed_diff) - tolerance
)

여기서는 평균 차이의 절댓값을 극단성 기준으로 정한 양측 검정이다. 방향이 있는 업무 질문이라면 자료를 보기 전에 단측 검정의 방향을 정해야 한다. 모든 통계량의 양측 검정이 같은 절댓값 규칙을 쓰는 것은 아니다.

재고 소진과 재고 부족을 섞는다

수요가 재고와 같은 날에는 남은 재고가 없지만 미충족 수요도 없다. 부족 확률을 구하면서 같은 값까지 포함하면 다른 사건을 계산하게 된다.

# 틀린 코드: 부족 사건에 같은 값까지 포함한다
bad_prob <- mean(demand >= stock)

# 고친 코드: 충족하지 못한 수요가 있는 날을 센다
good_prob <- mean(demand > stock)

# 별도 지표: 부족한 수량의 하루 평균
mean_unmet <- mean(pmax(demand - stock, 0))

부족한 날의 비율과 부족 수량의 평균도 서로 다르다. 발주 결정을 위해서는 폐기 비용, 보관 제약, 부족 수량까지 추가로 고려할 수 있지만 이번 계산의 목표는 사건의 확률을 구하는 데 있다.

한눈에 보기

함수와 결과의 의미를 함께 기억한다

도구	이 장의 역할	결과	해석할 때의 조건
set.seed()	난수 상태 지정	반복 가능한 계산	생성 방식과 호출 순서도 같아야 한다
sample.int()	관측 위치 복원 추출	재표집 표본	자료의 의존 구조를 고려한다
replicate()	통계량 반복 계산	평균 벡터	반복 안에서 시드를 초기화하지 않는다
quantile()	백분위 구간 계산	평균의 신뢰구간	개별 날짜의 예측구간과 구분한다
combn()	집단 배정 전수 열거	정확한 순열 p값	귀무가설 아래 교환 가능성이 필요하다
논리값의 mean()	부족 사건의 비율 계산	확률 추정값	정한 수요 모형에 조건부인 결과다

다른 코드 묶음과 연결할 때도 계산의 역할은 같다. 재표집 함수는 숫자 자료와 반복 횟수를 받고 통계량 벡터를 돌려준다. 저장과 그래프는 그 결과를 이용한다. 다음에는 이처럼 역할이 드러나는 함수를 묶고 검사하는 작업으로 이어진다.

다른 도구를 만났을 때 기능을 대응시키는 기준

현재 사용한 도구	다른 도구의 대응	공통 목적
vapply()	purrr의 map_dbl()	입력별 숫자 결과 하나를 모은다
data.frame()	tibble의 tibble()	결과를 열 단위로 정리한다
check(), stopifnot()	testthat의 기대값 검사	계산이 지켜야 할 조건을 확인한다

대응표의 도구는 이 프로그램의 실행에 필요하지 않다. 자료 구조와 검사의 의미를 이해한 뒤 같은 목적의 다른 표현을 읽는 데 사용하면 된다.

연습 문제

1. 평균 대신 중앙값의 부트스트랩 신뢰구간을 계산하라. 재표집 표본 크기와 복원 추출 조건을 유지하고 시드를 지정하라.
2. 완성 코드의 순열 통계량을 B 평균에서 A 평균을 뺀 값으로 바꾸라. 관측 차이의 부호와 양측 p값이 어떻게 변하는지 설명하라.
3. 같은 수요 벡터를 사용해 재고 19개와 21개의 부족 확률 및 평균 부족 수량을 계산하라. 각각의 이론값도 구하라.
4. 몬테카를로 반복 횟수를 100,000번에서 400,000번으로 늘리면 표준오차가 대략 어떻게 변하는지 설명하라. 이 변화가 수요 모형의 정확성을 높이는지도 답하라.

정답과 해설

1. 표본 추출은 그대로 두고 마지막 통계량을 median()으로 바꾼다. 중앙값은 가능한 값이 제한되어 재표집 분포에 같은 값이 많이 나타날 수 있다.

   set.seed(20261003)
   boot_medians <- replicate(5000L, {
     index <- sample.int(length(sales), length(sales),
                         replace = TRUE)
     median(sales[index])
   })
   median_ci <- quantile(boot_medians, c(0.025, 0.975), type = 7)
   print(median_ci)

   반복마다 원래 자료와 같은 개수를 뽑는 이유는 같은 표본 크기에서 중앙값이 어떻게 변하는지 살펴보기 위해서다.

2. 새 관측 차이는 4다. 모든 순열 통계량의 부호도 반대로 바뀐다. 절댓값은 같으므로 양측 p값은 2/70으로 유지된다.

   reversed_observed <- mean(branch_b) - mean(branch_a)
   reversed_permutations <- -perm_diffs
   reversed_p <- mean(
     abs(reversed_permutations) >=
       abs(reversed_observed) - tolerance
   )
   stopifnot(abs(reversed_p - perm_p) < tolerance)

   단측 검정이라면 비교 방향도 함께 검토해야 한다. 이름의 순서를 바꾸는 것과 업무 가설의 방향을 바꾸는 것은 구분해야 한다.

3. 재고 19개에서는 수요 20, 21, 22가 부족하므로 이론적 부족 확률은 0.6이다. 부족 수량은 다섯 수요에서 각각 0, 0, 1, 2, 3이므로 평균은 1.2개다. 재고 21개에서는 수요 22만 부족하므로 확률은 0.2이고 평균 부족 수량은 0.2개다.

   exercise_stock <- c(19L, 21L)
   exercise_result <- data.frame(
     stock = exercise_stock,
     probability = vapply(exercise_stock, function(s) {
       mean(demand > s)
     }, numeric(1)),
     mean_unmet = vapply(exercise_stock, function(s) {
       mean(pmax(demand - s, 0))
     }, numeric(1))
   )
   print(exercise_result)

   시뮬레이션 값은 이론값 근처에 나타나며 유한한 반복 때문에 차이가 생긴다. 같은 수요를 사용하면 두 재고 정책의 비교에서 가상 날짜 구성이 같아진다.

4. 표준오차는 반복 횟수의 제곱근에 반비례하므로 대략 절반이 된다. 추정 확률도 조금 달라질 수 있어 실제 비율이 정확히 절반일 필요는 없다.

   set.seed(20261004)
   more_demand <- sample(demand_values, 400000L, replace = TRUE)
   more_prob <- mean(more_demand > stock)
   more_se <- sqrt(more_prob * (1 - more_prob) / 400000L)
   print(c(probability = more_prob, standard_error = more_se))

   반복을 늘리면 주어진 모형 안에서 계산이 더 정밀해진다. 수요가 실제로 균등한지, 요일에 따라 달라지는지, 품절 때문에 관측이 제한되었는지는 별도의 자료와 검토로 판단해야 한다.