NumPyro를 통한 베셀 선형 단일 회귀

개시하다

나는 확률 모델링에 푹 빠졌다.최근 페이로에서 MC를 수행하던 중 느린 것을 참지 못해 넘페이로에 입성했다.따라서 이 글에서 학습하는 동시에 NumPyro에서 베스탄 선형 단일 회귀를 실시한다.또한 가시화하려면 베이스 모델링의 가시화 라이브러리arviz를 사용하십시오.
소스 코드는 GiitHub에 공개됩니다.
https://github.com/kajyuuen/playground

장난감 데이터 만들기

이번에 사용한 데이터는 위조 데이터를 사용할 것이다.
구체적으로 우리는 다음과 같은 공식이 표현할 수 있는 데이터를 고려할 것이다.
y_i = 0.9x_i + 5 +\epsilon_{\rm real}
\epsilon_{\rm real}\sim {\mathcal N}(0, 0.5)
그럼 데이터 제작을 시작합시다.
우선 필요한 라이브러리 그룹을 준비합니다.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import arviz as az
import numpyro
import numpyro.distributions as dist
import jax
plt.style.use('seaborn-darkgrid')

다음에 Numpy를 사용하여 데이터를 생성합니다.

N = 100
alpha_real, beta_real = 0.9, 5
epsilon_real = np.random.normal(0, 0.5, N)

x = np.random.normal(10, 1, N)
y_real = alpha_real * x + beta_real
y = y_real + epsilon_real

생성 데이터를 그린 결과는 다음과 같다.

베일스 선형 단회귀 모델

생성 과정

그럼 베스탄의 선형 단회귀 모델을 고려해 보자.
선형 단일 회귀는 기울기\alpha이고 슬라이스가\beta일 때 변수 yi에 관해서는 다음과 같이 쓸 수 있다.
y_i =\alpha x_i +\beta
이 선형 단일 회귀의 확률을 사용하여 다음과 같이 나타낸다.
{\boldsymbol y}\sim {\mathcal N}(\mu =\alpha {\boldsymbol x} +\beta,\sigma =\epsilon)
즉, 데이터 서열 {\boldsymboly}이 평균값\alpha{\boldsymbolx]+\beta, 표준 편차\epsilon의 정적 분포를 따른다고 가정한다.그러나 각 매개 변수\alpha,\beta,\epsilon에 관해서는 실제 값을 모르기 때문에 다음과 같이 미리 분포한다고 가정한다.이 예비 분포는 적당한 분포라면 다 된다.
\alpha\sim {\mathcal N}(\mu_\alpha,\sigma_\alpha)
\beta\sim {\mathcal N}(\mu_\beta,\sigma_\beta)
\epsilon\sim {\rm HalfCauchy(\sigma_\epsilon)}
이러한 사전 분포와 확률 모델의 도형 모델은 다음과 같이 쓸 수 있다.

이것은 NumPyro로 구현됩니다.

def model(x, y = None):
    alpha = numpyro.sample("alpha", dist.Normal(0, 1))
    beta = numpyro.sample("beta", dist.Normal(0, 10))
    eplison = numpyro.sample("eplison", dist.HalfCauchy(5))
    
    mu = numpyro.deterministic("mu", alpha * x + beta)
    
    with numpyro.plate("data", size = len(x)):
        y_pred = numpyro.sample("y_pred", dist.Normal(mu, eplison), obs=y)

추론

이어 MCMC에 따라 샘플링을 진행한다.이번에는 NO-U-turn sampler(NUTS)를 사용했다.

kernel = numpyro.infer.NUTS(model)
mcmc = numpyro.infer.MCMC(kernel, num_samples=1000, num_warmup=300, num_chains=4, chain_method="parallel")
mcmc.run(jax.random.PRNGKey(0), x, y)

arviz를 사용하여 샘플링 결과를 가시화합니다.

