 Xem trên TensorFlow.org |  Chạy trong Google Colab |  Xem nguồn trên GitHub |  Tải xuống sổ ghi chép |
Vấn đề tám trường ( Rubin 1981 ) xem xét tính hiệu quả của các chương trình huấn luyện SAT tiến hành song song ở tám trường. Nó đã trở thành một vấn đề cổ điển ( phân tích dữ liệu Bayesian , Stan ) minh họa tính hữu ích của mô hình thứ bậc để chia sẻ thông tin giữa các nhóm trao đổi.
Các implemention dưới đây là một sự thích nghi của một Edward 1.0 hướng dẫn .
Nhập khẩu
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')
Dữ liệu
Từ Phân tích Dữ liệu Bayes, phần 5.5 (Gelman và cộng sự 2013):
Một nghiên cứu đã được thực hiện cho Cơ quan Khảo thí Giáo dục để phân tích tác động của các chương trình huấn luyện đặc biệt cho SAT-V (Scholastic Aptitude Test-Verbal) ở mỗi trường trong số tám trường trung học. Biến kết quả trong mỗi nghiên cứu là điểm số trong kỳ thi SAT-V đặc biệt, một bài kiểm tra trắc nghiệm tiêu chuẩn do Viện Khảo thí Giáo dục quản lý và được sử dụng để giúp các trường cao đẳng đưa ra quyết định tuyển sinh; điểm số có thể thay đổi từ 200 đến 800, với trung bình khoảng 500 và độ lệch chuẩn khoảng 100. Các kỳ thi SAT được thiết kế để chống lại những nỗ lực ngắn hạn hướng đặc biệt vào việc cải thiện thành tích trong bài kiểm tra; thay vào đó, chúng được thiết kế để phản ánh kiến thức có được và khả năng được phát triển trong nhiều năm giáo dục. Tuy nhiên, mỗi trường trong số tám trường trong nghiên cứu này đều coi chương trình huấn luyện ngắn hạn của mình là rất thành công trong việc tăng điểm SAT. Ngoài ra, không có lý do gì trước đó để tin rằng bất kỳ chương trình nào trong số tám chương trình này hiệu quả hơn chương trình nào khác hoặc một số chương trình có hiệu quả tương tự nhau hơn bất kỳ chương trình nào khác.
Đối với mỗi trong tám trường (\(J = 8\)), chúng ta có một ước tính hiệu quả điều trị \(y_j\) và sai số chuẩn của các ước tính ảnh hưởng \(\sigma_j\). Hiệu quả điều trị trong nghiên cứu thu được bằng hồi quy tuyến tính trên nhóm điều trị sử dụng điểm PSAT-M và PSAT-V làm biến kiểm soát. Như không có niềm tin trước rằng bất kỳ của các trường đã ít nhiều tương tự hoặc bất kỳ chương trình huấn luyện sẽ có hiệu quả hơn, chúng ta có thể xem xét các tác động điều trị như trao đổi .
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()
Mô hình
Để nắm bắt dữ liệu, chúng tôi sử dụng mô hình bình thường phân cấp. Nó tuân theo quy trình tổng hợp,
\[ \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*} \]
nơi \(\mu\) đại diện cho hiệu quả điều trị trung bình và các năm trước \(\tau\) điều khiển bao nhiêu sai có giữa các trường. Các \(y_i\) và \(\sigma_i\) được quan sát. Như \(\tau \rightarrow \infty\), mô hình tiếp cận các mô hình không-tổng hợp, tức là, mỗi người trong số các ước tính hiệu quả điều trị trường được phép độc lập hơn. Như \(\tau \rightarrow 0\), mô hình tiếp cận các mô hình hoàn chỉnh-tổng hợp, nghĩa là, tất cả các hiệu ứng điều trị trường là gần gũi hơn với các nhóm trung bình \(\mu\). Để hạn chế độ lệch chuẩn là tích cực, chúng ta rút ra \(\tau\) từ một phân phối lognormal (tương đương với vẽ \(log(\tau)\) từ một phân phối chuẩn).
Sau chẩn đoán Biased Suy luận với phân kỳ , chúng tôi chuyển đổi mô hình trên vào một mô hình phi tập trung tương đương:
\[ \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*} \]
Chúng tôi cụ thể hóa mô hình này như một JointDistributionSequential dụ:
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))
Suy luận Bayes
Dữ liệu được cung cấp, chúng tôi thực hiện Hamiltonian Monte Carlo (HMC) để tính toán phân phối sau trên các tham số của mô hình.
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()
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()
Chúng ta có thể quan sát sự thu hẹp về phía nhóm avg_effect trên.
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
Sự chỉ trích
Để có được sự phân bố tiên đoán sau, tức là một mô hình dữ liệu mới \(y^*\) cho các dữ liệu quan sát \(y\):
\[ p(y^*|y) \propto \int_\theta p(y^* | \theta)p(\theta |y)d\theta\]
chúng tôi ghi đè lên các giá trị của các biến ngẫu nhiên trong mô hình để đặt chúng vào giá trị trung bình của phân phối sau, và mẫu từ mô hình đó để tạo mới dữ liệu \(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()
# The mean predicted treatment effects for each of the eight schools.
prediction = np.mean(predictive_treatment_effects, axis=0)
Chúng ta có thể xem xét phần còn lại giữa dữ liệu hiệu quả điều trị và dự đoán của mô hình sau. Những điều này tương ứng với biểu đồ ở trên cho thấy sự thu hẹp của các tác động ước tính đối với mức trung bình của dân số.
treatment_effects - prediction
array([14.905351 , 1.2838383, -5.6966295, 0.8327627, -2.3356671,
-2.0363257, 5.997898 , 4.3731265], dtype=float32)
Bởi vì chúng ta có một phân phối các dự đoán cho mỗi trường, chúng ta có thể xem xét cả phân phối của phần dư.
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()
Sự nhìn nhận
Hướng dẫn này ban đầu được viết bằng Edward 1.0 ( nguồn ). Chúng tôi cảm ơn tất cả những người đóng góp để viết và sửa đổi phiên bản đó.
Người giới thiệu
- Donald B. Rubin. Ước lượng trong các thí nghiệm ngẫu nhiên song song. Tạp chí Thống kê Giáo dục, 6 (4): 377-401, 1981.
- Andrew Gelman, John Carlin, Hal Stern, David Dunson, Aki Vehtari và Donald Rubin. Phân tích dữ liệu Bayes, Ấn bản thứ ba. Chapman và Hall / CRC, 2013.