Chapter 2.3 - 예제

2.2.8. 예제: 우주 왕복선 챌린저호 참사¶

In [1]:

from IPython.core.pylabtools import figsize
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import scipy.stats as stats
import pymc3 as pm
import theano.tensor as tt

matplotlib.rc('font', family='Malgun Gothic') # 한글표시
plt.rc('axes', unicode_minus=False) # 마이너스 기호 표시

WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

In [40]:

figsize(8, 3.5)
np.set_printoptions(precision=3, suppress=True) # float 프린트 옵션 조절

challenger_data = np.genfromtxt('challenger_data.csv', skip_header=1, usecols=[1, 2], missing_values='NA', delimiter=',')

In [30]:

~np.isnan(challenger_data[:, 1])

Out[30]:

array([ True,  True,  True, False,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True, False])

In [31]:

# drop nan
challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])] 

In [32]:

plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color='k', alpha=.5)
plt.ylabel('사고: 1')
plt.xlabel('외부 기온(화씨)')
plt.show()

• 외부 온도가 낮을수록 피해 사고가 발생할 확률이 증가하는 경향이 보임. 기온과 발생할 사고 사이에 엄격한 구분이 보이지 않으므로 확률로 모델링.
• '온도 t에서 손실사고의 확률은 얼마일까?'에 대한 답을 확인 가능.

• $p(t)$라는 온도 함수가 필요함. 이 함숫값의 범위는 0과 1 사이이고, 온도를 높이면 1부터 0까지 변함. 이런 함수가 실제로 많고, 가장 대중적인 함수는 로지스틱 함수 $$p(t) = \frac 1 {1+e^{\beta t}}$$
• 이 모델에서 $\beta$는 확신이 없는 변수. $\beta=1,3,-5$에 대해 작성한 함수를 아래 그림으로 나타냄

In [6]:

figsize(12, 3)

def logistic(x, beta):
    return 1 / (1+np.exp(beta*x))

x = np.linspace(-4, 4, 100)
plt.plot(x, logistic(x, 1), label=r"$\beta = 1$")
plt.plot(x, logistic(x, 3), label=r"$\beta = 3$")
plt.plot(x, logistic(x, -5), label=r"$\beta = -5$")
plt.legend()

Out[6]:

<matplotlib.legend.Legend at 0x116dbad1760>

• 로지스틱 함수에서 확률은 0 근처에서 변하는데, 챌린저호 데이터에서는 화씨 65~70도 근처에서 변함. 로지스틱 함수의 bias를 추가해서 변형 가능. $$p(t)=\frac 1 {1+e^{\beta t + \alpha}}$$

In [7]:

def logistic(x, beta, alpha=0):
    return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha))

x = np.linspace(-4, 4, 100)

plt.plot(x, logistic(x, 1), label=r"$\beta = 1$", ls="--", lw=1)
plt.plot(x, logistic(x, 3), label=r"$\beta = 3$", ls="--", lw=1)
plt.plot(x, logistic(x, -5), label=r"$\beta = -5$", ls="--", lw=1)

plt.plot(x, logistic(x, 1, 1), label=r"$\beta = 1, \alpha = 1$",
         color="#348ABD")
plt.plot(x, logistic(x, 3, -2), label=r"$\beta = 3, \alpha = -2$",
         color="#A60628")
plt.plot(x, logistic(x, -5, 7), label=r"$\beta = -5, \alpha = 7$",
         color="#7A68A6")

plt.legend(loc="lower left");

• 상수항 $\alpha$를 추가하여 곡선을 왼쪽 또는 오른쪽으로 이동시킴(그래서 bias라고 함)
• pymc에서 이를 모델링. $\beta$, $\alpha$ 모수들은 양수이거나, 구간이 있거나, 상대적으로 클 필요는 없으므로 다음에 소개할 정규확률변수로 모델링하기 좋음

2.2.9. normal distribution¶

