今日のローカルTensorFlowEverywhereイベントの出欠確認！
このページは Cloud Translation API によって翻訳されました。
Switch to English

8つの学校

TensorFlow.orgで表示 GoogleColabで実行 GitHubでソースを表示ノートブックをダウンロード

8つの学校の問題（ Rubin 1981 ）は、8つの学校で並行して実施されるSATコーチングプログラムの有効性を考慮しています。これは、交換可能なグループ間で情報を共有するための階層モデリングの有用性を示す古典的な問題（ Bayesian Data Analysis 、 Stan ）になっています。

以下の実装は、Edward1.0チュートリアルを応用したものです。

輸入

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

import tensorflow.compat.v2 as tf
import tensorflow_probability as tfp
from tensorflow_probability import distributions as tfd
import warnings

tf.enable_v2_behavior()

plt.style.use("ggplot")
warnings.filterwarnings('ignore')

データ

Bayesian Data Analysis、セクション5.5（Gelman etal。2013）から：

8つの高校のそれぞれでSAT-V（Scholastic Aptitude Test-Verbal）の特別なコーチングプログラムの効果を分析するために、Educational TestingServiceの調査が実施されました。各研究の結果変数は、SAT-Vの特別管理のスコアでした。これは、Educational Testing Serviceによって管理され、大学が入学決定を下すのを支援するために使用される標準化された多肢選択式テストです。スコアは200から800の間で変動し、平均は約500、標準偏差は約100です。SA​​T試験は、特にテストのパフォーマンスの向上に向けられた短期間の取り組みに耐えるように設計されています。代わりに、長年の教育で得られた知識と開発された能力を反映するように設計されています。それにもかかわらず、この研究の8つの学校のそれぞれは、その短期コーチングプログラムがSATスコアの増加に非常に成功していると考えました。また、8つのプログラムのいずれかが他のプログラムよりも効果的である、または一部のプログラムが他のプログラムよりも効果が類似していると信じる事前の理由はありませんでした。

8つの学校（$ J = 8 $）のそれぞれについて、推定治療効果$ y_j $と効果推定$ \ sigma_j $の標準誤差があります。この研究における治療効果は、PSAT-MおよびPSAT-Vスコアを制御変数として使用した治療群の線形回帰によって得られました。いずれの学校も多かれ少なかれ類似している、またはいずれかのコーチングプログラムがより効果的であるという事前の信念がなかったため、治療効果は交換可能であると見なすことができます。

num_schools = 8  # number of schools
treatment_effects = np.array(
    [28, 8, -3, 7, -1, 1, 18, 12], dtype=np.float32)  # treatment effects
treatment_stddevs = np.array(
    [15, 10, 16, 11, 9, 11, 10, 18], dtype=np.float32)  # treatment SE

fig, ax = plt.subplots()
plt.bar(range(num_schools), treatment_effects, yerr=treatment_stddevs)
plt.title("8 Schools treatment effects")
plt.xlabel("School")
plt.ylabel("Treatment effect")
fig.set_size_inches(10, 8)
plt.show()

png

モデル

データをキャプチャするために、階層正規モデルを使用します。それは生成プロセスに従います、

$$ \begin{align*} \mu &\sim \text{Normal}(\text{loc}{=}0,\ \text{scale}{=}10) \\ \log\tau &\sim \text{Normal}(\text{loc}{=}5,\ \text{scale}{=}1) \\ \text{for } & i=1\ldots 8:\\ & \theta_i \sim \text{Normal}\left(\text{loc}{=}\mu,\ \text{scale}{=}\tau \right) \\ & y_i \sim \text{Normal}\left(\text{loc}{=}\theta_i,\ \text{scale}{=}\sigma_i \right) \end{align*} $$

ここで、$ \ mu $は以前の平均治療効果を表し、$ \ tau $は学校間の差異の大きさを制御します。 $ y_i $と$ \ sigma_i $が監視されます。 $ \ tau \ rightarrow \ infty $として、モデルはプールなしモデルに近づきます。つまり、各学校の治療効果の推定値は、より独立したものになります。 $ \ tau \ rightarrow 0 $として、モデルは完全プーリングモデルに近づきます。つまり、すべての学校の治療効果はグループ平均$ \ mu $に近くなります。標準偏差を正に制限するために、対数正規分布から$ \ tau $を描画します（これは、正規分布から$ log（\ tau）$を描画するのと同じです）。

発散によるバイアス推論の診断に続いて、上記のモデルを同等の非中心モデルに変換します。

