 TensorFlow.org पर देखें |  Google Colab में चलाएं |  GitHub पर स्रोत देखें |  नोटबुक डाउनलोड करें |
1 परिचय
इस कोलाब में हम एक लोकप्रिय, टॉय डेटासेट के लिए एक रेखीय मिश्रित-प्रभाव प्रतिगमन मॉडल फिट करेंगे। हम इस फिट तीन बार बनाने, अनुसंधान का उपयोग होगा lme4 , स्टेन का मिश्रित प्रभाव पैकेज, और TensorFlow संभावना (टीएफपी) पुरातन। हम तीनों को मोटे तौर पर एक ही फिट किए गए मापदंडों और पश्च वितरण को दिखाते हुए निष्कर्ष निकालते हैं।
हमारा मुख्य निष्कर्ष यह है कि TFP फिट करने के लिए मॉडल hlm-तरह और यह परिणाम है जो अन्य सॉफ्टवेयर संकुल के साथ संगत कर रहे हैं, अर्थात्।, पैदा करता है कि आवश्यक सामान्य टुकड़े है lme4 , rstanarm । यह कोलाब तुलना किए गए किसी भी पैकेज की कम्प्यूटेशनल दक्षता का सटीक प्रतिबिंब नहीं है।
%matplotlib inline
import os
from six.moves import urllib
import numpy as np
import pandas as pd
import warnings
from matplotlib import pyplot as plt
import seaborn as sns
from IPython.core.pylabtools import figsize
figsize(11, 9)
import tensorflow.compat.v1 as tf
import tensorflow_datasets as tfds
import tensorflow_probability as tfp
2 पदानुक्रमित रैखिक मॉडल
आर, स्टेन, और TFP के बीच हमारी तुलना के लिए, हम एक फिट होगा श्रेणीबद्ध लीनियर मॉडल के लिए (HLM) रेडॉन डाटासेट में लोकप्रिय बना दिया Gelman, एट द्वारा बायेसियन डेटा विश्लेषण। अल. (पृष्ठ 559, दूसरा संस्करण; पृष्ठ 250, तीसरा संस्करण।)।
हम निम्नलिखित जनरेटिव मॉडल मानते हैं:
\[\begin{align*} \text{for } & c=1\ldots \text{NumCounties}:\\ & \beta_c \sim \text{Normal}\left(\text{loc}=0, \text{scale}=\sigma_C \right) \\ \text{for } & i=1\ldots \text{NumSamples}:\\ &\eta_i = \underbrace{\omega_0 + \omega_1 \text{Floor}_i}_\text{fixed effects} + \underbrace{\beta_{ \text{County}_i} \log( \text{UraniumPPM}_{\text{County}_i}))}_\text{random effects} \\ &\log(\text{Radon}_i) \sim \text{Normal}(\text{loc}=\eta_i , \text{scale}=\sigma_N) \end{align*}\]
आर के दशक में lme4 "अंकन टिल्ड", इस मॉडल के बराबर है:
log_radon ~ 1 + floor + (0 + log_uranium_ppm | county)
हम के लिए MLE मिलेगा \(\omega, \sigma_C, \sigma_N\) का पिछला वितरण (साक्ष्य के वातानुकूलित) का उपयोग कर \(\{\beta_c\}_{c=1}^\text{NumCounties}\)।
अनिवार्य रूप से एक ही मॉडल के लिए, लेकिन एक यादृच्छिक अवरोधन के साथ, देखने के परिशिष्ट ए ।
HLMs की एक अधिक सामान्य विनिर्देश के लिए, देखें परिशिष्ट बी ।
3 डेटा मुंगिंग
इस भाग में हम प्राप्त radon डाटासेट और यह हमारे ग्रहण मॉडल के अनुरूप हो कुछ न्यूनतम preprocessing है।
def load_and_preprocess_radon_dataset(state='MN'):
"""Preprocess Radon dataset as done in "Bayesian Data Analysis" book.
We filter to Minnesota data (919 examples) and preprocess to obtain the
following features:
- `log_uranium_ppm`: Log of soil uranium measurements.
- `county`: Name of county in which the measurement was taken.
- `floor`: Floor of house (0 for basement, 1 for first floor) on which the
measurement was taken.
The target variable is `log_radon`, the log of the Radon measurement in the
house.
"""
ds = tfds.load('radon', split='train')
radon_data = tfds.as_dataframe(ds)
radon_data.rename(lambda s: s[9:] if s.startswith('feat') else s, axis=1, inplace=True)
df = radon_data[radon_data.state==state.encode()].copy()
# For any missing or invalid activity readings, we'll use a value of `0.1`.
df['radon'] = df.activity.apply(lambda x: x if x > 0. else 0.1)
# Make county names look nice.
df['county'] = df.county.apply(lambda s: s.decode()).str.strip().str.title()
# Remap categories to start from 0 and end at max(category).
county_name = sorted(df.county.unique())
df['county'] = df.county.astype(
pd.api.types.CategoricalDtype(categories=county_name)).cat.codes
county_name = list(map(str.strip, county_name))
df['log_radon'] = df['radon'].apply(np.log)
df['log_uranium_ppm'] = df['Uppm'].apply(np.log)
df = df[['idnum', 'log_radon', 'floor', 'county', 'log_uranium_ppm']]
return df, county_name
radon, county_name = load_and_preprocess_radon_dataset()
# We'll use the following directory to store our preprocessed dataset.
CACHE_DIR = os.path.join(os.sep, 'tmp', 'radon')
# Save processed data. (So we can later read it in R.)
if not tf.gfile.Exists(CACHE_DIR):
tf.gfile.MakeDirs(CACHE_DIR)
with tf.gfile.Open(os.path.join(CACHE_DIR, 'radon.csv'), 'w') as f:
radon.to_csv(f, index=False)
3.1 अपना डेटा जानें
इस भाग में हम पता लगाने radon क्यों प्रस्तावित मॉडल उचित हो सकता है की एक बेहतर समझ प्राप्त करने के लिए डाटासेट।
radon.head()
fig, ax = plt.subplots(figsize=(22, 5));
county_freq = radon['county'].value_counts()
county_freq.plot(kind='bar', color='#436bad');
plt.xlabel('County index')
plt.ylabel('Number of radon readings')
plt.title('Number of radon readings per county', fontsize=16)
county_freq = np.array(zip(county_freq.index, county_freq.values)) # We'll use this later.
fig, ax = plt.subplots(ncols=2, figsize=[10, 4]);
radon['log_radon'].plot(kind='density', ax=ax[0]);
ax[0].set_xlabel('log(radon)')
radon['floor'].value_counts().plot(kind='bar', ax=ax[1]);
ax[1].set_xlabel('Floor');
ax[1].set_ylabel('Count');
fig.subplots_adjust(wspace=0.25)
निष्कर्ष:
- 85 काउंटियों की एक लंबी पूंछ है। (जीएलएमएम में एक सामान्य घटना।)
- दरअसल \(\log(\text{Radon})\) स्वेच्छापूर्ण है। (इसलिए रैखिक प्रतिगमन समझ में आ सकता है।)
- रीडिंग सबसे पर किया जाता है \(0\)वें मंजिल; कोई पढ़ने मंजिल से ऊपर बनाया गया था \(1\)। (इसलिए हमारे निश्चित प्रभावों के केवल दो भार होंगे।)
4 एचएलएम आर . में
इस भाग में हम अनुसंधान का उपयोग lme4 पैकेज संभाव्य मॉडल ऊपर वर्णित फिट करने के लिए।
suppressMessages({
library('bayesplot')
library('data.table')
library('dplyr')
library('gfile')
library('ggplot2')
library('lattice')
library('lme4')
library('plyr')
library('rstanarm')
library('tidyverse')
RequireInitGoogle()
})
data = read_csv(gfile::GFile('/tmp/radon/radon.csv'))
Parsed with column specification: cols( log_radon = col_double(), floor = col_integer(), county = col_integer(), log_uranium_ppm = col_double() )
head(data)
# A tibble: 6 x 4
log_radon floor county log_uranium_ppm
<dbl> <int> <int> <dbl>
1 0.788 1 0 -0.689
2 0.788 0 0 -0.689
3 1.06 0 0 -0.689
4 0 0 0 -0.689
5 1.13 0 1 -0.847
6 0.916 0 1 -0.847
# https://github.com/stan-dev/example-models/wiki/ARM-Models-Sorted-by-Chapter
radon.model <- lmer(log_radon ~ 1 + floor + (0 + log_uranium_ppm | county), data = data)
summary(radon.model)
Linear mixed model fit by REML ['lmerMod']
Formula: log_radon ~ 1 + floor + (0 + log_uranium_ppm | county)
Data: data
REML criterion at convergence: 2166.3
Scaled residuals:
Min 1Q Median 3Q Max
-4.5202 -0.6064 0.0107 0.6334 3.4111
Random effects:
Groups Name Variance Std.Dev.
county log_uranium_ppm 0.7545 0.8686
Residual 0.5776 0.7600
Number of obs: 919, groups: county, 85
Fixed effects:
Estimate Std. Error t value
(Intercept) 1.47585 0.03899 37.85
floor -0.67974 0.06963 -9.76
Correlation of Fixed Effects:
(Intr)
floor -0.330
qqmath(ranef(radon.model, condVar=TRUE))
$county
write.csv(as.data.frame(ranef(radon.model, condVar = TRUE)), '/tmp/radon/lme4_fit.csv')
स्टेन में 5 एचएलएम
इस भाग में हम का उपयोग rstanarm के रूप में एक ही सूत्र / सिंटैक्स का उपयोग एक स्टेन मॉडल फिट करने के लिए lme4 ऊपर मॉडल।
विपरीत lme4 और नीचे TF मॉडल, rstanarm एक पूरी तरह से बायेसियन मॉडल, खुद को एक वितरण से तैयार मानकों के साथ एक सामान्य वितरण से तैयार यानी, सभी मापदंडों माना जाता है।
fit <- stan_lmer(log_radon ~ 1 + floor + (0 + log_uranium_ppm | county), data = data)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1, Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1, Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1, Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1, Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1, Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1, Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1, Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1, Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1, Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1, Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1, Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1, Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 7.73495 seconds (Warm-up)
2.98852 seconds (Sampling)
10.7235 seconds (Total)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2, Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2, Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2, Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2, Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2, Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2, Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2, Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2, Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2, Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2, Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2, Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2, Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 7.51252 seconds (Warm-up)
3.08653 seconds (Sampling)
10.5991 seconds (Total)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3, Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3, Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3, Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3, Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3, Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3, Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3, Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3, Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3, Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3, Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3, Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3, Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 8.14628 seconds (Warm-up)
3.01001 seconds (Sampling)
11.1563 seconds (Total)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
Chain 4, Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4, Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4, Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4, Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4, Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4, Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4, Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4, Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4, Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4, Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4, Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4, Iteration: 2000 / 2000 [100%] (Sampling)
Elapsed Time: 7.6801 seconds (Warm-up)
3.23663 seconds (Sampling)
10.9167 seconds (Total)
fit
stan_lmer(formula = log_radon ~ 1 + floor + (0 + log_uranium_ppm |
county), data = data)
Estimates:
Median MAD_SD
(Intercept) 1.5 0.0
floor -0.7 0.1
sigma 0.8 0.0
Error terms:
Groups Name Std.Dev.
county log_uranium_ppm 0.87
Residual 0.76
Num. levels: county 85
Sample avg. posterior predictive
distribution of y (X = xbar):
Median MAD_SD
mean_PPD 1.2 0.0
Observations: 919 Number of unconstrained parameters: 90
color_scheme_set("red")
ppc_dens_overlay(y = fit$y,
yrep = posterior_predict(fit, draws = 50))
color_scheme_set("brightblue")
ppc_intervals(
y = data$log_radon,
yrep = posterior_predict(fit),
x = data$county,
prob = 0.8
) +
labs(
x = "County",
y = "log radon",
title = "80% posterior predictive intervals \nvs observed log radon",
subtitle = "by county"
) +
panel_bg(fill = "gray95", color = NA) +
grid_lines(color = "white")
# Write the posterior samples (4000 for each variable) to a CSV.
write.csv(tidy(as.matrix(fit)), "/tmp/radon/stan_fit.csv")
with tf.gfile.Open('/tmp/radon/lme4_fit.csv', 'r') as f:
lme4_fit = pd.read_csv(f, index_col=0)
lme4_fit.head()
बाद में विज़ुअलाइज़ेशन के लिए lme4 से समूह यादृच्छिक प्रभावों के लिए बिंदु अनुमान और सशर्त मानक विचलन पुनर्प्राप्त करें।
posterior_random_weights_lme4 = np.array(lme4_fit.condval, dtype=np.float32)
lme4_prior_scale = np.array(lme4_fit.condsd, dtype=np.float32)
print(posterior_random_weights_lme4.shape, lme4_prior_scale.shape)
(85,) (85,)
lme4 अनुमानित साधनों और मानक विचलन का उपयोग करके काउंटी भार के लिए नमूने बनाएं।
with tf.Session() as sess:
lme4_dist = tfp.distributions.Independent(
tfp.distributions.Normal(
loc=posterior_random_weights_lme4,
scale=lme4_prior_scale),
reinterpreted_batch_ndims=1)
posterior_random_weights_lme4_final_ = sess.run(lme4_dist.sample(4000))
posterior_random_weights_lme4_final_.shape
(4000, 85)
हम स्टेन फिट से काउंटी भार के पीछे के नमूने भी प्राप्त करते हैं।
with tf.gfile.Open('/tmp/radon/stan_fit.csv', 'r') as f:
samples = pd.read_csv(f, index_col=0)
samples.head()
posterior_random_weights_cols = [
col for col in samples.columns if 'b.log_uranium_ppm.county' in col
]
posterior_random_weights_final_stan = samples[
posterior_random_weights_cols].values
print(posterior_random_weights_final_stan.shape)
(4000, 85)
यह स्टेन उदाहरण से पता चलता है कि कैसे एक TFP, यानी, के करीब सीधे संभाव्य मॉडल का उल्लेख करके शैली में LMER लागू करेगा।
टीएफ संभावना में 6 एचएलएम
इस भाग में हम निम्न स्तर TensorFlow संभावना पुरातन (का उपयोग करेगा Distributions हमारे श्रेणीबद्ध लीनियर मॉडल निर्दिष्ट करने के लिए) के साथ ही अज्ञात पैरामीटर के अनुरूप हैं।
# Handy snippet to reset the global graph and global session.
with warnings.catch_warnings():
warnings.simplefilter('ignore')
tf.reset_default_graph()
try:
sess.close()
except:
pass
sess = tf.InteractiveSession()
6.1 मॉडल निर्दिष्ट करें
इस भाग में हम निर्दिष्ट मिश्रित प्रभाव मॉडल रैखिक राडोण TFP पुरातन का उपयोग कर। ऐसा करने के लिए, हम दो फ़ंक्शन निर्दिष्ट करते हैं जो दो TFP वितरण उत्पन्न करते हैं:
-
make_weights_prior: यादृच्छिक वजन (जो से गुणा कर रहे हैं के लिए एक मल्टीवेरिएट सामान्य पूर्व \(\log(\text{UraniumPPM}_{c_i})\) रैखिक भविष्यवक्ता compue के लिए)। -
make_log_radon_likelihood: का एक बैचNormalप्रत्येक मनाया से अधिक वितरण \(\log(\text{Radon}_i)\) निर्भर चर।
हम इन वितरण हम TF चर (यानी, का उपयोग करना चाहिए में से प्रत्येक के मापदंडों फिटिंग हो जाएगा के बाद से tf.get_variable )। हालाँकि, चूंकि हम अप्रतिबंधित अनुकूलन का उपयोग करना चाहते हैं, इसलिए हमें आवश्यक शब्दार्थ प्राप्त करने के लिए वास्तविक-मूल्यों को बाधित करने का एक तरीका खोजना चाहिए, उदाहरण के लिए, मानक विचलन का प्रतिनिधित्व करने वाले पोस्टिव्स।
inv_scale_transform = lambda y: np.log(y) # Not using TF here.
fwd_scale_transform = tf.exp
निम्नलिखित समारोह हमारे पूर्व, निर्माण करती \(p(\beta|\sigma_C)\) जहां \(\beta\) यादृच्छिक प्रभाव वजन और अर्थ है \(\sigma_C\) मानक विचलन।
हम का उपयोग tf.make_template सुनिश्चित करना है कि यह कार्य करने के लिए पहली कॉल TF चर इसे इस्तेमाल करता है और चर के वर्तमान मूल्य का पुन: उपयोग बाद के सभी कॉल को दर्शाता है।
def _make_weights_prior(num_counties, dtype):
"""Returns a `len(log_uranium_ppm)` batch of univariate Normal."""
raw_prior_scale = tf.get_variable(
name='raw_prior_scale',
initializer=np.array(inv_scale_transform(1.), dtype=dtype))
return tfp.distributions.Independent(
tfp.distributions.Normal(
loc=tf.zeros(num_counties, dtype=dtype),
scale=fwd_scale_transform(raw_prior_scale)),
reinterpreted_batch_ndims=1)
make_weights_prior = tf.make_template(
name_='make_weights_prior', func_=_make_weights_prior)
निम्नलिखित समारोह हमारे संभावना, निर्माण करती \(p(y|x,\omega,\beta,\sigma_N)\) जहां \(y,x\) निरूपित प्रतिक्रिया और सबूत, \(\omega,\beta\) निरूपित फिक्स्ड और यादृच्छिक प्रभाव भार, और \(\sigma_N\) मानक विचलन।
यहाँ फिर से उपयोग हम tf.make_template सुनिश्चित करने के लिए TF चर कॉल भर में पुन: उपयोग किया जाता है।
def _make_log_radon_likelihood(random_effect_weights, floor, county,
log_county_uranium_ppm, init_log_radon_stddev):
raw_likelihood_scale = tf.get_variable(
name='raw_likelihood_scale',
initializer=np.array(
inv_scale_transform(init_log_radon_stddev), dtype=dtype))
fixed_effect_weights = tf.get_variable(
name='fixed_effect_weights', initializer=np.array([0., 1.], dtype=dtype))
fixed_effects = fixed_effect_weights[0] + fixed_effect_weights[1] * floor
random_effects = tf.gather(
random_effect_weights * log_county_uranium_ppm,
indices=tf.to_int32(county),
axis=-1)
linear_predictor = fixed_effects + random_effects
return tfp.distributions.Normal(
loc=linear_predictor, scale=fwd_scale_transform(raw_likelihood_scale))
make_log_radon_likelihood = tf.make_template(
name_='make_log_radon_likelihood', func_=_make_log_radon_likelihood)
अंत में हम संयुक्त लॉग-घनत्व के निर्माण के लिए पूर्व और संभावना जनरेटर का उपयोग करते हैं।
def joint_log_prob(random_effect_weights, log_radon, floor, county,
log_county_uranium_ppm, dtype):
num_counties = len(log_county_uranium_ppm)
rv_weights = make_weights_prior(num_counties, dtype)
rv_radon = make_log_radon_likelihood(
random_effect_weights,
floor,
county,
log_county_uranium_ppm,
init_log_radon_stddev=radon.log_radon.values.std())
return (rv_weights.log_prob(random_effect_weights)
+ tf.reduce_sum(rv_radon.log_prob(log_radon), axis=-1))
6.2 प्रशिक्षण (उम्मीदों को अधिकतम करने का स्टोकेस्टिक अनुमान)
हमारे रैखिक मिश्रित-प्रभाव प्रतिगमन मॉडल को फिट करने के लिए, हम एक्सपेक्टेशन मैक्सिमाइज़ेशन एल्गोरिथम (SAEM) के एक स्टोकेस्टिक सन्निकटन संस्करण का उपयोग करेंगे। मूल विचार अपेक्षित संयुक्त लॉग-घनत्व (ई-स्टेप) का अनुमान लगाने के लिए पश्च से नमूनों का उपयोग करना है। फिर हम उन मापदंडों को ढूंढते हैं जो इस गणना (एम-स्टेप) को अधिकतम करते हैं। कुछ अधिक ठोस रूप से, निश्चित-बिंदु पुनरावृत्ति द्वारा दिया गया है:
\[\begin{align*} \text{E}[ \log p(x, Z | \theta) | \theta_0] &\approx \frac{1}{M} \sum_{m=1}^M \log p(x, z_m | \theta), \quad Z_m\sim p(Z | x, \theta_0) && \text{E-step}\\ &=: Q_M(\theta, \theta_0) \\ \theta_0 &= \theta_0 - \eta \left.\nabla_\theta Q_M(\theta, \theta_0)\right|_{\theta=\theta_0} && \text{M-step} \end{align*}\]
जहां \(x\) सबूत, को दर्शाता है \(Z\) कुछ अव्यक्त चर जो बाहर हाशिए पर किए जाने की जरूरत है, और \(\theta,\theta_0\) संभव parameterizations।
एक अधिक विस्तृत विवरण के लिए, देखें बर्नार्ड Delyon, मार्क Lavielle, एरिक, Moulines (एन। सांख्यिकीविद।, 1999) द्वारा ईएम एल्गोरिदम के एक स्टोकेस्टिक सन्निकटन संस्करण का अभिसरण ।
ई-स्टेप की गणना करने के लिए, हमें पश्च से नमूना लेना होगा। चूंकि हमारे पोस्टीरियर से नमूना लेना आसान नहीं है, इसलिए हम हैमिल्टनियन मोंटे कार्लो (HMC) का उपयोग करते हैं। एचएमसी एक मोंटे कार्लो मार्कोव चेन प्रक्रिया है जो नए नमूनों का प्रस्ताव करने के लिए असामान्य पोस्टीरियर लॉग-घनत्व के ग्रेडियेंट (wrt राज्य, पैरामीटर नहीं) का उपयोग करती है।
असामान्य पोस्टीरियर लॉग-घनत्व को निर्दिष्ट करना सरल है - यह केवल संयुक्त लॉग-घनत्व "पिन किया गया" है जिस पर हम शर्त लगाना चाहते हैं।
# Specify unnormalized posterior.
dtype = np.float32
log_county_uranium_ppm = radon[
['county', 'log_uranium_ppm']].drop_duplicates().values[:, 1]
log_county_uranium_ppm = log_county_uranium_ppm.astype(dtype)
def unnormalized_posterior_log_prob(random_effect_weights):
return joint_log_prob(
random_effect_weights=random_effect_weights,
log_radon=dtype(radon.log_radon.values),
floor=dtype(radon.floor.values),
county=np.int32(radon.county.values),
log_county_uranium_ppm=log_county_uranium_ppm,
dtype=dtype)
अब हम एचएमसी ट्रांजिशन कर्नेल बनाकर ई-स्टेप सेटअप को पूरा करते हैं।
टिप्पणियाँ:
हम का उपयोग
state_stop_gradient=Trueके माध्यम से backpropping से एम-कदम एमसीएमसी से ड्रॉ को रोकने के लिए। (याद है, हम के माध्यम से backprop की जरूरत नहीं है क्योंकि हमारे ई-कदम जानबूझकर पिछले सबसे अच्छा ज्ञात आकलनकर्ता पर parameterized है।)हम का उपयोग
tf.placeholderइतनी है कि जब हम अंततः हमारी TF ग्राफ निष्पादित, हम अगले चरण की श्रृंखला के मूल्य के रूप में पिछले यात्रा के यादृच्छिक एमसीएमसी नमूना फ़ीड कर सकते हैं।हम TFP के अनुकूली का उपयोग
step_sizeअनुमानी,tfp.mcmc.hmc_step_size_update_fn।
# Set-up E-step.
step_size = tf.get_variable(
'step_size',
initializer=np.array(0.2, dtype=dtype),
trainable=False)
hmc = tfp.mcmc.HamiltonianMonteCarlo(
target_log_prob_fn=unnormalized_posterior_log_prob,
num_leapfrog_steps=2,
step_size=step_size,
step_size_update_fn=tfp.mcmc.make_simple_step_size_update_policy(
num_adaptation_steps=None),
state_gradients_are_stopped=True)
init_random_weights = tf.placeholder(dtype, shape=[len(log_county_uranium_ppm)])
posterior_random_weights, kernel_results = tfp.mcmc.sample_chain(
num_results=3,
num_burnin_steps=0,
num_steps_between_results=0,
current_state=init_random_weights,
kernel=hmc)
अब हम एम-स्टेप सेट-अप करते हैं। यह अनिवार्य रूप से TF में किए जा सकने वाले अनुकूलन के समान है।
# Set-up M-step.
loss = -tf.reduce_mean(kernel_results.accepted_results.target_log_prob)
global_step = tf.train.get_or_create_global_step()
learning_rate = tf.train.exponential_decay(
learning_rate=0.1,
global_step=global_step,
decay_steps=2,
decay_rate=0.99)
optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate)
train_op = optimizer.minimize(loss, global_step=global_step)
हम कुछ हाउसकीपिंग कार्यों के साथ समाप्त करते हैं। हमें TF को बताना चाहिए कि सभी वेरिएबल्स को इनिशियलाइज़ किया गया है। हम भी हमारे TF चर के हैंडल बनाने के तो हम कर सकते हैं print प्रक्रिया के प्रत्येक चरण पर उनके मूल्यों।
# Initialize all variables.
init_op = tf.initialize_all_variables()
# Grab variable handles for diagnostic purposes.
with tf.variable_scope('make_weights_prior', reuse=True):
prior_scale = fwd_scale_transform(tf.get_variable(
name='raw_prior_scale', dtype=dtype))
with tf.variable_scope('make_log_radon_likelihood', reuse=True):
likelihood_scale = fwd_scale_transform(tf.get_variable(
name='raw_likelihood_scale', dtype=dtype))
fixed_effect_weights = tf.get_variable(
name='fixed_effect_weights', dtype=dtype)
6.3 निष्पादित
इस खंड में हम अपने एसएईएम टीएफ ग्राफ को निष्पादित करते हैं। यहां मुख्य चाल एचएमसी कर्नेल से हमारे अंतिम ड्रॉ को अगले पुनरावृत्ति में फीड करना है। इस के बारे में हमारी प्रयोग के माध्यम से हासिल की है feed_dict में sess.run कॉल।
init_op.run()
w_ = np.zeros([len(log_county_uranium_ppm)], dtype=dtype)
%%time
maxiter = int(1500)
num_accepted = 0
num_drawn = 0
for i in range(maxiter):
[
_,
global_step_,
loss_,
posterior_random_weights_,
kernel_results_,
step_size_,
prior_scale_,
likelihood_scale_,
fixed_effect_weights_,
] = sess.run([
train_op,
global_step,
loss,
posterior_random_weights,
kernel_results,
step_size,
prior_scale,
likelihood_scale,
fixed_effect_weights,
], feed_dict={init_random_weights: w_})
w_ = posterior_random_weights_[-1, :]
num_accepted += kernel_results_.is_accepted.sum()
num_drawn += kernel_results_.is_accepted.size
acceptance_rate = num_accepted / num_drawn
if i % 100 == 0 or i == maxiter - 1:
print('global_step:{:>4} loss:{: 9.3f} acceptance:{:.4f} '
'step_size:{:.4f} prior_scale:{:.4f} likelihood_scale:{:.4f} '
'fixed_effect_weights:{}'.format(
global_step_, loss_.mean(), acceptance_rate, step_size_,
prior_scale_, likelihood_scale_, fixed_effect_weights_))
global_step: 0 loss: 1966.948 acceptance:1.0000 step_size:0.2000 prior_scale:1.0000 likelihood_scale:0.8529 fixed_effect_weights:[ 0. 1.] global_step: 100 loss: 1165.385 acceptance:0.6205 step_size:0.2040 prior_scale:0.9568 likelihood_scale:0.7611 fixed_effect_weights:[ 1.47523439 -0.66043079] global_step: 200 loss: 1149.851 acceptance:0.6766 step_size:0.2081 prior_scale:0.7465 likelihood_scale:0.7665 fixed_effect_weights:[ 1.48918796 -0.67058587] global_step: 300 loss: 1163.464 acceptance:0.6811 step_size:0.2040 prior_scale:0.8445 likelihood_scale:0.7594 fixed_effect_weights:[ 1.46291411 -0.67586178] global_step: 400 loss: 1158.846 acceptance:0.6808 step_size:0.2081 prior_scale:0.8377 likelihood_scale:0.7574 fixed_effect_weights:[ 1.47349834 -0.68823022] global_step: 500 loss: 1154.193 acceptance:0.6766 step_size:0.1961 prior_scale:0.8546 likelihood_scale:0.7564 fixed_effect_weights:[ 1.47703862 -0.67521363] global_step: 600 loss: 1163.903 acceptance:0.6783 step_size:0.2040 prior_scale:0.9491 likelihood_scale:0.7587 fixed_effect_weights:[ 1.48268366 -0.69667786] global_step: 700 loss: 1163.894 acceptance:0.6767 step_size:0.1961 prior_scale:0.8644 likelihood_scale:0.7617 fixed_effect_weights:[ 1.4719094 -0.66897118] global_step: 800 loss: 1153.689 acceptance:0.6742 step_size:0.2123 prior_scale:0.8366 likelihood_scale:0.7609 fixed_effect_weights:[ 1.47345769 -0.68343043] global_step: 900 loss: 1155.312 acceptance:0.6718 step_size:0.2040 prior_scale:0.8633 likelihood_scale:0.7581 fixed_effect_weights:[ 1.47426116 -0.6748783 ] global_step:1000 loss: 1151.278 acceptance:0.6690 step_size:0.2081 prior_scale:0.8737 likelihood_scale:0.7581 fixed_effect_weights:[ 1.46990883 -0.68891817] global_step:1100 loss: 1156.858 acceptance:0.6676 step_size:0.2040 prior_scale:0.8716 likelihood_scale:0.7584 fixed_effect_weights:[ 1.47386014 -0.6796245 ] global_step:1200 loss: 1166.247 acceptance:0.6653 step_size:0.2000 prior_scale:0.8748 likelihood_scale:0.7588 fixed_effect_weights:[ 1.47389269 -0.67626756] global_step:1300 loss: 1165.263 acceptance:0.6636 step_size:0.2040 prior_scale:0.8771 likelihood_scale:0.7581 fixed_effect_weights:[ 1.47612262 -0.67752427] global_step:1400 loss: 1158.108 acceptance:0.6640 step_size:0.2040 prior_scale:0.8748 likelihood_scale:0.7587 fixed_effect_weights:[ 1.47534692 -0.6789524 ] global_step:1499 loss: 1161.030 acceptance:0.6638 step_size:0.1941 prior_scale:0.8738 likelihood_scale:0.7589 fixed_effect_weights:[ 1.47624075 -0.67875224] CPU times: user 1min 16s, sys: 17.6 s, total: 1min 33s Wall time: 27.9 s
ऐसा लगता है कि ~1500 कदमों के बाद, मापदंडों के हमारे अनुमान स्थिर हो गए हैं।
6.4 परिणाम
अब जबकि हम मापदंडों को फिट कर चुके हैं, आइए बड़ी संख्या में पश्च नमूने तैयार करें और परिणामों का अध्ययन करें।
%%time
posterior_random_weights_final, kernel_results_final = tfp.mcmc.sample_chain(
num_results=int(15e3),
num_burnin_steps=int(1e3),
current_state=init_random_weights,
kernel=tfp.mcmc.HamiltonianMonteCarlo(
target_log_prob_fn=unnormalized_posterior_log_prob,
num_leapfrog_steps=2,
step_size=step_size))
[
posterior_random_weights_final_,
kernel_results_final_,
] = sess.run([
posterior_random_weights_final,
kernel_results_final,
], feed_dict={init_random_weights: w_})
CPU times: user 1min 42s, sys: 26.6 s, total: 2min 8s Wall time: 35.1 s
print('prior_scale: ', prior_scale_)
print('likelihood_scale: ', likelihood_scale_)
print('fixed_effect_weights: ', fixed_effect_weights_)
print('acceptance rate final: ', kernel_results_final_.is_accepted.mean())
prior_scale: 0.873799 likelihood_scale: 0.758913 fixed_effect_weights: [ 1.47624075 -0.67875224] acceptance rate final: 0.7448
अब हम एक बॉक्स और गलमुच्छा के चित्र का निर्माण \(\beta_c \log(\text{UraniumPPM}_c)\) यादृच्छिक प्रभाव। हम काउंटी आवृत्ति को कम करके यादृच्छिक प्रभावों का आदेश देंगे।
x = posterior_random_weights_final_ * log_county_uranium_ppm
I = county_freq[:, 0]
x = x[:, I]
cols = np.array(county_name)[I]
pw = pd.DataFrame(x)
pw.columns = cols
fig, ax = plt.subplots(figsize=(25, 4))
ax = pw.boxplot(rot=80, vert=True);
इस बॉक्स और गलमुच्छा चित्र से, हम उस काउंटी स्तरीय के विचरण का निरीक्षण \(\log(\text{UraniumPPM})\) काउंटी के रूप में यादृच्छिक प्रभाव बढ़ जाती है कम डेटासेट में प्रतिनिधित्व किया है। सहज रूप से यह समझ में आता है - अगर हमारे पास इसके लिए कम सबूत हैं तो हमें एक निश्चित काउंटी के प्रभाव के बारे में कम निश्चित होना चाहिए।
7 अगल-बगल तुलना
अब हम तीनों प्रक्रियाओं के परिणामों की तुलना करते हैं। ऐसा करने के लिए, हम स्टैन और टीएफपी द्वारा उत्पन्न पश्च नमूनों के गैर-पैरामीटर अनुमानों की गणना करेंगे। हम यह भी parameteric (अनुमानित) आर के द्वारा उत्पन्न अनुमान के खिलाफ की तुलना करेंगे lme4 पैकेज।
निम्नलिखित कथानक मिनेसोटा में प्रत्येक काउंटी के लिए प्रत्येक भार के पश्च वितरण को दर्शाता है। हम स्टेन (लाल), TFP (नीला), और अनुसंधान के लिए परिणाम बताते हैं lme4 (नारंगी)। हम स्टेन और टीएफपी के परिणामों को छायांकित करते हैं, इसलिए जब दोनों सहमत होते हैं तो बैंगनी रंग देखने की उम्मीद करते हैं। सादगी के लिए हम R से परिणामों को छायांकित नहीं करते हैं। प्रत्येक सबप्लॉट एक एकल काउंटी का प्रतिनिधित्व करता है और रेखापुंज स्कैन क्रम में अवरोही आवृत्ति में क्रमबद्ध होता है (यानी, बाएं से दाएं फिर ऊपर से नीचे)।
nrows = 17
ncols = 5
fig, ax = plt.subplots(nrows, ncols, figsize=(18, 21), sharey=True, sharex=True)
with warnings.catch_warnings():
warnings.simplefilter('ignore')
ii = -1
for r in range(nrows):
for c in range(ncols):
ii += 1
idx = county_freq[ii, 0]
sns.kdeplot(
posterior_random_weights_final_[:, idx] * log_county_uranium_ppm[idx],
color='blue',
alpha=.3,
shade=True,
label='TFP',
ax=ax[r][c])
sns.kdeplot(
posterior_random_weights_final_stan[:, idx] *
log_county_uranium_ppm[idx],
color='red',
alpha=.3,
shade=True,
label='Stan/rstanarm',
ax=ax[r][c])
sns.kdeplot(
posterior_random_weights_lme4_final_[:, idx] *
log_county_uranium_ppm[idx],
color='#F4B400',
alpha=.7,
shade=False,
label='R/lme4',
ax=ax[r][c])
ax[r][c].vlines(
posterior_random_weights_lme4[idx] * log_county_uranium_ppm[idx],
0,
5,
color='#F4B400',
linestyle='--')
ax[r][c].set_title(county_name[idx] + ' ({})'.format(idx), y=.7)
ax[r][c].set_ylim(0, 5)
ax[r][c].set_xlim(-1., 1.)
ax[r][c].get_yaxis().set_visible(False)
if ii == 2:
ax[r][c].legend(bbox_to_anchor=(1.4, 1.7), fontsize=20, ncol=3)
else:
ax[r][c].legend_.remove()
fig.subplots_adjust(wspace=0.03, hspace=0.1)
8 निष्कर्ष
इस कोलाब में हम रेडॉन डेटासेट के लिए एक रेखीय मिश्रित-प्रभाव प्रतिगमन मॉडल फिट करते हैं। हमने तीन अलग-अलग सॉफ़्टवेयर पैकेज आज़माए: R, Stan, और TensorFlow Probability। हमने तीन अलग-अलग सॉफ्टवेयर पैकेजों द्वारा गणना के अनुसार 85 पश्च वितरण की साजिश रचकर निष्कर्ष निकाला।
परिशिष्ट ए: वैकल्पिक रेडॉन एचएलएम (रैंडम इंटरसेप्ट जोड़ें)
इस खंड में हम एक वैकल्पिक एचएलएम का वर्णन करते हैं जिसमें प्रत्येक काउंटी से जुड़ा एक यादृच्छिक अवरोधन भी होता है।
\[\begin{align*} \text{for } & c=1\ldots \text{NumCounties}:\\ & \beta_c \sim \text{MultivariateNormal}\left(\text{loc}=\left[ \begin{array}{c} 0 \\ 0 \end{array}\right] , \text{scale}=\left[\begin{array}{cc} \sigma_{11} & 0 \\ \sigma_{12} & \sigma_{22} \end{array}\right] \right) \\ \text{for } & i=1\ldots \text{NumSamples}:\\ & c_i := \text{County}_i \\ &\eta_i = \underbrace{\omega_0 + \omega_1\text{Floor}_i \vphantom{\log( \text{CountyUraniumPPM}_{c_i}))} }_{\text{fixed effects} } + \underbrace{\beta_{c_i,0} + \beta_{c_i,1}\log( \text{CountyUraniumPPM}_{c_i}))}_{\text{random effects} } \\ &\log(\text{Radon}_i) \sim \text{Normal}(\text{loc}=\eta_i , \text{scale}=\sigma) \end{align*}\]
आर के दशक में lme4 "अंकन टिल्ड", इस मॉडल के बराबर है:
log_radon ~ 1 + floor + (1 + log_county_uranium_ppm | county)
परिशिष्ट बी: सामान्यीकृत रैखिक मिश्रित-प्रभाव मॉडल
इस खंड में हम मुख्य निकाय में उपयोग किए जाने वाले की तुलना में पदानुक्रमित रैखिक मॉडल का अधिक सामान्य लक्षण वर्णन देते हैं। यह और अधिक सामान्य मॉडल एक के रूप में जाना जाता है सामान्यीकृत रैखिक मिश्रित प्रभाव मॉडल (GLMM)।
GLMMs का सामान्यीकरण कर रहे हैं सामान्यीकृत रैखिक मॉडल (GLMs)। जीएलएमएम अनुमानित रैखिक प्रतिक्रिया में नमूना विशिष्ट यादृच्छिक शोर को शामिल करके जीएलएम का विस्तार करते हैं। यह आंशिक रूप से उपयोगी है क्योंकि यह दुर्लभ रूप से देखी जाने वाली सुविधाओं को अधिक सामान्य रूप से देखी जाने वाली सुविधाओं के साथ जानकारी साझा करने की अनुमति देता है।
एक जनरेटिव प्रक्रिया के रूप में, एक सामान्यीकृत रैखिक मिश्रित-प्रभाव मॉडल (GLMM) की विशेषता है:
\प्रारंभ{संरेखण} \पाठ{के लिए} और r = 1\ldots R: \hspace{2.45cm}\text \ और \ begin {गठबंधन} और \ ईटा मैं \ underbrace {\ vphantom {\ योग {r = 1} ^ R} = एक्स मैं ^ \ शीर्ष \ ओमेगा} \ text {तय प्रभाव} + \ underbrace {\ योग {r = 1} ^ आर जेड {r, मैं} ^ \ शीर्ष \ beta_ {r, सी आर (i)}} \ text {यादृच्छिक प्रभाव} \ और Y_i | x मैं, \ ओमेगा, {z {r, मैं}, \ बीटा आर} {r = 1} ^ आर \ सिम \ text {वितरण} (\ text {मतलब} = छ ^ {- 1} (\ eta_i )) \ अंत {गठबंधन} \ अंत {संरेखण}
कहाँ पे:
\प्रारंभ{संरेखण} R &= \text{यादृच्छिक-प्रभाव समूहों की संख्या}\ |C_r| और = \ text {समूह के लिए श्रेणियों की संख्या \(r\)\ N & = \ text {प्रशिक्षण नमूनों की संख्या} \ x_i, \ ओमेगा और \ में \ mathbb {R} ^ {D_0} \ D_0 & = \ text {} निश्चित प्रभाव की संख्या} \ सी आर (i) और = \ text {श्रेणी (समूह के तहत \(r\)का) \(i\)वें नमूना} \ z {r, मैं} और \ में \ mathbb {R} ^ { D_r} \ डी आर एंड = \ text {समूह के साथ जुड़े यादृच्छिक प्रभाव की संख्या \(r\) ^ {D_r \ बार D_r}} \ सिग्मा {r} और \ {S \ में \ mathbb {R में \}: एस \ succ 0 }\ \eta_i\mapsto g^{-1}(\eta_i) &= \mu_i, \text{उलटा लिंक फ़ंक्शन}\ \text{वितरण} &=\text{कुछ वितरण केवल इसके माध्य से पैरामीटर करने योग्य} \ अंत {संरेखित करें}
शब्दों में, यह कहा गया है कि प्रत्येक समूह के हर वर्ग एक आईआईडी MVN, साथ जुड़ा हुआ है \(\beta_{rc}\)। हालांकि \(\beta_{rc}\) ड्रॉ हमेशा वे केवल indentically एक समूह के लिए वितरित कर रहे हैं स्वतंत्र हैं, \(r\); नोटिस वहाँ वास्तव में एक है \(\Sigma_r\) प्रत्येक के लिए \(r\in\{1,\ldots,R\}\)।
Affinely एक नमूना के समूह के फ़ीचर, के साथ संयुक्त जब \(z_{r,i}\), परिणाम पर नमूना-विशिष्ट शोर है \(i\)वें रैखिक प्रतिक्रिया (जो अन्यथा है की भविष्यवाणी की \(x_i^\top\omega\))।
जब हम अनुमान \(\{\Sigma_r:r\in\{1,\ldots,R\}\}\) हम अनिवार्य रूप से शोर की मात्रा एक यादृच्छिक प्रभाव समूह में किया जाता है जो अन्यथा में संकेत मौजूद दबनी हैं का अनुमान लगा रहे \(x_i^\top\omega\)।
वहाँ के लिए विकल्प की एक किस्म है \(\text{Distribution}\) और उलटा लिंक समारोह , \(g^{-1}\)। आम विकल्प हैं:
- \(Y_i\sim\text{Normal}(\text{mean}=\eta_i, \text{scale}=\sigma)\),
- \(Y_i\sim\text{Binomial}(\text{mean}=n_i \cdot \text{sigmoid}(\eta_i), \text{total_count}=n_i)\), और,
- \(Y_i\sim\text{Poisson}(\text{mean}=\exp(\eta_i))\)।
अधिक संभावनाएं के लिए, देखें tfp.glm मॉड्यूल।