ヘルプKaggleにTensorFlowグレートバリアリーフを保護チャレンジに参加

{TF確率、R、スタン}における線形混合効果回帰

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

1はじめに

このコラボでは、線形混合効果回帰モデルを人気のあるおもちゃのデータセットに適合させます。私たちは、Rの使用、このフィット感を3回行いますlme4 、スタンの混合効果パッケージ、およびTensorFlow確率(TFP)プリミティブを。 3つすべてがほぼ同じ適合パラメーターと事後分布を与えることを示すことによって結論を下します。

私たちの主な結論は、TFPがHLM-のようなモデルをフィットするために必要な一般的な部分を持っており、それが他のソフトウェアパッケージと一致している結果、すなわち、生成することであるlme4rstanarm 。このコラボは、比較したパッケージの計算効率を正確に反映したものではありません。

%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階層線形モデル

R、スタン、とTFPの間に私達の比較のために、私たちは収まる階層線形モデルに(HLM)のラドンデータセットに普及したゲルマンらによるベイズデータ分析。 al。 (559ページ、第2版、250ページ、第3版)。

次の生成モデルを想定しています。

\[\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*}\]

Rさんにlme4 「チルダ表記」、このモデルは同等です:

log_radon ~ 1 + floor + (0 + log_uranium_ppm | county)

私たちは、のためにMLEを見つけるでしょう \(\omega, \sigma_C, \sigma_N\) の(証拠に条件付け)事後分布使っ \(\{\beta_c\}_{c=1}^\text{NumCounties}\)。

本質的には同じモデルのためではなく、ランダム切片、参照付録Aを

HLMsのより一般的な仕様については、付録Bを

3データの変更

このセクションでは、得られるradonデータセットを、それが私たちの想定モデルに準拠させるためにいくつかの最低限の前処理を行います。

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.

png

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)

png

結論:

  • 85の郡の長い尾があります。 (GLMMでよく発生します。)
  • 実際、 \(\log(\text{Radon})\) 拘束されません。 (したがって、線形回帰は理にかなっているかもしれません。)
  • 読みがほとんどで作られてい \(0\)番目の階。なし読み取りは、床の上に作られた \(1\)。 (したがって、固定効果には2つの重みしかありません。)

4 HLM In R

このセクションでは、Rの使用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

png

write.csv(as.data.frame(ranef(radon.model, condVar = TRUE)), '/tmp/radon/lme4_fit.csv')

5 HLM In Stan

このセクションでは、使用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))

png

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")

png

# 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)

このスタンの例では、 1が直接確率モデルを指定することにより、すなわち、近いTFPのスタイルでLMERを実装する方法を示しています。

6 HLM InTF確率

このセクションでは、低レベル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プリミティブを使用します。これを行うには、2つのTFP分布を生成する2つの関数を指定します。

  • make_weights_prior :A多変量正規(乗算される、ランダムな重みのための従来 \(\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)の確率的近似バージョンを使用します。基本的な考え方は、後方からのサンプルを使用して、予想される関節の対数密度(Eステップ)を概算することです。次に、この計算を最大化するパラメーターを見つけます(Mステップ)。もう少し具体的には、不動点反復は次の式で与えられます。

\[\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\) 可能パラメータ化します。

より完全な説明については、バーナードDelyon、マルク・Lavielle、エリック、Moulines(アン。Statist。、1999)によるEMアルゴリズムの確率的近似バージョンの収束を

Eステップを計算するには、後方からサンプリングする必要があります。後部のサンプリングは簡単ではないため、ハミルトニアンモンテカルロ(HMC)を使用します。 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)

これで、HMC遷移カーネルを作成してEステップのセットアップを完了します。

ノート:

  • 私たちは、使用state_stop_gradient=True MCMCから描画してbackproppingからM-ステップを防ぐために。 (リコールは、私たちはE-ステップは意図的に以前の最もよく知られている推定でパラメータ化されているのでを通じてバックプロパゲーションする必要はありません。)

  • 私たちは、使用tf.placeholder我々は最終的に私たちのTFグラフを実行するとき、我々は次の反復のチェーンの値として、前の反復のランダムMCMCサンプルを供給することができるように。

  • 私たちは、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)

ここで、Mステップを設定します。これは基本的に、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実行

このセクションでは、SAEMTFグラフを実行します。ここでの主なトリックは、HMCカーネルからの最後の描画を次の反復にフィードすることです。これが私たちの使用によって達成されるfeed_dictsess.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);

png

このボックスおよびウィスカー図から、我々は、郡レベルの分散ことを観察 \(\log(\text{UraniumPPM})\) 郡等のランダム効果増加が少ないデータセットに示されています。直感的には、これは理にかなっています。証拠が少ない場合は、特定の郡の影響について確信が持てないはずです。