$$ \begin{align*} \mu &\sim \text{Normal}(\text{loc}{=}0,\ \text{scale}{=}10) \\ \log\tau &\sim \text{Normal}(\text{loc}{=}5,\ \text{scale}{=}1) \\ \text{for } & i=1\ldots 8:\\ & \theta_i' \sim \text{Normal}\left(\text{loc}{=}0,\ \text{scale}{=}1 \right) \\ & \theta_i = \mu + \tau \theta_i' \\ & y_i \sim \text{Normal}\left(\text{loc}{=}\theta_i,\ \text{scale}{=}\sigma_i \right) \end{align*} $$

このモデルをJointDistributionSequentialインスタンスとして具体化します。

model = tfd.JointDistributionSequential([
  tfd.Normal(loc=0., scale=10., name="avg_effect"),  # `mu` above
  tfd.Normal(loc=5., scale=1., name="avg_stddev"),  # `log(tau)` above
  tfd.Independent(tfd.Normal(loc=tf.zeros(num_schools),
                             scale=tf.ones(num_schools),
                             name="school_effects_standard"),  # `theta_prime` 
                  reinterpreted_batch_ndims=1),
  lambda school_effects_standard, avg_stddev, avg_effect: (
      tfd.Independent(tfd.Normal(loc=(avg_effect[..., tf.newaxis] +
                                      tf.exp(avg_stddev[..., tf.newaxis]) *
                                      school_effects_standard),  # `theta` above
                                 scale=treatment_stddevs),
                      name="treatment_effects",  # `y` above
                      reinterpreted_batch_ndims=1))
])

def target_log_prob_fn(avg_effect, avg_stddev, school_effects_standard):
  """Unnormalized target density as a function of states."""
  return model.log_prob((
      avg_effect, avg_stddev, school_effects_standard, treatment_effects))

ベイズ推定

データが与えられた場合、ハミルトニアンモンテカルロ（HMC）を実行して、モデルのパラメーターの事後分布を計算します。

num_results = 5000
num_burnin_steps = 3000

# Improve performance by tracing the sampler using `tf.function`
# and compiling it using XLA.
@tf.function(autograph=False, experimental_compile=True)
def do_sampling():
  return tfp.mcmc.sample_chain(
      num_results=num_results,
      num_burnin_steps=num_burnin_steps,
      current_state=[
          tf.zeros([], name='init_avg_effect'),
          tf.zeros([], name='init_avg_stddev'),
          tf.ones([num_schools], name='init_school_effects_standard'),
      ],
      kernel=tfp.mcmc.HamiltonianMonteCarlo(
          target_log_prob_fn=target_log_prob_fn,
          step_size=0.4,
          num_leapfrog_steps=3))

states, kernel_results = do_sampling()

avg_effect, avg_stddev, school_effects_standard = states

school_effects_samples = (
    avg_effect[:, np.newaxis] +
    np.exp(avg_stddev)[:, np.newaxis] * school_effects_standard)

num_accepted = np.sum(kernel_results.is_accepted)
print('Acceptance rate: {}'.format(num_accepted / num_results))

Acceptance rate: 0.5974


fig, axes = plt.subplots(8, 2, sharex='col', sharey='col')
fig.set_size_inches(12, 10)
for i in range(num_schools):
  axes[i][0].plot(school_effects_samples[:,i].numpy())
  axes[i][0].title.set_text("School {} treatment effect chain".format(i))
  sns.kdeplot(school_effects_samples[:,i].numpy(), ax=axes[i][1], shade=True)
  axes[i][1].title.set_text("School {} treatment effect distribution".format(i))
axes[num_schools - 1][0].set_xlabel("Iteration")
axes[num_schools - 1][1].set_xlabel("School effect")
fig.tight_layout()
plt.show()

png

print("E[avg_effect] = {}".format(np.mean(avg_effect)))
print("E[avg_stddev] = {}".format(np.mean(avg_stddev)))
print("E[school_effects_standard] =")
print(np.mean(school_effects_standard[:, ]))
print("E[school_effects] =")
print(np.mean(school_effects_samples[:, ], axis=0))

E[avg_effect] = 5.57183933258
E[avg_stddev] = 2.47738981247
E[school_effects_standard] =
0.08509017
E[school_effects] =
[15.0051     7.103311   2.4552586  6.2744603  1.3364682  3.1125953
 12.762501   7.743602 ]


# Compute the 95% interval for school_effects
school_effects_low = np.array([
    np.percentile(school_effects_samples[:, i], 2.5) for i in range(num_schools)
])
school_effects_med = np.array([
    np.percentile(school_effects_samples[:, i], 50) for i in range(num_schools)
])
school_effects_hi = np.array([
    np.percentile(school_effects_samples[:, i], 97.5)
    for i in range(num_schools)
])

