این صفحه به‌وسیله ‏Cloud Translation API‏ ترجمه شده است.

هشت مدرسه

مشاهده در TensorFlow.org

در Google Colab اجرا شود

مشاهده منبع در GitHub

دانلود دفترچه یادداشت

مشکل هشت مدرسه ( روبین 1981 ) اثر بخشی برنامه های مربیگری شنبه به صورت موازی در هشت مدرسه انجام می داند. این تبدیل به یک مشکل کلاسیک ( بیزی تجزیه و تحلیل داده ، استن ) که نشان می دهد سودمندی مدل سازی سلسله مراتبی برای به اشتراک گذاری اطلاعات بین دو گروه تعویض.

implemention زیر یک اقتباس از ادوارد 1.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')

داده

از تجزیه و تحلیل داده های بیزی، بخش 5.5 (گلمن و همکاران 2013):

مطالعه ای برای سرویس تست آموزشی به منظور تجزیه و تحلیل اثرات برنامه های مربیگری ویژه برای SAT-V (آزمون استعداد تحصیلی-کلامی) در هر یک از هشت دبیرستان انجام شد. متغیر نتیجه در هر مطالعه، امتیاز یک مدیریت ویژه SAT-V بود، یک آزمون استاندارد چند گزینه‌ای که توسط سرویس تست آموزشی اجرا می‌شد و برای کمک به کالج‌ها در تصمیم‌گیری پذیرش استفاده می‌شد. نمرات می تواند بین 200 تا 800، با میانگین حدود 500 و انحراف استاندارد حدود 100 متغیر باشد. در عوض آنها به گونه ای طراحی شده اند که منعکس کننده دانش به دست آمده و توانایی های توسعه یافته در طی سالیان متمادی آموزش باشند. با این وجود، هر یک از هشت مدرسه در این مطالعه، برنامه مربیگری کوتاه مدت خود را در افزایش نمرات SAT بسیار موفق در نظر گرفتند. همچنین، هیچ دلیل قبلی برای این باور وجود نداشت که هر یک از هشت برنامه مؤثرتر از هر برنامه دیگری بوده یا اینکه برخی از نظر تأثیر به یکدیگر بیش از هر برنامه دیگری شبیه هستند.

برای هر یک از هشت مدرسه (\(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, jit_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

سپاسگزاریها

این آموزش در ابتدا در ادوارد 1.0 (نوشته شده بود منبع ). ما از همه مشارکت کنندگان در نوشتن و بازنگری آن نسخه تشکر می کنیم.

منابع

دونالد بی روبین. برآورد در آزمایش های تصادفی موازی. مجله آمار آموزشی، 6(4):377-401، 1360.
اندرو گلمن، جان کارلین، هال استرن، دیوید دانسون، آکی وهتری و دونالد روبین. تجزیه و تحلیل داده های بیزی، ویرایش سوم. چپمن و هال/CRC، 2013.