7サイドバイサイド比較

ここで、3つの手順すべての結果を比較します。これを行うために、StanとTFPによって生成された事後サンプルの非パラメーター推定値を計算します。また、Rのによって生成推定parameteric(概算)と比較しますlme4パッケージ。

次のプロットは、ミネソタ州の各郡の各重みの事後分布を示しています。私たちは、スタン(赤)、TFP(青)、およびRのための結果を示してlme4 (オレンジ)。スタンとTFPの結果をシェーディングするため、2つが一致すると紫色になると予想されます。簡単にするために、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)

png

8結論

このコラボでは、線形混合効果回帰モデルをラドンデータセットに適合させます。 R、Stan、TensorFlowProbabilityの3つの異なるソフトウェアパッケージを試しました。最後に、3つの異なるソフトウェアパッケージによって計算された85の事後分布をプロットしました。

付録A:代替ラドンHLM(ランダムインターセプトの追加)

このセクションでは、各郡にランダムな切片が関連付けられている代替HLMについて説明します。

\[\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*}\]

Rさんにlme4 「チルダ表記」、このモデルは同等です:

log_radon ~ 1 + floor + (1 + log_county_uranium_ppm | county)

付録B:一般化線形混合効果モデル

このセクションでは、本体で使用されているものよりも、階層線形モデルのより一般的な特性を示します。このより一般的なモデルとして知られている一般化線形混合効果モデル(GLMM)。

GLMMsはの一般化したものである一般化線形モデル(GLMS)。 GLMMは、サンプル固有のランダムノイズを予測線形応答に組み込むことにより、GLMを拡張します。これは、めったに見られない機能がより一般的に見られる機能と情報を共有できるようにするため、部分的には便利です。

生成プロセスとして、一般化線形混合効果モデル(GLMM)は次の特徴があります。

\ begin {align} \ text {for}&r = 1 \ ldots R:\ hspace {2.45cm} \ text {#各ランダム効果グループ} \&\ begin {aligned} \ text {for}&c = 1 \ ldots | CのR |:\ HSPACE {1.3センチメートル} \ {グループの各カテゴリ( "レベル")のため#テキスト \(r\)} \&\開始{整列} \ベータ{RC}&\シム\テキスト{MultivariateNormal }(\ text {loc} = 0_ {D_r}、\ text {scale} = \ Sigma_r ^ {1/2})\ end {aligned} \ end {aligned} \\ \ text {for}&i = 1 \ ldots N:\ HSPACE {2.45センチメートル} \ {各サンプルのため#}テキスト\&\ \ ETA iが\ underbrace {\ vphantom {\和{R = 1} ^ Rを} = {整列}&開始X I ^ \上部\オメガ} \テキスト{固定効果} + \ underbrace {\和{R = 1} ^ R zの{R、I} ^ \上部\ beta_ {R、C R(I)}} \テキスト{ランダム効果} \ &Y_I | x iは、\オメガ、{Z {R、I}、\ベータR} {R = 1} ^ R \ SIM \テキスト{分布}(\テキスト{平均} = G ^ { - 1}(\ eta_i ))\ end {aligned} \ end {align}

どこ:

\ begin {align} R&= \ text {ランダム効果グループの数} \ | C_r | &= \テキスト{グループのカテゴリの数 \(r\)} \ N&= \テキスト{訓練サンプルの数} \ X_I、\オメガ&\ \ mathbb {R} ^ {D_0} \ D_0&= \テキスト{で固定効果の数}(グループ下\ C R(I)&= \テキスト{カテゴリ \(r\)の) \(i\)番目のサンプル} \ Z {R、I}&\で\ mathbb {R} ^ { D_R} \ R&D = \テキスト{グループに関連付けられたランダム効果の数 \(r\) ^ {D_R \回D_R} {S \で\ mathbb {Rにおける} \ \シグマ{R}&\}:S \ succ 0} \ \ eta_i \ mapsto g ^ {-1}(\ eta_i)&= \ mu_i、\ text {逆リンク関数} \ \ text {分布}&= \ text {平均によってのみパラメータ化可能な分布} \ end {align}

言葉では、これは、各グループのすべてのカテゴリがIID MVN、関連付けられていることを述べている \(\beta_{rc}\)。けれども \(\beta_{rc}\) 描く彼らだけindenticallyグループのために配布され、常に独立している \(r\)。通知は正確に一つあり \(\Sigma_r\) ごとに、 \(r\in\{1,\ldots,R\}\)。

アフィンサンプルのグループの特徴と組み合わされたとき \(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モジュールを。