fig, ax = plt.subplots(nrows=1, ncols=1, sharex=True)
ax.scatter(np.array(range(num_schools)), school_effects_med, color='red', s=60)
ax.scatter(
    np.array(range(num_schools)) + 0.1, treatment_effects, color='blue', s=60)

plt.plot([-0.2, 7.4], [np.mean(avg_effect),
                       np.mean(avg_effect)], 'k', linestyle='--')

ax.errorbar(
    np.array(range(8)),
    school_effects_med,
    yerr=[
        school_effects_med - school_effects_low,
        school_effects_hi - school_effects_med
    ],
    fmt='none')

ax.legend(('avg_effect', 'HMC', 'Observed effect'), fontsize=14)

plt.xlabel('School')
plt.ylabel('Treatment effect')
plt.title('HMC estimated school treatment effects vs. observed data')
fig.set_size_inches(10, 8)
plt.show()

png

上記のグループavg_effectへの縮小を観察できます。

print("Inferred posterior mean: {0:.2f}".format(
    np.mean(school_effects_samples[:,])))
print("Inferred posterior mean se: {0:.2f}".format(
    np.std(school_effects_samples[:,])))

Inferred posterior mean: 6.97
Inferred posterior mean se: 10.41

批判

事後予測分布、つまり、観測データ$ y $が与えられた場合の新しいデータ$ y ^ * $のモデルを取得するには：

$$ p(y^*|y) \propto \int_\theta p(y^* | \theta)p(\theta |y)d\theta$$

モデル内の確率変数の値をオーバーライドして事後分布の平均に設定し、そのモデルからサンプリングして新しいデータ$ y ^ * $を生成します。

sample_shape = [5000]

_, _, _, predictive_treatment_effects = model.sample(
    value=(tf.broadcast_to(np.mean(avg_effect, 0), sample_shape),
           tf.broadcast_to(np.mean(avg_stddev, 0), sample_shape),
           tf.broadcast_to(np.mean(school_effects_standard, 0),
                           sample_shape + [num_schools]),
           None))

fig, axes = plt.subplots(4, 2, sharex=True, sharey=True)
fig.set_size_inches(12, 10)
fig.tight_layout()
for i, ax in enumerate(axes):
  sns.kdeplot(predictive_treatment_effects[:, 2*i].numpy(),
              ax=ax[0], shade=True)
  ax[0].title.set_text(
      "School {} treatment effect posterior predictive".format(2*i))
  sns.kdeplot(predictive_treatment_effects[:, 2*i + 1].numpy(),
              ax=ax[1], shade=True)
  ax[1].title.set_text(
      "School {} treatment effect posterior predictive".format(2*i + 1))
plt.show()

png

# The mean predicted treatment effects for each of the eight schools.
prediction = np.mean(predictive_treatment_effects, axis=0)

治療効果データとモデル事後予測の間の残差を見ることができます。これらは、母集団の平均に対する推定効果の縮小を示す上記のプロットに対応しています。

treatment_effects - prediction

array([14.905351 ,  1.2838383, -5.6966295,  0.8327627, -2.3356671,
       -2.0363257,  5.997898 ,  4.3731265], dtype=float32)

学校ごとに予測の分布があるため、残差の分布も考慮することができます。

residuals = treatment_effects - predictive_treatment_effects

fig, axes = plt.subplots(4, 2, sharex=True, sharey=True)
fig.set_size_inches(12, 10)
fig.tight_layout()
for i, ax in enumerate(axes):
  sns.kdeplot(residuals[:, 2*i].numpy(), ax=ax[0], shade=True)
  ax[0].title.set_text(
      "School {} treatment effect residuals".format(2*i))
  sns.kdeplot(residuals[:, 2*i + 1].numpy(), ax=ax[1], shade=True)
  ax[1].title.set_text(
      "School {} treatment effect residuals".format(2*i + 1))
plt.show()

png

謝辞

このチュートリアルは、もともとEdward 1.0（ソース）で書かれました。そのバージョンの作成と改訂に貢献してくれたすべての人に感謝します。

参考文献

  1. ドナルド・B・ルービン。並行ランダム化実験での推定。 Journal of Educational Statistics、6（4）：377-401、1981。
  2. アンドリュー・ゲルマン、ジョン・カーリン、ハル・スターン、デビッド・ダンソン、アキ・ヴェタリ、ドナルド・ルービン。ベイジアンデータ分析、第3版。チャップマンアンドホール/ CRC、2013年。