Bayesian P- Value
- Bayesian P-Value는 **모델의 적합도(goodness-of-fit)**를 평가하는 데 사용됩니다.
- 사후 분포(posterior distribution)를 이용하여 실제 데이터와 모델이 생성한 예상 데이터를 비교함으로써, 관측된 데이터가 모델에 의해 얼마나 잘 설명되는지를 평가합니다.
- 빈도주의 p-값은 "관찰된 데이터보다 극단적인 데이터가 나올 확률"을 계산하지만, Bayesian P-Value는 "모델이 실제 데이터를 얼마나 잘 설명하는지"를 평가하는 데 중점을 둡니다.
계산 방법
- Bayesian P-Value는 주로 두 가지 값을 비교합니다:
- 관측된 통계량(Observed Test Statistic): 실제 데이터에서 계산한 통계량입니다
- 예측된 통계량(Predicted Test Statistic): 사후 분포를 기반으로 모델이 예측한 데이터에서 계산한 통계량입니다.
- 이 두 통계량을 비교하여 모델이 실제 데이터를 얼마나 잘 설명하는지 평가할 수 있습니다. $$P_B = P(T(\tilde{y}) \geq T(y) \mid y)$$
- $T(\tilde{y})$: 모델에 의해 예측된 데이터 $\tilde{y}$에서 계산된 통계량
- $T(y)$: 실제 관측 데이터 $y$에서 계산된 통계량
- 베이지안 p-값은 예측된 데이터의 통계량이 실제 관측 데이터보다 크거나 같은 확률로 정의됩니다.
장단점
장점
- 베이지안 방식의 장점 반영: 불확실성을 자연스럽게 모델링하고, 사전 정보(prior knowledge)를 반영할 수 있습니다.
- 적합도 평가: 모델이 데이터를 얼마나 잘 설명하는지 평가하는 유연한 도구로 활용될 수 있습니다.
단점
- 직관성 부족: 빈도주의적 p-값처럼 단순히 임계값(예: 0.05)을 넘으면 가설을 기각하는 명확한 기준이 없습니다. 대신, 적합도를 상대적으로 평가해야 합니다.
- 계산 복잡성: 사후 분포를 샘플링하고, 이를 바탕으로 통계량을 계산하는 과정이 빈도주의적 방법에 비해 더 복잡하고 시간이 소요될 수 있습니다.
코드
- P-value가 0.75라는 것은, 모델이 예측한 값들이 실제 관측 데이터와 매우 유사하거나 더 나은 설명을 제공할 수 있는 확률이 75%라는 의미입니다. 즉, 모델이 우리가 관측한 데이터를 비교적 잘 설명하고 있다는 신호로 볼 수 있습니다.
import pymc as pm
import numpy as np
import arviz as az
# 가상의 데이터 생성
true_mu = 5
observed_data = np.random.normal(true_mu, 1, size=100)
# 모델 정의
with pm.Model() as model:
# 사전 분포 정의
mu = pm.Normal('mu', mu=0, sigma=10)
sigma = pm.HalfNormal('sigma', sigma=1)
# 가능도 정의
y = pm.Normal('y', mu=mu, sigma=sigma, observed=observed_data)
# 사후 분포에서 샘플링
trace = pm.sample(1000)
# 사후 예측 분포에서 샘플링
posterior_predictive = pm.sample_posterior_predictive(trace, var_names=['y'])
# Bayesian P-value 계산
discrepancy = lambda x: np.mean(x)
T_obs = discrepancy(observed_data)
# posterior_predictive에서 'y' 추출
pp_y = posterior_predictive.posterior_predictive['y'].values
# T_rep 계산
T_rep = np.array([discrepancy(y) for y in pp_y])
p_value = np.mean(T_rep >= T_obs)
print(f"Bayesian P-value: {p_value}")
# 추가 정보 출력
print(f"\nPosterior predictive shape: {pp_y.shape}")
print(f"T_rep shape: {T_rep.shape}")
print(f"Observed data shape: {observed_data.shape}")
