아이스 커피 주문량 예측
Kong, Seokkyu
2022-06-04
Executive Summary
나는 아이스 커피를 판매하는 가상의 가게를 운영하고 있다. 더운 여름에는 최고 기온에 따른 아이스 커피의 판매량이 달라진다.
해결하고자 하는 문제는 최고 기온에 따른 아이스 커피 판매량을 예측하고 싶다. 이를 위헤 단순 선형 회귀식을 모델링하며 베이지안 추론 방식으로 접근한다.
데이터는 ‘만화로 쉽게 배우는 회귀분석’ ([1] page 66) 의 최고 기온에 따른 주문 수를 참고로 한다. 온도는 섭씨를 따른다.
Introduction
더운 여름날, 아이스 커피의 판매량은 최고 기온와 상관관계가 있다고 생각한다. 최고 기온이 높을 수록 사람들은 아이스 커피를 더욱 많이 찾게 된다.
데이터 수집 계획 및 방법으로 2주 동안 일자별 최고 기온과 아이스 커피 판매량을 기록한다. 최고 기온이 높을수록 아이스 커피 판매량은 선형적으로 증가한다고 가정한다. 요소가 하나인 단순 선형 회귀식을 모델로 작성한다.
선형 회귀식을 모델링한 이후, 최고 기온에 따른 판매량을 예측한다.
Data
데이터는 2주 동안 최고 기온과 주문 수량을 일자별로 기록한다. 데이터를 취득하는 과정에서 결측치나 누락치는 없다.
dat <- read.csv('dat_icetea_order.csv')
# 데이터의 앞부분 확인
head(dat)## date temperature orders
## 1 22 29 77
## 2 23 28 62
## 3 24 34 93
## 4 25 31 84
## 5 26 25 59
## 6 27 29 64summary(dat)## date temperature orders
## Min. : 1.00 Min. :24.00 Min. :51.00
## 1st Qu.: 8.50 1st Qu.:26.50 1st Qu.:62.50
## Median :24.50 Median :29.50 Median :74.00
## Mean :19.64 Mean :29.14 Mean :72.57
## 3rd Qu.:27.75 3rd Qu.:31.00 3rd Qu.:83.00
## Max. :31.00 Max. :34.00 Max. :93.00# 최고 기온에 따른 주문 수량 확인
plot(orders ~ temperature, data = dat)
Model
최고 기온에 따른 주문 수량에 대한 데이터 플롯을 보면, 선형 회귀식 모델로 충분히 설명 가능하다고 본다.
추정하고자 하는 선형 회귀 모델은 \(y_i = \beta_0 + \beta_1 x_1 + e_i, e_i \sim N(0, \alpha^2), i = 1,...,n\) 이며, 아래와 같이 우도(likelihood) 및 계수에 대한 사전분포(prior)로 표현할 수 있다.
\(y_i | x_1, \beta, \alpha^2 \sim N(\beta_0 + \beta_1 x_1, \alpha^2)\)
\(\beta_i \sim N(0, 1e6)\)
\(\alpha^2 \sim ~ Gamma(a, b)\)
\(y_i\) 는 최고 기온 \(x_1\), \(\beta_0\), \(\beta_1\), \(\alpha^2\) 이 주어졌을 때, 정규분포를 따른다고 가정한다.
이 모델의 \(y_i\) 는 \(E(\beta_0 + \beta_1 x_1)\) 이 되며, \(\beta_1\) 계수는 최고 기온이 1도 오를 때 주문량 기대치의 증감량이 된다고 해석할 수 있다. \(\beta_0\) 는 절편이며, \(\beta_1\) 이 0 일 때의 값이 된다. 하지만, 이것을 최고 기온이 0 도 일경우의 판매량이라고 해석하기는 어렵다.
위와 같은 선형 회귀 모델을 베이지안으로 검토하면 아래와 같다.
library(rjags) # load library
mode1_string = " model {
# likelihood 구현
for (i in 1:n) {
y[i] ~ dnorm(mu[i], prec)
mu[i] = b[1] + b[2] * temperature[i]
}
# 계수에 대한 사전분포
for (i in 1:2) {
b[i] ~ dnorm(0.0, 1.0/1.0e6)
}
# 분산에 대한 사전분포는 역감마분포가 된다.
# 따라서 정밀도(=1/분산)는 감마분포를 따른다.
prec ~ dgamma(5/2.0, 5*10.0/2.0)
sig2 = 1.0 / prec
sig = sqrt(sig2)
} "
set.seed(2022) # 실험 재현성을 위해서 seed 를 설정한다.
# 데이터 명시
data1_jags = list(y=dat$orders, n=nrow(dat),
temperature=dat$temperature)
# 파라미터 명시
params1 = c("b", "sig")
# 초기함수 명시
inits1 = function() {
inits = list("b" = rnorm(2, 0.0, 10.0),
"prec" = rgamma(1, 1.0, 1.0))
}
# 모델 명시, 체인은 3개 구동
mod = jags.model(textConnection(mode1_string), data=data1_jags, inits=inits1, n.chains=3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 14
## Unobserved stochastic nodes: 3
## Total graph size: 64
##
## Initializing modelupdate(mod, 5000) # burn-in 실행
# 샘플링 수행
mod_sim = coda.samples(model=mod, variable.names=params1, n.iter=20000)
mod_csim = do.call(rbind, mod_sim) # 체인을 행 단위로 결합한다.Model checking
시뮬레이션 모델에 대한 traceplot 을 그려본다. 여기서 \(\beta_0\), \(\beta_1\) 에 대한 플롯을 보면, 3개의 체인은 장기적 추세를 파악하기 어렵고 일정한 값으로 수렴하고 있음을 볼 수 있다.
plot(mod_sim)
Autocorrelation 을 보자. 절편 항과 기울기 항에서 초기 지연에 대해 높은 자기 상관을 보이고 있다. 따라서 번인으로 최소 500 번 이상 수행해야만 함을 알 수 있다.
autocorr.diag(mod_sim)## b[1] b[2] sig
## Lag 0 1.0000000 1.0000000 1.00000000
## Lag 1 0.9896460 0.9895272 0.10997601
## Lag 5 0.9493934 0.9492563 0.04912665
## Lag 10 0.9012823 0.9010427 0.04457158
## Lag 50 0.5871947 0.5867184 0.01521664autocorr.plot(mod_sim, lag.max=1000)
Gelman and Rubin’s Diagostics 를 검토해보자. 3개의 서로 다른 파라미터 \(\beta_0\), \(\beta_1\), \(\sigma\) 에 대해서 Potential Scale Reduction Factors 가 1000번 이후로 상당히 1에 가깝다는 것을 볼 수 있다. 그것은 우리가 수렴했음을 알려준다.
gelman.diag(mod_sim) # 3개의 서로 다른 파라미터에 대한 잠재 스케일 감소 인자를 본다. ## Potential scale reduction factors:
##
## Point est. Upper C.I.
## b[1] 1.02 1.07
## b[2] 1.02 1.07
## sig 1.00 1.00
##
## Multivariate psrf
##
## 1.02gelman.plot(mod_sim)
효과표본크기를 보자. 예측을 위해 사후분포 평균을 이용한다면, 효과표본크기가 작아도 괜찮다.
effectiveSize(mod_sim) # 사후분포 평균이라면, 효과표본크기가 작아도 okay.## b[1] b[2] sig
## 312.2282 314.8890 11673.3202하지만, 95% 확률 구간과 같은 정보를 구하기 위해서는 얼마만큼의 효과표본크기가 필요한지는 아래처럼 알아보자.
Results
모델에 대한 사후분포 요약을 보자. 우리는 처음 5000 개의 번인을 사용했다. 그리고 20000 개를 표본으로 저장했다. 3개의 체인을 실행했고 여기에 파라미터에 대한 결과 추정치가 있다. 사후분포 평균에 기반해서 회귀식의 절편 b[1]은 약 -37.7 이 되고, 계수 b[2] 는 약 3.8 정도 된다.
summary(mod_sim)##
## Iterations = 5001:25000
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 20000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## b[1] -35.568 14.4207 0.058872 0.816104
## b[2] 3.711 0.4920 0.002009 0.027699
## sig 5.344 0.9772 0.003990 0.009204
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## b[1] -63.327 -45.525 -35.555 -25.970 -6.586
## b[2] 2.722 3.382 3.711 4.052 4.659
## sig 3.834 4.654 5.208 5.885 7.621참조 베이지안 분석(Reference Bayesian Analysis)을 해보자. \(\beta_0\) 의 계수는 -36.3 으로 참조 모델과 베이지안 모델이 비슷하며, \(\beta_1\) 의 계수도 3.7 로 비슷함을 알 수 있다. 그리고 잔차 표준오차(Residual standard error) 가 약 5.7이고 \(R^2\) 이 약 0.82 로 나타났다.
mod_lm <- lm(orders ~ temperature, data=dat)
summary(mod_lm)##
## Call:
## lm(formula = orders ~ temperature, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.037 -5.693 2.094 4.409 8.225
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.3612 14.6873 -2.476 0.0292 *
## temperature 3.7379 0.5012 7.457 7.66e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.709 on 12 degrees of freedom
## Multiple R-squared: 0.8225, Adjusted R-squared: 0.8077
## F-statistic: 55.61 on 1 and 12 DF, p-value: 7.661e-06# 잔차에 있어서 특별한 패턴이 보이지 않는다.
#plot(resid(mod_lm))
#plot(predict(mod_lm), resid(mod_lm))사후분포 표본을 이용한 잔차를 비교해보자. 잔차 플롯에서 특별한 패턴이 보이지 않는다.
X <- cbind(rep(1.0, data1_jags$n), data1_jags$temperature)
pm_params1 = colMeans(mod_csim) # posterior mean: 각 컬럼(파라미터에 대한 평균값)
pm_params1## b[1] b[2] sig
## -35.567616 3.710871 5.343723# 파라미터에 대한 사후분포 평균과 관측치를 이용해서 예측된 값을 계산한다.
yhat1 = drop(X %*% pm_params1[1:2])
# 잔차 플롯 그리기
resid1 = data1_jags$y - yhat1
plot(resid1) # against data index
plot(yhat1, resid1)
#sd(resid1) # 잔차에 대한 표준편차 확인Conclusions
모델 검토를 통해서, 참조 모델과 베이지안 선형 회귀 모델이 거의 비슷한 파라미터 값을 추정함을 알았다. 또한, 잔차를 검토해서 데이터의 비선형성을 확인하지 못했다.
계수에 대한 사후분포 추정치를 사용해서 최고 기온이 27 도 일 때 아이스 커피의 주문 수는 어떻게 될까? 약 65잔의 주문 수를 예측할 수 있다.
pm_params1['b[1]'] + pm_params1['b[2]'] * 27## b[1]
## 64.6259사후분포를 이용해서 시뮬레이션 해보자. 약 95% 확률로 62 ~ 68 잔이 주문 수량으로 예측된다.
yhat2 <- mod_csim[,'b[1]'] + mod_csim[,'b[2]'] * 27
hist(yhat2, breaks = 100)
quantile(yhat2, c(0.025, 0.975))## 2.5% 97.5%
## 61.09826 68.22106최고 기온이 27도일 때, 70잔 이상 주문될 확률은 얼마일까? 0.2% 정도가 됨을 알 수 있다.
mean (yhat2 > 70)## [1] 0.002633333References
[1] Takahashi Shin, 2011, 만화로 쉽게 배우는 회귀분석, 성안당