az.plot_trace(mcmc, var_names=["alpha", "beta", "eplison"])

왼쪽은 내핵 밀도를 추정하는 도표다.나는 중심극제한리를 통해 정규 분포와 비슷한 도표임을 확인할 수 있다고 생각한다.
이것은 오른쪽에 쓴 견본에서 얻은 값이다.규칙성과 주기성이 보이지 않는 도표로 잘 섞인 모습을 엿볼 수 있다.
또 print_summary()를 사용하면 주요 통계량을 볼 수 있다.

mcmc.print_summary()

                mean       std    median      5.0%     95.0%     n_eff     r_hat
     alpha      0.95      0.05      0.95      0.86      1.04   1090.27      1.01
      beta      4.56      0.54      4.56      3.71      5.44   1089.63      1.01
   eplison      0.51      0.04      0.51      0.45      0.57   1546.22      1.00

Number of divergences: 0

예측 분포
다음은 후방 분포의 샘플링에 따라 예측 분포를 생성한다.
먼저 get_samples()를 사용하고 후분포에서 각 매개 변수\alpha,\beta,\epsilon을 샘플링한다.

posterior_samples = mcmc.get_samples()

다음에 샘플을 채취한 매개 변수\alpha,\beta,\epsilon을 사용하여 {boldsymboly}의 예측 분포를 구한다.

posterior_predictive = numpyro.infer.Predictive(model, posterior_samples)(
    jax.random.PRNGKey(1), x
)

실제 데이터를 그리는 내핵밀도 추정과 사후분포 계산의 내핵밀도 추정.

_, ax = plt.subplots()
for y_pred in posterior_predictive["y_pred"][0:50]:
    sns.kdeplot(y_pred, alpha=0.1, c='r', ax=ax)
# real data
sns.kdeplot(y, linewidth=3, color='k', ax=ax)
    
plt.xlabel('$y$', fontsize=16);

그림의 도표가 실제 데이터와 매우 가깝다는 것을 확인할 수 있다.
마지막으로 샘플링된 직선과\pm\sigma,\pm3\sigma의 불확실성을 테이프로 표현해 봤다.

def summary(samples):
    s_mean, s_std = samples.mean(0), samples.std(0)
    site_stats = {
        "mean": s_mean,
        "std": s_std,
        "upper_sigma1": s_mean + s_std,
        "lower_sigma1": s_mean - s_std,
        "upper_sigma3": s_mean + 3 * s_std,
        "lower_sigma3": s_mean - 3 * s_std
    }
    return site_stats
pred_summary = summary(posterior_predictive["y_pred"])

plt.plot(x, y, 'C0.');

# 直線
alpha_m = posterior_samples['alpha'].mean()
beta_m = posterior_samples['beta'].mean()
plt.plot(x, alpha_m * x + beta_m,
         c='k', label='y = {:.2f} * x + {:.2f}'.format(alpha_m, beta_m))

# 範囲
idx = np.argsort(x)
x_ord = x[idx]
plt.fill_between(x_ord, pred_summary["upper_sigma1"][idx], pred_summary["lower_sigma1"][idx], alpha=0.3, color="b")
plt.fill_between(x_ord, pred_summary["upper_sigma3"][idx], pred_summary["lower_sigma3"][idx], alpha=0.1, color="b")



plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16, rotation=0)
plt.legend(loc=2, fontsize=14)

좋은 느낌이 나오지 않나요?
넘페로는 파이로에 비해 MCMC의 속도가 빠르다. 앞으로 MCMC를 쓰면 넘페로, 변분 베이즈법, 신경망 쪽에서 여러 가지 신경 쓰는 일을 하고 싶다면 각각 파이로를 사용할 수 있지 않을까... 하는 생각이 든다. 어떨까?자세한 사람 있으면 알려주세요.

참고 자료

파이톤 기반 베이스 통계 모델링
https://www.kyoritsu-pub.co.jp/bookdetail/9784320113374