• normal random variable(정규확률변수)는 $X\sim N(\mu, 1/\tau)$로 나타내며 두 개의 모수(평균 $\mu$와 정밀도(precision)을 나타내는 $\tau$로 나타냄
• 정규분포에 친숙하다면 $\tau ^{-1}$ 대신 $\sigma ^2$으로 알겠지만, 둘은 사실 역수관계임
• $\tau$가 작을수록 variance가 높아지고, 반대의 경우 variance가 작아짐. $\tau$는 항상 양수

$$f(x | \mu, \tau) = \sqrt{\frac{\tau}{2\pi}} \exp\left( -\frac{\tau}{2} (x-\mu)^2 \right)$$

In [8]:

import scipy.stats as stats

nor = stats.norm
x = np.linspace(-8, 7, 150)
mu = (-2, 0, 3)
tau = (.7, 1, 2.8)
colors = ["#348ABD", "#A60628", "#7A68A6"]
parameters = zip(mu, tau, colors)

for _mu, _tau, _color in parameters:
    plt.plot(x, nor.pdf(x, _mu, scale=1./_tau),
             label="$\mu = %d,\;\\tau = %.1f$" % (_mu, _tau), color=_color)
    plt.fill_between(x, nor.pdf(x, _mu, scale=1./_tau), color=_color,
                     alpha=.33)

plt.legend(loc="upper right")
plt.xlabel("$x$")
plt.ylabel("$x$ 밀도함수")
plt.title("정규확률변수 3개에 대한 확률분포");

• 우주 왕복선 챌린저호에 대한 모델링 계속 진행

In [9]:

temperature = challenger_data[:, 0]
D = challenger_data[:, 1] # 검출되었는지, 아닌지

# value를 주목할 것.
with pm.Model() as model:
    beta = pm.Normal("beta", mu=0, tau=0.001, testval=0)
    alpha = pm.Normal("alpha", mu=0, tau=0.001, testval=0)
    p = pm.Deterministic("p", 1.0/(1. + tt.exp(beta*temperature + alpha)))

• 우리는 확률을 가지고 있지만, 확률을 어떻게 관찰된 데이터에 연결할 수 있을까?
• 이 때, 베르누이 확률변수를 활용하여, 모델을 아래와 같이 만들 수 있다

$$ \text{결함사고, $D_i$} \sim \text{Ber}( \;p(t_i)\; ), \;\; i=1..N$$

• $p(t)$는 우리의 로지스틱함수이고, $t_i$는 관측한 온도이다. 코드에서 beta와 alpha의 값을 0으로 설정한 것에 주목
• beta와 alpha가 너무 크면 $p$는 1 또는 0이 되기 때문임. 계수값(beta와 alpha)를 0으로 설정하여 p값을 합리적인 시작값으로 설정함
• 이는 결과에 영향을 주지 않고, 사전확률에 어떤 부가적인 정보를 포함한다는 의미도 아님

In [11]:

# 베르누이 확률변수를 통해 확률과 관측치를 연결
with model:
    observed = pm.Bernoulli('bernoulli_obs', p, observed=D)
    
    # 3장에서 설명
    start = pm.find_MAP()
    step = pm.Metropolis()
    trace = pm.sample(120000, step=step, start=start)
    burned_trace = trace[100000::2]

100.00% [26/26 00:00<00:00 logp = -19.024, ||grad|| = 9.9071]

<ipython-input-11-1c6e439a85bf>:8: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.
  trace = pm.sample(120000, step=step, start=start)
Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>Metropolis: [alpha]
>Metropolis: [beta]

100.00% [484000/484000 02:52<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 120_000 draw iterations (4_000 + 480_000 draws total) took 194 seconds.
The number of effective samples is smaller than 10% for some parameters.

In [20]:

alpha_samples = burned_trace['alpha'][:, None] # 1차원으로 만드는게 좋음
beta_samples = burned_trace['beta'][:, None]

figsize(6, 3)

# histogram of the samples:
plt.subplot(211)
plt.title(r"모델의 모수 $\alpha$와 $\beta$의 사후확률 분포")
plt.hist(beta_samples, histtype='stepfilled', bins=35, alpha=0.85, label=r"posterior of $\beta$", color="#7A68A6", density=True)
plt.legend()

plt.subplot(212)
plt.hist(alpha_samples, histtype='stepfilled', bins=35, alpha=0.85, label=r"posterior of $\alpha$", color="#A60628", density=True)
plt.legend()
plt.show();

• $\beta$의 모든 표본값은 0보다 큼. 사후확률이 0 주변에 집중되었다면 우리는 $\beta=0$을 의심할 것이며, 이는 온도가 결함확률에 영향을 주지않는다는 것을 암시함
• 마찬가지로 $\alpha$ 사후확률 값은 모두 음수이고, 0과는 거리가 멈. 이는 $\alpha$가 0보다 상당히 작다는 것을 의미.
• 분포가 넓은 것을 통해, 모수가 무엇인지는 확신할 수 없음.

• 온도의 특정 값에 대한 기대확률을 살펴봄. 이를 위해 사후확률분포에서 구한 모든 표본에 대해 평균을 내어 $p(t_i)$의 추정값을 얻음

In [23]:

t = np.linspace(temperature.min() - 5, temperature.max()+5, 50)[:, None]
p_t = logistic(t.T, beta_samples, alpha_samples)

mean_prob_t = p_t.mean(axis=0)

In [33]:

plt.plot(t, mean_prob_t, lw=3, label="결함에 대한 사후확률")
plt.plot(t, p_t[0, :], ls="--", label="사후 실현")
plt.plot(t, p_t[-2, :], ls="--", label="사후 실현")
plt.scatter(temperature, D, color="k", s=50, alpha=0.5)
plt.title("두 가지 실현값을 포함한 결함확률을 사후 기댓값")
plt.legend(loc="lower left")
plt.ylim(-0.1, 1.1)
plt.xlim(t.min(), t.max())
plt.ylabel("확률")
plt.xlabel("기온")
plt.show();

• 또한, 구현 가능한 실제 내부 시스템 두 가지를 나타냄. 두 개 확률은 같음.
• 파란선은 점선 20000개를 모두 평균 낼 때 발생

• 결함확률에 대해 어떤 온도에서 가장 확신이 없을지에 대해 알기 위해, 95% 신뢰구간을 설정함

In [41]:

from scipy.stats.mstats import mquantiles

# 신뢰구간을 위해 맨 아래부터 상위 2.5%를 벡터화함
qs = mquantiles(p_t, [0.025, 0.975], axis=0)
plt.fill_between(t[:, 0], *qs, alpha=0.7, color="#7A68A6")

plt.plot(t[:, 0], qs[0], label="95% CI", color="#7A68A6", alpha=0.7)
plt.plot(t, mean_prob_t, lw=1, ls="--", color="k", label="결함 확률에 대한 평균 사후확률")

plt.xlim(t.min(), t.max())
plt.ylim(-0.02, 1.02)
plt.legend(loc="lower left")
plt.scatter(temperature, D, color="k", s=50, alpha=0.5)
plt.xlabel("기온, $t$")

plt.ylabel("확률 추정치")
plt.title("Posterior probability estimates given temp. $t$")
plt.show();

• 65도에서는 25~80% 확률로 결함확률이 있다고 95% 확신할 수 있음.

• 일반적으로 온도가 60도에 가까워지면 CI가 빠르게 퍼지고, 70도를 지나면 다시 조여짐. 이는 다음에 진행하는 방법에 대한 통찰력을 얻을 수 있음.
• 해당 범위의 확률을 더 잘 추정하려면 60-65도 온도 주변에서 더 많은 테스트가 수행되어야 함.

2.2.10. 챌린저호 참사 당일¶

• 챌린저호 참사 당일 외부 온도는 화씨 31도였음. 이 온도에서 결함이 발생할 사후확률분포는 어떠한가?

In [42]:

prob_31 = logistic(31, beta_samples, alpha_samples)

plt.xlim(0.995, 1)
plt.hist(prob_31, bins=1000, density=True, histtype='stepfilled')
plt.title("$t = 31$일 때 결함확률의 사후확률분포")
plt.xlabel("O-ring에서 발생하는 결함확률")
plt.show();

2.3 우리 모델은 적절할까?¶

• $p(t)$에 대해 의도적으로 로지스틱 함수를 선택했고 특정 사전확률을 선택했고, 다른 함수 또는 다른 사전확률은 다른 결과를 줄 것이라고 생각할 수 있음.
• 우리의 모델이 데이터를 표현한다는 것을 어떻게 알 수 있을까? 이를 위해 모델이 우리의 관측과 얼마나 적합한지를 측정함

• 이를 측정하는 방법은 관측 데이터(즉, 고정된 확률변수)와 시뮬레이션하는 인위적인 데이터셋을 비교하는 것.
• 시뮬레이션된 데이터셋이 관측데이터셋과 통계적으로 유사해보이지 않는다면 우리의 모델은 관측 데이터를 정확히 나타내지 못함

• 사후확률분포에서 표본을 추출해야함. 다행히 베이지안 프레임워크가 표본추출을 쉽게 해결해줌.
• 우리는 새 stochastic 변수를 만들면 됨.

In [43]:

# 1만번 시뮬레이션 해보자
N = 10000
with pm.Model() as model:
    beta = pm.Normal('beta', mu=0, tau=0.001, testval=0)
    alpha = pm.Normal('alpha', mu=0, tau=0.001, testval=0)
    p = pm.Deterministic('p', 1.0/(1.+tt.exp(beta*temperature+alpha)))
    observed = pm.Bernoulli('bernoulli_obs', p, observed=D)
    
    simulated = pm.Bernoulli('bernoulli_sim', p, shape=p.tag.test_value.shape)
    step = pm.Metropolis(vars=[p])
    trace = pm.sample(N, step=step)

<ipython-input-43-8248ce71447d>:11: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.
  trace = pm.sample(N, step=step)
Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>CompoundStep
>>Metropolis: [beta]
>>Metropolis: [alpha]
>BinaryGibbsMetropolis: [bernoulli_sim]

100.00% [44000/44000 00:23<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 10_000 draw iterations (4_000 + 40_000 draws total) took 46 seconds.
The rhat statistic is larger than 1.4 for some parameters. The sampler did not converge.
The estimated number of effective samples is smaller than 200 for some parameters.

In [65]:

simulations = trace["bernoulli_sim"]
print(simulations.shape)

plt.title("Simulated dataset using posterior parameters")
figsize(12.5, 6)
for i in range(4):
    ax = plt.subplot(4, 1, i+1)
    plt.scatter(temperature, simulations[1000*i, :], color="k", s=50, alpha=0.6)
    
plt.ylabel('사고?')
plt.xlabel('외부 기온(화씨)')
plt.show()

(40000, 23)

• 기반 데이터가 다르므로 각 도표도 다름. 그러나 데이터셋은 모두 같은 기본모델에서 파생됨
• 무작위성으로 인해 데이터셋들은 다르게 보이며, 통계는 데이터셋을 하나로 묶어줌.

• 모델이 잘 만들어졌는지 그래픽적으로 작업함. 그래서 덜 객관적인 방법처럼 보일 수 있수 있음.
• 다른 방법은 통계적으로 요약한 베이지안-p 값을 사용하는 것. 빈도주의 통계에서의 p-value와 비슷하지만, 여전히 주관적임(좋고 나쁨 사이의 적절한 기준이)

2.3.1. 분리도표¶

• 다음의 그래픽 테스트는 로지스틱 회귀분석법을 위한 새로운 데이터 시각화 방법
• 이런 도표를 분리도표라고 함.
• 분리도표는 사용자가 비교하고 싶은 서로 다른 모델을 그래픽으로(그래프로) 비교.

• 각 모델에 대해 사후확률 시뮬레이션이 특정 온도일 때 값이 1인 횟수의 비율을 계산
• 즉, 반환된 모든 시뮬레이션의 평균, $P(Defect=1|t)$을 계산. 이렇게 해서 우리 데이터셋의 각 데이터 포인트에서 결함의 사후확률을 얻음
• 예를 들어 우리가 사용한 모델의 경우 코드는 다음과 같음

In [79]:

posterior_probability = simulations.mean(axis=0)

print("결함의 사후확률  | 실제 결함 ")
for i in range(len(D)):
    print("%.2f             |   %d" % (posterior_probability[i], D[i]))

결함의 사후확률  | 실제 결함 
0.42             |   0
0.24             |   1
0.27             |   0
0.32             |   0
0.37             |   0
0.18             |   0
0.15             |   0
0.24             |   0
0.76             |   1
0.57             |   1
0.24             |   1
0.08             |   0
0.37             |   0
0.83             |   1
0.37             |   0
0.12             |   0
0.24             |   0
0.06             |   0
0.10             |   0
0.07             |   0
0.12             |   1
0.10             |   0
0.74             |   1

• 다음으로 사후확률에 따라 각 columns를 정렬함

In [80]:

ix = np.argsort(posterior_probability)
print("probb | defect ")
for i in range(len(D)):
    print("%.2f  |   %d" % (posterior_probability[ix[i]], D[ix[i]]))

probb | defect 
0.06  |   0
0.07  |   0
0.08  |   0
0.10  |   0
0.10  |   0
0.12  |   0
0.12  |   1
0.15  |   0
0.18  |   0
0.24  |   0
0.24  |   0
0.24  |   1
0.24  |   1
0.27  |   0
0.32  |   0
0.37  |   0
0.37  |   0
0.37  |   0
0.42  |   0
0.57  |   1
0.74  |   1
0.76  |   1
0.83  |   1

• 이 데이터를 그림으로 더 잘 표현할 수 있음. 아래 그림처럼 데이터를 separation_plot 함수로 나타냄

In [ ]: