confidence_interval_demo
Sang woo Ham
October 10, 2017
Confidence interval 신뢰 구간
예시를 들어서 설명하자
모집단: 데이터 \(x_i\)의 모분포를 알고 있다고 가정하자. \(x_i\sim \text{Gamma}(\alpha,\beta)=\text{Gamma}(6,1)\).
샘플: 이 모집단으로 부터 20개의 샘플 데이터를 얻었다고 하자 \(x_i^{[1]}=\{x_1^{[1]},x_2^{[1]},...x_{20}^{[1]} \}\). 이때, \(x_i^{[k]}\)은 \(k\) 번째 샘플의 i 번째 데이터를 의미한다.
Central limit theorem (CLT) 중심 극한 정리
Central limit theorem에 따라서 샘플의 평균은 다음과 같은 분포를 따른다.
\[\bar{\textbf{x}}^{[k]}\sim \mathcal{N}\bigg(\mu_p,\frac{\sigma_p^2}{n^{[k]}}\bigg)\]
Confidence interval of sample mean \(\bar{\textbf{x}}\) (샘플평균의 신뢰구간)
Sample mean의 confidence interval은 다음과 같이 표현된다. 이때, 모분포의 표준편차 \(\sigma_p\)를 모를 경우, t 분포와 sample mean의 standard deviation (standard error, \(s^{[k]}\))를 이용해서 근사한다.
\[\bar{\textbf{x}}^{[k]}\pm z_{\alpha/2}\sqrt{\frac{\sigma^2_p}{n^{[k]}}}\approx\bar{\textbf{x}}^{[k]}\pm t_{\alpha/2}\sqrt{\frac{s^{[k]}}{n^{[k]}}}\] 여기서 \(s^{[k]}\)는 estimated standard deviation (standard error). 또한, 모집단의 \(\sigma_p\)를 모르기 때문에, 우리는 샘플의 평균이 t분포를 따른다고 가정한다 (central limit theorem) \(t(\bar{\textbf{x}}^{[k]},\,s^{[k]}/\sqrt{n})\) with \(n^{[k]}-1\) degree of freedom. 여기서 \(s^{[k]}=\frac{1}{n^{[k]}-1}\sum_{i=1}^{n^{[k]}} \big(x_i^{[k]}-\bar{\textbf{x}}^{[k]}\big)^2\) (see this link).
R로 계산을해보자
우선 모집단 분포를 생성하자. \(\text{Gamma}(6,1)\) 감마 분포의 평균과 표준편차는 다음과 같이 구한다 \(\mu_p=\alpha/\beta=\)shape/rate and \(\sigma_p=\alpha/\beta^2\) shape/rate\(^2\)
shape = 6
rate = 1
mu_p = shape/rate # population mean
sd_p = sqrt(shape/(rate^2)) # poplulation std
set.seed(123) # seed
pdist <- rgamma(1e+06, shape, rate) # generate population data50개의 샘플(num_sample)을 관측하는데 이 과정을 1000번(num_repeat) 반복한다고 하자.
그리고 num_repeat x \(3\) 매트릭스를 만들어서 각 행에는 반복 시행(\([k]\))을 첫번째 열에는 신뢰구간의 lower bound \(\bar{\textbf{x}}^{[k]}- t_{\alpha/2}\sqrt{\frac{s^{[k]}}{n^{[k]}}}\)을 두번째 열에는 샘플의 평균 \(\text{mean}(\bar{\textbf{x}}^{[k]})\), 그리고 세번째 열에는 신뢰구간의 upper bound \(\bar{\textbf{x}}^{[k]}+ t_{\alpha/2}\sqrt{\frac{s^{[k]}}{n^{[k]}}}\)를 저장하자.
num_repeat = 10000 # number of sampling
num_sample = 20 # number of samples in each sample
output <- data.frame(matrix(0, ncol = 3, nrow = num_repeat)) # 1st col: lower bound, 2nd: mean, 3rd: upper
alpha = 0.95
for (i in 1:num_repeat) {
set.seed(i * 40 + 200)
vec <- sample(pdist, num_sample, replace = T) # sample from population
output[i, 2] <- mean(vec) # sample mean
ssd <- sd(vec) # sample stdanrd deviation
output[i, 3] <- mean(vec) + t(qt(alpha, df = num_sample - 1) * ssd/sqrt(num_sample)) # alpha=95% (0.025-0.975)
output[i, 1] <- mean(vec) - t(qt(alpha, df = num_sample - 1) * ssd/sqrt(num_sample))
}우선 CLT를 그려보자. 이때, \(\alpha=\mu_p^2/(\sigma_p^2/n_{\text{sample}})\) and \(\beta=\mu_p/(\sigma_p^2/n_{\text{sample}})\)
hist(output[, 2], breaks = 100, freq = F, main = expression("Red line: Gamma" *
(mu[p] * "," * sigma[p]/sqrt(n)))) # histogram of sample mean
points(seq(3, 9, 0.01), dgamma(seq(3, 9, 0.01), shape = mu_p^2/(sd_p^2/num_sample),
rate = mu_p/(sd_p^2/num_sample)), type = "l", col = "red")
즉, 감마 분포에서 추출한 샘플의 평균이 정규 분포를 잘 따르는 것을 알 수 있다.
다음은 신뢰구간을 알아 보자.
신뢰구간의 기본 의미를 보면 (90% 기준), 각 샘플의 신뢰구간에 모평균이 전체 반복 횟수 중 90%가 포함 되어야 한다.
아래 그래프를 보면 not include 라고 표시된 4개 (10%)가 모평균을 포함하지 않는 것을 확인할 수 있다.
n_plot = 40
iixx <- which((output[1:n_plot, 1] > mu_p) | (output[1:n_plot, 3] < mu_p)) # index for plot
plot(output[1:n_plot, 2], seq(n_plot), xlim = c((mu_p - 1.5 * (sd_p)), (mu_p +
1.5 * (sd_p))), xlab = expression(x), ylab = "Replication index")
arrows(output[1:n_plot, 1], seq(n_plot), output[1:n_plot, 3], seq(n_plot), length = 0.02,
angle = 0, code = 3)
abline(v = mu_p, col = "red") #True mean
for (i in 1:length(iixx)) {
text(min(output[iixx[i], ]), iixx[i] - 0.2, "Not include", col = "blue")
}
그리고 전체 반복에 대해서 결과를 구해보자. 즉, (각 샘플의 신뢰구간이 모평균을 포함한 총 횟수)/(전체 반복횟수)=\(\alpha\)=0.9가 나와야 한다.
print(sum((output[, 1] <= mu_p) & (output[, 3] >= mu_p))/num_repeat)## [1] 0.89990.9에 근사한 값을 얻을 수 있다.
이것이 가지는 의미는 무엇일까? 우리가 모집단의 분포를 모른다고 할 때, 우리는 샘플 데이터만을 수집해서 분석할 수 있으며, CLT에 따라서 샘플의 신뢰구간을 구하면 그 작업을 여러번 반복 하였을 때 90% 정도는 모집단의 평균이 그 안에 들어 온다는 것이다. 즉, 모집단의 평균이 포함될 수도 있는 구간을 대략적으로나마 알 수 있다는 것이다.
하지만, 주의해야할 점이 있다. 예를 들어서 우리가 구하고자 하는 것이 어떤 사건이 일어날 횟수라고 하자. 샘플로 부터 구한 신뢰구간이 \(3\pm2\) 라고 한다면, “우리는 주로 이 사건이 1-5회정도 90% 확률로 일어난다” 라고 말하고 싶지만, 사실을 그렇지 않다. “여러개 신뢰구간의 90%는 모평균을 포함한다” 라고 말하는 것이 더 맞을 것이다. (아마도 그래서 “모평균이 신뢰구간에 들어올 것은 90% 정도 자신한다” 라는 표현을 쓰는 것 같다).
말이 길지만 중요한 것은 하나다. 신뢰구간에는 모평균이 들어올 수도 있고 안들어 올 수도 있다는 것. 샘플 과정을 반복하면 모평균이 들어올 것을 \(\alpha\)% 확신할 수 있다는 것.
모집단의 신뢰구간?
그렇다면, 표현이 정확하진 않지만, 모집단의 신뢰구간에는 샘플의 평균이 들어갈 것인가? 사실 정확한 식을 찾지 못해서 유추하여 작성하였다.
앞서 구한 CLT의 분포를 생각해본다면, \(\bar{\textbf{x}}^{[k]}\sim \mathcal{N}\bigg(\mu_p,\frac{\sigma_p^2}{n}\bigg)\)
\[\mu_p\pm z_{\alpha/2}\sqrt{\frac{\sigma^2_p}{n^{[k]}}}\approx\mu_p\pm t_{\alpha/2}\sqrt{\frac{s^{[k]}}{n^{[k]}}}\]
사실 빈도주의 통계학에서는 모평균, 모분포라 개념을 구하기가 어렵기 때문에 이러한 구간을 설정하는 예시를 많이 보지 못했다. 그도 그럴것이 모평균을 알면 누가 샘플 조사를 하겠는가? (물론 샘플 조사의 타당성 같은 것을 증명하는 경우를 제외하곤 말이다). 여하튼 범위를 계산하여 구해보자.
\(\mu_p\pm z_{\alpha/2}\sqrt{\frac{\sigma^2_p}{n^{[k]}}}\)의 경우에는 약 0.9가 나타나는 것을 알 수 있다.
mu_p_lower = mu_p - sd_p * qnorm(alpha, 0, 1)/sqrt(num_sample)
mu_p_upper = mu_p + sd_p * qnorm(alpha, 0, 1)/sqrt(num_sample)
print(sum((output[, 2] >= mu_p_lower) & (output[, 2] <= mu_p_upper))/num_repeat)## [1] 0.9047\(\mu_p\pm t_{\alpha/2}\sqrt{\frac{s^{[k]}}{n^{[k]}}}\)의 경우에도 약 0.9가 나타나는 것을 알 수 있다.
mu_p_s_lower = mu_p - (output[, 2] - output[, 1])
mu_p_s_upper = mu_p + (output[, 3] - output[, 2])
print(sum((output[, 2] >= mu_p_s_lower) & (output[, 2] <= mu_p_s_upper))/num_repeat)## [1] 0.8999여러개의 샘플 평균이 모평균의 신뢰구간에 90% 정도는 존재한다 라고 결론을 내려야 할까? 그렇게 해석할 수도 있겠지만, 앞서 말했듯이, 우리가 구할 수 있는 것은 샘플과 샘플의 신뢰구간이다. 이 정보를 가지고 모평균의 신뢰구간을 사실 우리는 알 수 없다. 따라서 그냥 앞서 내린 것 과 같이 “모평균이 신뢰구간에 존재한다고 90% 확신한다” 와 같은 결론을 내릴 수 밖에 없는 것이다. 사실 위의 그래프를 생각해본다면, 이는 각각 샘플의 신뢰구간만을 잘라다가 모평균에 붙혀놓은 것과 다를바 없다.
n_plot = 40
plot(output[1:n_plot, 2], seq(n_plot), xlim = c((mu_p - 1.5 * (sd_p)), (mu_p +
1.5 * (sd_p))), xlab = expression(x), ylab = "Replication index")
arrows(mu_p - output[1:n_plot, 2] + output[1:n_plot, 1], seq(n_plot), mu_p -
output[1:n_plot, 2] + output[1:n_plot, 3], seq(n_plot), length = 0.02, angle = 0,
code = 3, col = "red")
abline(v = mu_p, col = "red") #True mean
마지막으로 어떤 사건이 1-5회 일어날 확률이 90%다 같은 결론을 내릴 수는 없는 것일까? 어차피 모평균을 추정하는 방식의 하나이긴 하지만, 베이지안 기법을 사용하면 된다.
이를 Credible interval이라 한다.
Credible Interval
정리가 힘들어서 수식관련해서는 예전 정리 자료를 복사 붙혀넣기 하였다.
Let’s say our population is \(x_i^{[k]}\sim \text{Gamma}(\alpha,\beta)\)
The mean (\(\mu_p\)) and standard deviation\(\sigma_p\) of population distribution is \(\mu_p=\alpha/\beta=\)shape/rate and \(\sigma_p=\alpha/\beta^2\) shape/rate\(^2\)
Now come back to Bayesian. That is actually the starting point. I want to discuss the credible interval actually.
Here, I am not going to cover all the basics of Bayes theorem.
From Arthur Charpentier’s site, parameters are random and data are fixed. Given the data \(\textbf{x}\), we infer the posterior’s distribution (\(\text{Pr}(\mu_p,\sigma_p|\textbf{x})\)).
And then, we finally get posterior predictive distribution \[\text{Pr}(\tilde{\textbf{x}}|\textbf{x})=\int_{\theta}\text{Pr}(\tilde{\textbf{x}}|\textbf{x},\theta)\text{Pr}(\theta|\textbf{x})\] This might be confusing, but we can simply sample posterior \(\text{Pr}(\theta|\textbf{x})\) and then replicate the \(\tilde{\textbf{x}}\) from the given the posterior.
From the previous analysis, the likelihood is \[\prod_i^{n^{[k]}}\text{Pr}(x^{[k]}_i|\alpha,\beta)\]
Let’s calculate MLE solution. Let’s assume we have 50 samples. We will use same shape and rate (i.e., \(\alpha\) and \(\beta\)).
set.seed(123)
x = rgamma(50, shape = shape, rate = rate)The mle solution will maximize the log-likelihood.
alpha_grid <- seq(1, 10, 0.01)
beta_grid <- seq(1, 10, 0.01)
mat_df <- matrix(0, nrow = length(alpha_grid), ncol = length(beta_grid))
for (i in 1:length(alpha_grid)) {
for (j in 1:length(beta_grid)) {
mat_df[i, j] = sum(dgamma(x, shape = alpha_grid[i], rate = beta_grid[j],
log = T))
}
}
iiix <- which(mat_df == max(mat_df), arr.ind = T)
alpha_mle = alpha_grid[iiix[1]]
beta_mle = beta_grid[iiix[2]]
print(paste0("alpha_mle: ", alpha_mle, " and beta_mle: ", beta_mle))## [1] "alpha_mle: 6.43 and beta_mle: 1.2"This is close to population parameter \(\alpha=6\) and \(\beta=1\). Actually, if you have lots of data, MLE solution will converge to the population parameter.
In Bayesian approach, we define the posterior distribution
\[\text{Pr}(\theta|x)\propto\text{Pr}(x|\theta)\text{Pr}(\theta)\] In other words, \[\text{Pr}(\alpha,\beta|\textbf{x})\propto\text{Pr}(\textbf{x}^{[k]})|\alpha,\beta)\text{Pr}(\alpha,\beta)\]
Let’s do Bayesian inference. We assume \(\alpha\) and \(\beta\) are indepedent. \(\text{Pr}(\alpha,\beta)=\text{Pr}(\alpha)\text{Pr}(\beta)\).
Let’s assign priors on them. \(\alpha\sim\text{Unif}(4,8)\) and \(\beta\sim\text{Unif}(0.5,2)\). This is equivalent to \(\mu_p=[2,16]\) and \(\sigma_p=[1,32]\) I will use Stan for Bayesian inference.
library(rstan)
stan_code <- "
data{
int<lower=0> N; //number of sample observation
real<lower=0> x[N];
}
parameters{
real<lower=0> alpha;
real<lower=0> beta;
}
model{
alpha ~ uniform(4,8);
beta ~ uniform(0.5,2);
x ~ gamma(alpha,beta);
}
"
dat <- list(N = length(x), x = x)
bayes_out = stan(model_code = stan_code, data = dat, iter = 2000, chains = 1,
refresh = 1000)## In file included from C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from filed2c6ecb5904.cpp:8:
## C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## In file included from C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/multi_array/base.hpp:28:0,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/multi_array.hpp:21,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/numeric/odeint.hpp:61,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:13,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/arr.hpp:38,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/mat.hpp:298,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:12,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from filed2c6ecb5904.cpp:8:
## C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<N>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index_range index_range;
## ^
## C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index index;
## ^
## C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<0ull>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index_range index_range;
## ^
## C:/Users/SangWooHam/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index index;
## ^
## In file included from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:42:0,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from filed2c6ecb5904.cpp:8:
## C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:
## C:/Users/SangWooHam/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^
##
## SAMPLING FOR MODEL '0b019b52ac31ab4919fdf476f8e0d7c0' NOW (CHAIN 1).
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 0.059 seconds (Warm-up)
## 0.045 seconds (Sampling)
## 0.104 seconds (Total)summary(bayes_out)## $summary
## mean se_mean sd 2.5% 25%
## alpha 6.388872 0.09196131 0.9463190 4.6228479 5.627258
## beta 1.195294 0.01818227 0.1902520 0.8345222 1.047120
## lp__ -104.511720 0.05656687 0.8475645 -106.7398497 -104.866248
## 50% 75% 97.5% n_eff Rhat
## alpha 6.437943 7.182256 7.919058 105.8924 1.006255
## beta 1.208613 1.335963 1.534708 109.4869 1.003320
## lp__ -104.242075 -103.889396 -103.682307 224.5025 1.006150
##
## $c_summary
## , , chains = chain:1
##
## stats
## parameter mean sd 2.5% 25% 50%
## alpha 6.388872 0.9463190 4.6228479 5.627258 6.437943
## beta 1.195294 0.1902520 0.8345222 1.047120 1.208613
## lp__ -104.511720 0.8475645 -106.7398497 -104.866248 -104.242075
## stats
## parameter 75% 97.5%
## alpha 7.182256 7.919058
## beta 1.335963 1.534708
## lp__ -103.889396 -103.682307res <- rstan::extract(bayes_out)
alpha_map <- mean(res$alpha)
beta_map <- mean(res$beta)
print(paste0("alpha_map: ", alpha_map, " and beta_map: ", beta_map))## [1] "alpha_map: 6.38887239615815 and beta_map: 1.19529366261899"Thinking that \(\frac{\alpha}{\beta}=\mu_p\), the distribution of \(\sigma_p\) from the posterior is
hist(res$alpha/res$beta, breaks = 50, main = expression("Posterior distribution of " *
mu[p]), xlab = expression(mu[p]^s))
print("credible interval is")## [1] "credible interval is"print(quantile(res$alpha/res$beta, probs = c(0.05, 0.5, 0.95)))## 5% 50% 95%
## 4.909728 5.352753 5.845120# Confidence interval
lower_q <- mean(x) - t(qt(0.9, df = length(x) - 1) * sd(x)/sqrt(length(x)))
upper_q <- mean(x) + t(qt(0.9, df = length(x) - 1) * sd(x)/sqrt(length(x)))
round(lower_q, 2)## [,1]
## [1,] 4.97print(paste0("Confidence interval is [", round(lower_q, 2), ", ", round(mean(x),
2), ", ", round(upper_q, 2), "]"))## [1] "Confidence interval is [4.97, 5.36, 5.76]"The two value looks similar, but the thing is credible interval can be interpreted as “the population mean will be in the credible interval with 90% probability”.
Of course, we have parameter distributions, so we can generate lots of replicated samples. However there is no guarantee that this analysis is truely same with the real population distribution.