{TF Probability, R, Stan} में लीनियर मिक्स्ड-इफेक्ट रिग्रेशन

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 मॉड्यूल।