ホーム » 「TensorFlow Probability」タグがついた投稿

タグアーカイブ: TensorFlow Probability

TensorFlow Probability : Tutorials : 確率的 PCA

TensorFlow Probability : Tutorials : 確率的 PCA (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 11/19/2017

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです:

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : 確率的 PCA

確率的主成分分析 (PCA) は低次元潜在空間を通してデータを解析する次元削減テクニックです (Tipping and Bishop 1999)。それはデータに欠落値があるときやマルチ次元スケーリングのためにしばしば使用されます。

 

インポート

!pip install -q tensorflow-probability
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

import tensorflow as tf
import tensorflow_probability as tfp
from tensorflow_probability import edward2 as ed
import warnings

plt.style.use("ggplot")
warnings.filterwarnings('ignore')

 

モデル

$N$ データポイントのデータ集合 $\mathbf{X} = \{\mathbf{x}_n\}$ を考えます、そこでは各データポイントは $D$-次元で、$\mathbf{x}_n \in \mathbb{R}^D$ です。各 $\mathbf{x}_n$ をより低次元 $K < D$ である潜在変数 $\mathbf{z}_n \in \mathbb{R}^K$ で表わすことを狙います。主軸の集合 $\mathbf{W}$ は潜在変数をデータに関係づけます。

特に、各潜在変数は正規分布に従うものと仮定します :
$$
\begin{equation*}
\mathbf{z}_n \sim N(\mathbf{0}, \mathbf{I}).
\end{equation*}
$$

対応するデータポイントは射影を通して生成されます :
$$
\begin{equation*}
\mathbf{x}_n \mid \mathbf{z}_n
\sim N(\mathbf{W}\mathbf{z}_n, \sigma^2\mathbf{I}),
\end{equation*}
$$

ここで行列 $\mathbf{W}\in\mathbb{R}^{D\times K}$ は主軸として知られています。確率的 PCA では、典型的には主軸 $\mathbf{W}$ とノイズ項 $\sigma^2$ を推定することに興味があります。

確率的 PCA は古典的 PCA を一般化します。潜在変数を周辺化すると、各データポイントの分布は :
$$
\begin{equation*}
\mathbf{x}_n \sim N(\mathbf{0}, \mathbf{W}\mathbf{W}^\top + \sigma^2\mathbf{I}).
\end{equation*}
$$

ノイズの共分散が無限に小さくなるとき、$\sigma^2 \to 0$、古典的 PCA は確率的 PCA の特定のケースです。

下でモデルをセットアップします。私達の分析では、$\sigma$ は既知と仮定し、モデルパラメータとして $\mathbf{W}$ を点推定する代わりに、主軸に渡る分布を推論するためにそれに渡る事前分布を置きます。

def probabilistic_pca(data_dim, latent_dim, num_datapoints, stddv_datapoints): # (unmodeled) data
  w = ed.Normal(loc=tf.zeros([data_dim, latent_dim]),
                scale=2.0 * tf.ones([data_dim, latent_dim]),
                name="w")  # parameter
  z = ed.Normal(loc=tf.zeros([latent_dim, num_datapoints]),
                scale=tf.ones([latent_dim, num_datapoints]), 
                name="z")  # parameter
  x = ed.Normal(loc=tf.matmul(w, z),
                scale=stddv_datapoints * tf.ones([data_dim, num_datapoints]),
                name="x")  # (modeled) data
  return x, (w, z)

log_joint = ed.make_log_joint_fn(probabilistic_pca)

 

データ

Edward2 モデルをデータを生成するために使用できます。

num_datapoints = 5000
data_dim = 2
latent_dim = 1
stddv_datapoints = 0.5

model = probabilistic_pca(data_dim=data_dim,
                          latent_dim=latent_dim,
                          num_datapoints=num_datapoints,
                          stddv_datapoints=stddv_datapoints)

with tf.Session() as sess:
  x_train, (actual_w, actual_z) = sess.run(model)

print("Principal axes:")
print(actual_w)
Principal axes:
[[-0.10648464]
 [ 1.3358579 ]]

データセットを可視化します。

plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1)
plt.axis([-20, 20, -20, 20])
plt.title("Data set")
plt.show()

 

最大事後確率推論

最初に事後確率密度を最大化する潜在変数の点推定を探索します。これは最大事後確率 (MAP) 推論として知られ、事後密度 $p(\mathbf{W}, \mathbf{Z} \mid \mathbf{X}) \propto p(\mathbf{W}, \mathbf{Z}, \mathbf{X})$ を最大化する$\mathbf{W}$ と $\mathbf{Z}$ の値を計算することで成されます。

tf.reset_default_graph()

w = tf.Variable(np.ones([data_dim, latent_dim]), dtype=tf.float32)
z = tf.Variable(np.ones([latent_dim, num_datapoints]), dtype=tf.float32)


def target(w, z):
  """Unnormalized target density as a function of the parameters."""
  return log_joint(data_dim=data_dim,
                   latent_dim=latent_dim,
                   num_datapoints=num_datapoints,
                   stddv_datapoints=stddv_datapoints,
                   w=w, z=z, x=x_train)


energy = -target(w, z)

optimizer = tf.train.AdamOptimizer(learning_rate=0.05)
train = optimizer.minimize(energy)

init = tf.global_variables_initializer()

t = []

num_epochs = 200

with tf.Session() as sess:
  sess.run(init)

  for i in range(num_epochs):
    sess.run(train)
    if i % 5 == 0:
      cE, cw, cz = sess.run([energy, w, z])
      t.append(cE)

  w_inferred_map = sess.run(w)
  z_inferred_map = sess.run(z)
x = range(1, num_epochs, 5)
plt.plot(x, t)

$\mathbf{W}$ と $\mathbf{Z}$ のために推論された値のためのデータをサンプリングするために Edward2 モデルを使用し、そして (その上で条件付けした) 実際のデータセットと比較することができます。

print("MAP-estimated axes:")
print(w_inferred_map)

def replace_latents(w=actual_w, z=actual_z):

  def interceptor(rv_constructor, *rv_args, **rv_kwargs):
    """Replaces the priors with actual values to generate samples from."""
    name = rv_kwargs.pop("name")
    if name == "w":
      rv_kwargs["value"] = w
    elif name == "z":
      rv_kwargs["value"] = z
    return rv_constructor(*rv_args, **rv_kwargs)

  return interceptor

with ed.interception(replace_latents(w_inferred_map, z_inferred_map)):
  generate = probabilistic_pca(
      data_dim=data_dim, latent_dim=latent_dim,
      num_datapoints=num_datapoints, stddv_datapoints=stddv_datapoints)

with tf.Session() as sess:
  x_generated, _ = sess.run(generate)

plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1, label='Actual data')
plt.scatter(x_generated[0, :], x_generated[1, :], color='red', alpha=0.1, label='Simulated data (MAP)')
plt.legend()
plt.axis([-20, 20, -20, 20])
plt.show()
MAP-estimated axes:
[[-0.20085734]
 [ 2.7221382 ]]

 

変分推論

MAP は事後分布のモード (or モードの一つ) を見つけるために使用できますが、それについて他の洞察は一切提供しません。次に、変分推論を使用します、そこでは事後分布 $p(\mathbf{W}, \mathbf{Z} \mid \mathbf{X})$ は $\boldsymbol{\lambda}$ によりパラメータ化された変分分布 $q(\mathbf{W}, \mathbf{Z})$ を使用して近似されます。目的は、q と事後分布 $\mathrm{KL}(q(\mathbf{W}, \mathbf{Z}) \mid\mid p(\mathbf{W}, \mathbf{Z} \mid \mathbf{X}))$ の間の KL ダイバージェンスを最小化する、あるいは同値にエビデンス下界 $\mathbb{E}_{q(\mathbf{W},\mathbf{Z};\boldsymbol{\lambda})}\left[ \log p(\mathbf{W},\mathbf{Z},\mathbf{X}) – \log q(\mathbf{W},\mathbf{Z}; \boldsymbol{\lambda}) \right]$ を最大化する、変分パラメータ $\boldsymbol{\lambda}$ を見つけることです。

tf.reset_default_graph()

def variational_model(qw_mean, qw_stddv, qz_mean, qz_stddv):
  qw = ed.Normal(loc=qw_mean, scale=qw_stddv, name="qw")
  qz = ed.Normal(loc=qz_mean, scale=qz_stddv, name="qz")
  return qw, qz


log_q = ed.make_log_joint_fn(variational_model)

def target_q(qw, qz):
  return log_q(qw_mean=qw_mean, qw_stddv=qw_stddv,
               qz_mean=qz_mean, qz_stddv=qz_stddv,
               qw=qw, qz=qz)


qw_mean = tf.Variable(np.ones([data_dim, latent_dim]), dtype=tf.float32)
qz_mean = tf.Variable(np.ones([latent_dim, num_datapoints]), dtype=tf.float32)
qw_stddv = tf.nn.softplus(tf.Variable(-4 * np.ones([data_dim, latent_dim]), dtype=tf.float32))
qz_stddv = tf.nn.softplus(tf.Variable(-4 * np.ones([latent_dim, num_datapoints]), dtype=tf.float32))

qw, qz = variational_model(qw_mean=qw_mean, qw_stddv=qw_stddv,
                           qz_mean=qz_mean, qz_stddv=qz_stddv)

energy = target(qw, qz)
entropy = -target_q(qw, qz)

elbo = energy + entropy


optimizer = tf.train.AdamOptimizer(learning_rate = 0.05)
train = optimizer.minimize(-elbo)

init = tf.global_variables_initializer()

t = []

num_epochs = 100

with tf.Session() as sess:
  sess.run(init)

  for i in range(num_epochs):
    sess.run(train)
    if i % 5 == 0:
      t.append(sess.run([elbo]))

  w_mean_inferred = sess.run(qw_mean)
  w_stddv_inferred = sess.run(qw_stddv)
  z_mean_inferred = sess.run(qz_mean)
  z_stddv_inferred = sess.run(qz_stddv)
print("Inferred axes:")
print(w_mean_inferred)
print("Standard Deviation:")
print(w_stddv_inferred)

plt.plot(range(1, num_epochs, 5), t)
plt.show()

with ed.interception(replace_latents(w_mean_inferred, z_mean_inferred)):
  generate = probabilistic_pca(
      data_dim=data_dim, latent_dim=latent_dim,
      num_datapoints=num_datapoints, stddv_datapoints=stddv_datapoints)

with tf.Session() as sess:
  x_generated, _ = sess.run(generate)

plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1, label='Actual data')
plt.scatter(x_generated[0, :], x_generated[1, :], color='red', alpha=0.1, label='Simulated data (VI)')
plt.legend()
plt.axis([-20, 20, -20, 20])
plt.show()

 

References

  • [1]: Michael E. Tipping and Christopher M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3): 611-622, 1999.
 

以上






TensorFlow Probability : Tutorials : 線形混合効果モデル

TensorFlow Probability : Tutorials : 線形混合効果モデル (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 11/03/2017

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです:

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : 線形混合効果モデル

線形混合効果モデルは構造化された線形関係をモデル化するための単純なアプローチです (Harville, 1997; Laird and Ware, 1982)。各データポイントは — グループにカテゴリー分けされた — 様々なタイプと実数値の出力から成ります。線形混合効果モデルは階層モデルです: それは任意の個々のデータポイントについて推論を改良するためにグループに渡り統計的強度を共有します。

このチュートリアルでは、TensorFlow Probability で現実世界の例で線形混合効果モデルを示します。Edward2 (tfp.edward2) とマルコフ連鎖モンテカルロ (tfp.mcmc) モジュールを使用します。

%matplotlib inline
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import csv
import IPython
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
import tensorflow as tf
import tensorflow_probability as tfp
import warnings

from tensorflow_probability import edward2 as ed

plt.style.use('ggplot')

 

データ

R のポピュラーな lme4 パッケージ からの InstEval データセットを使用します (Bates et al., 2015)。それは課程とそれらの評価格付けのデータセットです。各課程は生徒、講師そして学部のようなメタデータを含みそして関心のある応答変数は評価格付けです。

def load_insteval():
  """Loads the InstEval data set.

  It contains 73,421 university lecture evaluations by students at ETH
  Zurich with a total of 2,972 students, 2,160 professors and
  lecturers, and several student, lecture, and lecturer attributes.
  Implementation is built from the `observations` Python package.

  Returns:
    Tuple of np.darray `x_train` with 73,421 rows and 7 columns and
    dictionary `metadata` of column headers (feature names).
  """
  url = ('https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/'
         'lme4/InstEval.csv')
  with requests.Session() as s:
    download = s.get(url)
    f = download.content.decode().splitlines()

  iterator = csv.reader(f)
  columns = next(iterator)[1:]
  x_train = np.array([row[1:] for row in iterator], dtype=np.int)
  metadata = {'columns': columns}
  return x_train, metadata

データセットをロードして前処理します。データの 20% を取り置きますので、未見のデータポイント上で適合されたモデルを評価できます。下で最初の数行を可視化します。

data, metadata = load_insteval()
data = pd.DataFrame(data, columns=metadata['columns'])
data = data.rename(columns={'s': 'students',
                            'd': 'instructors',
                            'dept': 'departments',
                            'y': 'ratings'})
data['students'] -= 1  # start index by 0
# Remap categories to start from 0 and end at max(category).
data['instructors'] = data['instructors'].astype('category').cat.codes
data['departments'] = data['departments'].astype('category').cat.codes

train = data.sample(frac=0.8)
test = data.drop(train.index)

train.head()

生徒
(students)
講師
(instructors)
studage lectage service 学部
(departments)
格付け
(ratings)
40736 1651 17 4 1 1 12 5
23953 967 226 6 2 0 11 4
4252 159 897 6 5 1 1 5
64541 2616 561 8 6 1 7 1
59310 2404 838 4 2 1 1 3

入力の特徴辞書と格付けに相当するラベル出力の観点からデータセットをセットアップします。各特徴は整数としてエンコードされて各ラベル (評価格付け) は浮動小数点数としてエンコードされます。


get_value = lambda dataframe, key, dtype: dataframe[key].values.astype(dtype)
features_train = {
    k: get_value(train, key=k, dtype=np.int32)
    for k in ['students', 'instructors', 'departments', 'service']}
labels_train = get_value(train, key='ratings', dtype=np.float32)

features_test = {k: get_value(test, key=k, dtype=np.int32)
                 for k in ['students', 'instructors', 'departments', 'service']}
labels_test = get_value(test, key='ratings', dtype=np.float32)
num_students = max(features_train['students']) + 1
num_instructors = max(features_train['instructors']) + 1
num_departments = max(features_train['departments']) + 1
num_observations = train.shape[0]

print("Number of students:", num_students)
print("Number of instructors:", num_instructors)
print("Number of departments:", num_departments)
print("Number of observations:", num_observations)
Number of students: 2972
Number of instructors: 1128
Number of departments: 14
Number of observations: 58737

 

モデル

典型的な線形モデルは独立性を仮定します、そこではデータポイントの任意のペアは定数線形関係を持ちます。InstEval データセットでは、観測はグループ内で発生してその各々は様々な傾きと切片を持つかもしれません。階層線形モデルまたはマルチレベル線形モデルとしても知られる、線形混合効果モデルはこの現象を捕捉します (Gelman & Hill, 2006)。

現象の例は次を含みます :

  • 生徒 (= Students): 生徒からの観測は独立ではありません: 何人かの生徒はシステマチックに低い (or 高い) 講義格付けを与えるかもしれません。
  • 講師 (= Instructors): 講師からの観測は独立ではありません: 良い講師は一般に良い格付けを持ちそして悪い講師は一般に悪い格付けを持つことを期待します。
  • 学部 (= Departments): 学部からの観測は独立ではありません: ある学部は一般に無味乾燥なマテリアルやより厳しい等級付け (= grading) を持ち、そしてそのため他よりも低く格付けされるかもしれません。

これを捕捉するために、$N\times D$ 特徴 $\mathbf{X}$ と $N$ ラベル $\mathbf{y}$ のデータセットについて、線形回帰は次のモデルを仮定することを思い出してください :

$$
\begin{equation*}
\mathbf{y} = \mathbf{X}\beta + \alpha + \epsilon,
\end{equation*}
$$

ここで傾きベクトル $\beta\in\mathbb{R}^D$、切片 intercept $\alpha\in\mathbb{R}$ そしてランダム・ノイズ $\epsilon\sim\text{Normal}(\mathbf{0}, \mathbf{I})$ があります。

$\beta$ と $\alpha$ は「固定効果」であると言います : それらはデータポイント $(x, y)$ の母集団に渡り定数として保持されるは効果です。尤度としての等式の同値な定式化は $\mathbf{y} \sim \text{Normal}(\mathbf{X}\beta + \alpha, \mathbf{I})$ です。この尤度は、データに適合する $\beta$ と $\alpha$ の点推定を見つけるために推論の間に最大化されます。

線形混合効果モデルは線形回帰を次として拡張します :
$$
\begin{align*}
\eta &\sim \text{Normal}(\mathbf{0}, \sigma^2 \mathbf{I}), \\
\mathbf{y} &= \mathbf{X}\beta + \mathbf{Z}\eta + \alpha + \epsilon.
\end{align*}
$$

ここで依然として傾きベクトル $\beta\in\mathbb{R}^P$、切片 $\alpha\in\mathbb{R}$ そしてランダム・ノイズ $\epsilon\sim\text{Normal}(\mathbf{0}, \mathbf{I})$ があります。加えて、項 $\mathbf{Z}\eta$ があります、そこでは $\mathbf{Z}$ は特徴行列で $\eta\in\mathbb{R}^Q$ はランダム傾きのベクトルです ; $\eta$ は分散成分パラメータ $\sigma^2$ を持ち正規分布に従います。$\mathbf{Z}$ は新しい$N\times P$ 行列 $\mathbf{X}$ と $N\times Q$ 行列 $\mathbf{Z}$ の観点で元の $N\times D$ 特徴行列を分割することにより形成されます、ここで $P + Q=D$ です : この分割は固定効果 $\beta$ と潜在変数 $\eta$ をそれぞれ使用して特徴を別々にモデル化することを可能にします。

潜在変数 $\eta$ を「ランダム効果 (変量効果)」と言います: それらは母数に渡り変化する効果です (それらは部分母集団に渡り定数かもしれませんが)。特に、ランダム効果 $\eta$ は平均 0 を持ちますので、データラベルの平均は $\mathbf{X}\beta + \alpha$ で捕捉されます。ランダム効果成分はデータの変分 $\mathbf{Z}\eta$ を捕捉します: 例えば、「講師 #54 は平均よりも 1.4 ポイント高く格付けされます」

このチュートリアルでは以下の効果を仮定します :

  • 固定効果: service (サービス)。service は課程が講師の主要な学部に属するか否かに対応する二値共変量です。どれほど追加データを集めようとも、それは値 $0$ と $1$ を取るだけです。
  • ランダム効果: students (生徒)、instructors (講師) と departments (学部)。課程評価格付けの母集団からより多くの観測が与えられたとき、新しい生徒、先生 (= teachers) または学部を見ているかもしれません。

R の lme4 パッケージ (Bates et al., 2015) のシンタックスでは、モデルは次のように要約できます :

ratings ~ service + (1|students) + (1|instructors) + (1|departments) + 1

ここで x は固定効果を表し、(1|x) は x のためのランダム効果を表し、そして 1 は切片項を表します。

このモデルを下で Edward プログラムとして実装します。

def linear_mixed_effects_model(features):
  # Set up fixed effects and other parameters.
  intercept = tf.get_variable("intercept", [])  # alpha in eq
  effect_service = tf.get_variable("effect_service", [])  # beta in eq
  stddev_students = tf.exp(
      tf.get_variable("stddev_unconstrained_students", []))  # sigma in eq
  stddev_instructors = tf.exp(
      tf.get_variable("stddev_unconstrained_instructors", [])) # sigma in eq
  stddev_departments = tf.exp(
      tf.get_variable("stddev_unconstrained_departments", [])) # sigma in eq

  # Set up random effects.
  effect_students = ed.MultivariateNormalDiag(
      loc=tf.zeros(num_students),
      scale_identity_multiplier=stddev_students,
      name="effect_students")
  effect_instructors = ed.MultivariateNormalDiag(
      loc=tf.zeros(num_instructors),
      scale_identity_multiplier=stddev_instructors,
      name="effect_instructors")
  effect_departments = ed.MultivariateNormalDiag(
      loc=tf.zeros(num_departments),
      scale_identity_multiplier=stddev_departments,
      name="effect_departments")

  # Set up likelihood given fixed and random effects.
  # Note we use `tf.gather` instead of matrix-multiplying a design matrix of
  # one-hot vectors. The latter is memory-intensive if there are many groups.
  ratings = ed.Normal(
      loc=(effect_service * features["service"] +
           tf.gather(effect_students, features["students"]) +
           tf.gather(effect_instructors, features["instructors"]) +
           tf.gather(effect_departments, features["departments"]) +
           intercept),
      scale=1.,
      name="ratings")
  return ratings

# Wrap model in a template. All calls to the model template will use the same
# TensorFlow variables.
model_template = tf.make_template("model", linear_mixed_effects_model)

Edward プログラムとして、その計算グラフの視点からモデルの構造も可視化できます。このグラフはプログラムの確率変数に渡るデータフローを、グラフィカルモデルの視点からそれらの関係を明確にしてエンコードします (Jordan, 2003)。

統計ツールとして、例えば、intercept と effect_service が与えられた格付けに条件付きで依存していることをより良く見るためにグラフを見るかもしれません; これはプログラムがクラス、モジュールに渡る交差参照、and/or サブルーチンで書かれている場合ソースコードから見ることはより困難かもしれません。計算ツールとして、tf.gather ops を通して ratings 変数になる潜在変数に気付くかもしれません。これは indexing Tensor がコスト高であればあるハードウェア・アクセラレータではボトルネックかもしれません; グラフの可視化はこれを容易に明白にします。

def strip_consts(graph_def, max_const_size=32):
  """Strip large constant values from graph_def."""
  strip_def = tf.GraphDef()
  for n0 in graph_def.node:
    n = strip_def.node.add()
    n.MergeFrom(n0)
    if n.op == 'Const':
      tensor = n.attr['value'].tensor
      size = len(tensor.tensor_content)
      if size > max_const_size:
        tensor.tensor_content = bytes("<stripped %d bytes>"%size, 'utf-8')
  return strip_def

def draw_graph(model, *args, **kwargs):
  """Visualize TensorFlow graph."""
  graph = tf.Graph()
  with graph.as_default():
    model(*args, **kwargs)
  graph_def = graph.as_graph_def()
  strip_def = strip_consts(graph_def, max_const_size=32)
  code = """
      <script>
        function load() {{
          document.getElementById("{id}").pbtxt = {data};
        }}
      </script>
      <link rel="import" href="https://tensorboard.appspot.com/tf-graph-basic.build.html" onload=load()>
      <div style="height:600px">
        <tf-graph-basic id="{id}"></tf-graph-basic>
      </div>
  """.format(data=repr(str(strip_def)), id='graph'+str(np.random.rand()))

  iframe = """
      <iframe seamless style="width:1200px;height:620px;border:0" srcdoc="{}"></iframe>
  """.format(code.replace('"', '"'))
  IPython.display.display(IPython.display.HTML(iframe))


draw_graph(linear_mixed_effects_model, features_train)

 

パラメータ推定

データが与えられたとき、推論のゴールはモデルの固定効果傾き $\beta$、切片 $\alpha$ 、そして分散成分パラメータ $\sigma^2$ を適合させることです。最大尤度原理はこのタスクを次のように形式化します :
$$
\max_{\beta, \alpha, \sigma}~\log p(\mathbf{y}\mid \mathbf{X}, \mathbf{Z}; \beta, \alpha, \sigma) = \max_{\beta, \alpha, \sigma}~\log \int p(\eta; \sigma) ~p(\mathbf{y}\mid \mathbf{X}, \mathbf{Z}, \eta; \beta, \alpha)~d\eta.
$$

このチュートリアルでは、この周辺密度を最大化するためにモンテカルロ EM アルゴリズムを使用します (Dempster et al., 1977; Wei and Tanner, 1990)。ランダム効果 (“E-step”) に関する条件付き尤度の期待値を計算するためにマルコフ連鎖モンテカルロを遂行し、そしてパラメータ (“M-step”) に関する期待値を最大化するために勾配降下を遂行します :

  • E-step については、ハミルトニアン・モンテカルロ (HMC) をセットアップします。それは現在の状態— student, instructor と department 効果 — を取りそして新しい状態を返します。その新しい状態を TensorFlow 変数に割り当て、これは HMC chain の状態を表します。
  • M-step については、定数までの周辺尤度の公平な (= unbiased) 推定を計算するために HMC からの事後分布サンプルを使用します。それから関心のあるパラメータに関する勾配を適用します。これは周辺尤度上の公平な確率的勾配ステップを生成します。それを Adam TensorFlow optimizer で実装して周辺の負を最小化します。
log_joint = ed.make_log_joint_fn(model_template)

def target_log_prob_fn(effect_students, effect_instructors, effect_departments):
  """Unnormalized target density as a function of states."""
  return log_joint(  # fix `features` and `ratings` to the training data
    features=features_train,
    effect_students=effect_students,
    effect_instructors=effect_instructors,
    effect_departments=effect_departments,
    ratings=labels_train)

tf.reset_default_graph()

# Set up E-step (MCMC).
effect_students = tf.get_variable(  # `trainable=False` so unaffected by M-step
    "effect_students", [num_students], trainable=False)
effect_instructors = tf.get_variable(
    "effect_instructors", [num_instructors], trainable=False)
effect_departments = tf.get_variable(
    "effect_departments", [num_departments], trainable=False)

hmc = tfp.mcmc.HamiltonianMonteCarlo(
    target_log_prob_fn=target_log_prob_fn,
    step_size=0.015,
    num_leapfrog_steps=3)

current_state = [effect_students, effect_instructors, effect_departments]
next_state, kernel_results = hmc.one_step(
      current_state=current_state,
      previous_kernel_results=hmc.bootstrap_results(current_state))

expectation_update = tf.group(
    effect_students.assign(next_state[0]),
    effect_instructors.assign(next_state[1]),
    effect_departments.assign(next_state[2]))

# Set up M-step (gradient descent).
# The following should work. However, TensorFlow raises an error about taking
# gradients through IndexedSlices tensors. This may be a TF bug. For now,
# we recompute the target's log probability at the current state.
# loss = -kernel_results.accepted_results.target_log_prob
with tf.control_dependencies([expectation_update]):
  loss = -target_log_prob_fn(effect_students,
                             effect_instructors,
                             effect_departments)
  optimizer = tf.train.AdamOptimizer(learning_rate=0.01)
  minimization_update = optimizer.minimize(loss)

warm-up ステージを遂行します、これは訓練が事後分布の確率質量内で初期化されるように多くの反復のために一つの MCMC chain を実行します。それから訓練ループを実行します。それは結合的に E と M-steps を実行して訓練の間に値を記録します。

init = tf.global_variables_initializer()

num_warmup_iters = 1000
num_iters = 1500
num_accepted = 0
effect_students_samples = np.zeros([num_iters, num_students])
effect_instructors_samples = np.zeros([num_iters, num_instructors])
effect_departments_samples = np.zeros([num_iters, num_departments])
loss_history = np.zeros([num_iters])

sess = tf.Session()
sess.run(init)

# Run warm-up stage.
for t in range(num_warmup_iters):
  _, is_accepted_val = sess.run(
      [expectation_update, kernel_results.is_accepted])
  num_accepted += is_accepted_val
  if t % 500 == 0 or t == num_warmup_iters - 1:
    print("Warm-Up Iteration: {:>3} Acceptance Rate: {:.3f}".format(
        t, num_accepted / (t + 1)))

num_accepted = 0  # reset acceptance rate counter

# Run training.
for t in range(num_iters):
  for _ in range(5):  # run 5 MCMC iterations before every joint EM update
    _ = sess.run(expectation_update)
  [
      _,
      _,
      effect_students_val,
      effect_instructors_val,
      effect_departments_val,
      is_accepted_val,
      loss_val,
  ] = sess.run([
      expectation_update,
      minimization_update,
      effect_students,
      effect_instructors,
      effect_departments,
      kernel_results.is_accepted,
      loss,
  ])
  effect_students_samples[t, :] = effect_students_val
  effect_instructors_samples[t, :] = effect_instructors_val
  effect_departments_samples[t, :] = effect_departments_val
  num_accepted += is_accepted_val
  loss_history[t] = loss_val
  if t % 500 == 0 or t == num_iters - 1:
    print("Iteration: {:>4} Acceptance Rate: {:.3f} Loss: {:.3f}".format(
        t, num_accepted / (t + 1), loss_val))
Warm-Up Iteration:   0 Acceptance Rate: 1.000
Warm-Up Iteration: 500 Acceptance Rate: 0.806
Warm-Up Iteration: 999 Acceptance Rate: 0.736
Iteration:    0 Acceptance Rate: 1.000 Loss: 103901.820
Iteration:  500 Acceptance Rate: 0.647 Loss: 96190.906
Iteration: 1000 Acceptance Rate: 0.671 Loss: 96115.094
Iteration: 1499 Acceptance Rate: 0.674 Loss: 96050.664

上で、収束閾値が検出されるまでアルゴリズムを実行しませんでした。訓練に道理があったかを確認するために、損失関数が訓練反復に渡り実際に収束する傾向があるかを検証します。

plt.plot(loss_history)
plt.ylabel(r'Loss $-\log$ $p(y\mid\mathbf{x})$')
plt.xlabel('Iteration')
plt.show()

trace プロットもまた使用します、これは特定の潜在次元に渡るマルコフ連鎖モンテカルロ・アルゴリズムの軌道を示します。下では特定の講師効果が実際にそれらの初期状態から離れて意味深く推移して状態空間を探究することを見ます。trace プロットはまた (類似の混合した挙動を持ちながら)効果が講師に渡り異なることを示しています。

for i in range(7):
  plt.plot(effect_instructors_samples[:, i])

plt.legend([i for i in range(7)], loc='lower right')
plt.ylabel('Instructor Effects')
plt.xlabel('Iteration')
plt.show()

 

批評

上で、モデルを適合させました。今はデータを使用してその適合を批評する調査をします、これはモデルを探究してより良く理解させてくれます。そのようなテクニックの一つは残差プロットです、これは各データポイントに対するモデルの予測と正解の差異をプロットします。モデルが正しければ、それらの差異は標準正規分布に従います; プロットのこのパターンからの任意の偏差はモデルの不適合を示します。

最初に訓練データが与えられたときの事後分布を持つランダム効果上の格付けに渡り事後予想分布を形成することにより残差プロットを構築します。特にモデルを forward に実行して推論された事後平均² を持つ事前ランダム効果上の依存性をインターセプトします。

def interceptor(rv_constructor, *rv_args, **rv_kwargs):
  """Replaces prior on effects with empirical posterior mean from MCMC."""
  name = rv_kwargs.pop("name")
  if name == "effect_students":
    rv_kwargs["value"] = np.mean(effect_students_samples, 0)
  elif name == "effect_instructors":
    rv_kwargs["value"] = np.mean(effect_instructors_samples, 0)
  elif name == "effect_departments":
    rv_kwargs["value"] = np.mean(effect_departments_samples, 0)
  return rv_constructor(*rv_args, **rv_kwargs)

with ed.interception(interceptor):
  ratings_posterior = model_template(features=features_test)

ratings_prediction = ratings_posterior.distribution.mean()

視覚的な調査において、残差はある程度標準正規分布に従うように見えます。けれども適合は完全ではありません; 終わりの方で正規分布よりもより大きな確率質量があり、これはモデルがその normality assumption (訳注: 正規分布に従うという仮定) をその適合を改善するかもしれないことを示します。

特に、InstEval データセットで格付けをモデル化するために正規分布を使用することは殆ど一般的ですが、データを良く見ると課程評価格付けは実際に 1 から 5 の元の値であることを明らかにします。これは、十分なデータがあれば元の分布か Categorical でさえも使用するべきであることを提案します。これが上の Edward プログラムへの 1-行の変更です; 同じ推論コードが適用可能です。

ratings_pred = sess.run(ratings_prediction)
plt.title("Residuals for Predicted Ratings on Test Set")
plt.xlim(-4, 4)
plt.ylim(0, 800)
plt.hist(ratings_pred - labels_test, 75)
plt.show()

モデルが個々の予測をどのように行なうかを調査するために、students, instructors と departments のための効果のヒストグラムを見ます。これはデータポイントの特徴ベクトルの個々の成分がどのように出力結果に影響する傾向があるかを理解させてくれます。

驚くほどのことではありませんが、各生徒は講師の評価格付け上に典型的には殆ど効果を持たないことを下で見ます。興味深いことに、講師が所属する学部は大きな効果を持つことを見ます。

[
    effect_students_mean,
    effect_instructors_mean,
    effect_departments_mean,
] = sess.run([
    tf.reduce_mean(effect_students_samples, 0),
    tf.reduce_mean(effect_instructors_samples, 0),
    tf.reduce_mean(effect_departments_samples, 0),
])
plt.title("Histogram of Student Effects")
plt.hist(effect_students_mean, 75)
plt.show()

plt.title("Histogram of Instructor Effects")
plt.hist(effect_instructors_mean, 75)
plt.show()

plt.title("Histogram of Department Effects")
plt.hist(effect_departments_mean, 75)
plt.show()

 

References

  • Douglas Bates and Martin Machler and Ben Bolker and Steve Walker. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1):1-48, 2015.
  • Arthur P. Dempster, Nan M. Laird, and Donald B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B (Methodological), 1-38, 1977.
  • Andrew Gelman and Jennifer Hill. Data analysis using regression and multilevel/hierarchical models. Cambridge University Press, 2006.
  • David A. Harville. Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association, 72(358):320-338, 1977.
  • Michael I. Jordan. An Introduction to Graphical Models. Technical Report, 2003.
  • Nan M. Laird and James Ware. Random-effects models for longitudinal data. Biometrics, 963-974, 1982.
  • Greg Wei and Martin A. Tanner. A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms. Journal of the American Statistical Association, 699-704, 1990.
 

以上






TensorFlow Probability : Tutorials : ガウス分布の階乗混合

TensorFlow Probability : Tutorials : ガウス分布の階乗混合 (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 11/02/2017

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです:

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : ガウス分布の階乗混合

この notebook では次のように定義されるガウス分布の階乗混合分布からサンプリングするために TensorFlow Probability (TFP) をどのように使用するかを示します :

$$
p(x_1, …, x_n) = \prod_i p_i(x_i)
$$

ここで:
$$
\begin{align*} p_i &\equiv \frac{1}{K}\sum_{i=1}^K \pi_{ik}\,\text{Normal}\left(\text{loc}=\mu_{ik},\, \text{scale}=\sigma_{ik}\right)\\1&=\sum_{k=1}^K\pi_{ik}, \forall i.\hphantom{MMMMMMMMMMM}\end{align*}
$$

各変数 $x_i$ はガウス分布の混合としてモデル化され、そして総ての $n$ 変数に渡る結合分布はこれらの密度の積です。

データセット $x^{(1)}, …, x^{(T)}$ が与えられたとき、各データポイント $x^{(j)}$ をガウス分布の階乗混合としてモデル化します :
$$
p(x^{(j)}) = \prod_i p_i (x_i^{(j)})
$$

階乗混合は小さい数のパラメータと巨大な数のモードを持つ分布を作成する単純な方法です。

import tensorflow as tf
import numpy as np
import tensorflow_probability as tfp
import matplotlib.pyplot as plt
import seaborn as sns
tfd = tfp.distributions

# Use try/except so we can easily re-execute the whole notebook.
try:
  tf.enable_eager_execution()
except:
  pass

 

TFP を使用してガウス分布の階乗混合を構築する

num_vars = 2        # Number of variables (`n` in formula).
var_dim = 1         # Dimensionality of each variable `x[i]`.
num_components = 3  # Number of components for each mixture (`K` in formula).
sigma = 5e-2        # Fixed standard deviation of each component.

# Set seed. Remove this line to generate different mixtures!
tf.set_random_seed(77)

# Choose some random (component) modes.
component_mean = tfd.Uniform().sample([num_vars, num_components, var_dim])

factorial_mog = tfd.Independent(
   tfd.MixtureSameFamily(
       # Assume uniform weight on each component.
       mixture_distribution=tfd.Categorical(
           logits=tf.zeros([num_vars, num_components])),
       components_distribution=tfd.MultivariateNormalDiag(
           loc=component_mean, scale_diag=[sigma])),
   reinterpreted_batch_ndims=1)

tfd.Independent の使用に注意してください。この「メタ分布」は右端の reinterpreted_batch_ndims バッチ次元に渡る log_prob 計算内で reduce_sum を適用します。私達のケースでは、log_prob を計算するときこれはバッチ次元だけを残して variables 次元を総計します。

 

密度をプロットする

点のグリッド上で密度を計算して赤い星でモデルの位置を示します。階乗混合の各モードはガウス分布の基礎的な individual-variable 混合からのモードのペアに対応します。下のプロットでは 9 モードを見て取れますが、6 パラメータだけが必要でした ($x_1$ のモードの位置を指定するために 3 つと、$x_2$ のモードの位置を指定するために 3 つです)。対照的に、2d 空間 $(x_1, x_2)$ のガウス分布の混合は 9 モードを指定するために 2 * 9 = 18 パラメータを必要とするでしょう。

plt.figure(figsize=(6,5))

# Compute density.
nx = 250 # Number of bins per dimension.
x = np.linspace(-3 * sigma, 1 + 3 * sigma, nx).astype('float32')
vals = tf.reshape(tf.stack(np.meshgrid(x, x), axis=2), (-1, num_vars, var_dim))
probs = factorial_mog.prob(vals).numpy().reshape(nx, nx)

# Display as image.
from matplotlib.colors import ListedColormap
cmap = ListedColormap(sns.color_palette("Blues", 256))
p = plt.pcolor(x, x, probs, cmap=cmap)
ax = plt.axis('tight');

# Plot locations of means.
means_np = component_mean.numpy().squeeze()
for mu_x in means_np[0]:
  for mu_y in means_np[1]:
    plt.scatter(mu_x, mu_y, s=150, marker='*', c='r', edgecolor='none');
plt.axis(ax);

plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.title('Density of factorial mixture of Gaussians');

Plot samples and marginal density estimates

samples = factorial_mog.sample(1000).numpy()

g = sns.jointplot(
    x=samples[:, 0, 0],
    y=samples[:, 1, 0],
    kind="scatter",
    stat_func=None,
    marginal_kws=dict(bins=50))
g.set_axis_labels("$x_1$", "$x_2$");

 

以上






TensorFlow Probability : Edward2

TensorFlow Probability : Edward2 (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
更新日時 : 11/18/2018
作成日時 : 11/01/2018

* 本ページは、tensorflow/probability の Edward2 の README.md を翻訳した上で適宜、補足説明したものです:

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : Edward2

Edward2 は TensorFlow と Python の確率的プログラミング言語です。それは TensorFlow エコシステムを拡張し、その結果、モデルを確率的プログラムとして宣言して柔軟な訓練、潜在変数推論、そして予測のためにモデルの計算を操作できます。

貴方は Edward からアップグレードしていますか?ガイド Upgrading_from_Edward_to_Edward2.md をチェックしてください。

 

1. 確率的プログラムとしてのモデル

確率変数

Edward2 では、確率的モデルの構造を指定するために RandomVariables を使用します。確率変数 rv は確率分布 (rv.distribution) を運搬 (= carry) し、これは log_prob と sample のような確率変数のメソッドを管理する TensorFlow Distribution インスタンスです。

確率変数は TensorFlow Distributions のように形成されます。

from tensorflow_probability import edward2 as ed

normal_rv = ed.Normal(loc=0., scale=1.)
## <ed.RandomVariable 'Normal/' shape=() dtype=float32>
normal_rv.distribution.log_prob(1.231)
## <tf.Tensor 'Normal/log_prob/sub:0' shape=() dtype=float32>

dirichlet_rv = ed.Dirichlet(concentration=tf.ones([2, 10]))
## <ed.RandomVariable 'Dirichlet/' shape=(2, 10) dtype=float32>

デフォルトでは、確率変数 rv のインスタンス化は tensor rv.value ~ rv.distribution.sample() を形成するためにサンプリング op を作成します。サンプルのデフォルト数 (rv への sample_shape 引数を通して制御可能です) は一つで、そしてオプションの value 引数が提供されれば、サンプリング op は作成されません。確率変数は TensorFlow op と相互作用できます: TF ops はサンプル上で動作します。

x = ed.Normal(loc=tf.zeros(10), scale=tf.ones(10))
y = 5.
x + y, x / y
## (<tf.Tensor 'add:0' shape=(10,) dtype=float32>,
##  <tf.Tensor 'div:0' shape=(10,) dtype=float32>)
tf.tanh(x * y)
## <tf.Tensor 'Tanh:0' shape=(10,) dtype=float32>
x[2]  # 3rd normal rv
## <tf.Tensor 'strided_slice:0' shape=() dtype=float32>

 

確率モデル

Edward2 の確率モデルは一つまたはそれ以上の RandomVariables をインスタンス化する Python 関数として表現されます。典型的には、関数 (「プログラム」) は生成過程を実行してサンプルを返します。関数への入力はモデルがその上で条件付けられる値として考えることができます。

下でベイズロジスティック回帰を書きますが、そこでは特徴、係数、そして切片が与えられたとき二値出力結果が生成されます。係数と切片に渡る事前分布があります。関数の実行は TensorFlow グラフに演算を追加します、そして TensorFlow session で結果ノードを求めると事前分布から係数と切片をサンプリングし、そしてこれらのサンプルを出力結果を計算するために使用します。

def logistic_regression(features):
  """Bayesian logistic regression p(y | x) = int p(y | x, w, b) p(w, b) dwdb."""
  coeffs = ed.Normal(loc=tf.zeros(features.shape[1]), scale=1., name="coeffs")
  intercept = ed.Normal(loc=0., scale=1., name="intercept")
  outcomes = ed.Bernoulli(
      logits=tf.tensordot(features, coeffs, [[1], [0]]) + intercept,
      name="outcomes")
  return outcomes

num_features = 10
features = tf.random_normal([100, num_features])
outcomes = logistic_regression(features)

# Execute the model program, returning a sample np.ndarray of shape (100,).
with tf.Session() as sess:
  outcomes_ = sess.run(outcomes)

Edward2 プログラムはまたデータを直接的にモデル化したものを越えて分布を表わすこともできます。例えば、下では学習可能な分布をそれをロジスティック回帰事後分布に近似する意図で書きます。

import tensorflow_probability as tfp

def logistic_regression_posterior(num_features):
  """Posterior of Bayesian logistic regression p(w, b | {x, y})."""
  posterior_coeffs = ed.MultivariateNormalTriL(
      loc=tf.get_variable("coeffs_loc", [num_features]),
      scale_tril=tfp.trainable_distributions.tril_with_diag_softplus_and_shift(
          tf.get_variable("coeffs_scale", [num_features*(num_features+1) / 2])),
      name="coeffs_posterior")
  posterior_intercept = ed.Normal(
      loc=tf.get_variable("intercept_loc", []),
      scale=tfp.trainable_distributions.softplus_and_shift(
          tf.get_variable("intercept_scale", [])),
      name="intercept_posterior")
  return coeffs, intercept

coeffs, intercept = logistic_regression_posterior(num_features)

# Execute the program, returning a sample
# (np.ndarray of shape (55,), np.ndarray of shape ()).
with tf.Session() as sess:
  sess.run(tf.global_variables_initializer())
  posterior_coeffs_, posterior_intercept_ = sess.run(
      [posterior_coeffs, posterior_intercept])

上の例で tfp.trainable_distributions の使用方法に注意してください。それは (正規分布のスケールのように) 制約付き空間にあるパラメータを持つ分布のためのユティリティを提供します。

変分プログラムとしてのテクニックの完全な例については、確率的 PCA チュートリアル を見てください。

 

モデル計算を操作する

Interceptor

確率モデルの訓練とテストは典型的には生成過程からの単なるサンプル以上を必要とします。柔軟な訓練とテストを有効にするために、interceptor を使用してモデルの計算を操作します。

interceptor はもう一つの関数 f とその引数 *args, **kwargs 上で動作する関数です。それは出力を返す前に様々な計算を遂行します (典型的には f(*args, **kwargs) : その関数自身を適用する結果です)。ed.interception コンテキスト・マネージャは interceptor をスタック上にプッシュし、任意のインターセプト可能な関数はスタックによりインターセプトされます。総ての確率変数コンストラクタはインターセプト可能です。

下ではロジスティック回帰モデルの生成過程をインターセプトします。特に、その事前分布ではなく学習された事後分布平均で予測を行ないます。

def set_prior_to_posterior_mean(f, *args, **kwargs):
  """Forms posterior predictions, setting each prior to its posterior mean."""
  name = kwargs.get("name")
  if name == "coeffs":
    return posterior_coeffs.distribution.mean()
  elif name == "intercept":
    return posterior_intercept.distribution.mean()
  return f(*args, **kwargs)

with ed.interception(set_prior_to_posterior_mean):
  predictions = logistic_regression(features)

training_accuracy = (
    tf.reduce_sum(tf.cast(tf.equal(predictions, outcomes), tf.float32)) /
    tf.cast(tf.shape(outcomes), tf.float32))

 

プログラム変換

interceptor を使用して、プログラム変換を適用することもできます、これはモデルの一つの表現から他のものへマップします。これは (ダウンストリームなユースケースに依拠する) 異なるモデル特性への便利なアクセスを提供します。

例えば、マルコフ連鎖モンテカルロ・アルゴリズムはしばしば入力としてモデルの対数同時確率関数を必要とします。下では生成過程を指定するベイズロジスティック回帰問題を取りそしてその対数同時確率関数を得るために組み込みの ed.make_log_joint 変換を適用します。対数同時関数は入力として生成プログラムの元の入力とプログラムの確率変数を取ります。それは確率変数対数確率に渡り総計したスカラー Tensor を返します。

私達の例では、特徴と出力結果は固定され、coeffs と intercept の事後分布からサンプルをドローするためにハミルトニアン・モンテカルロを使用することを望みます。この使用のために、target_log_prob_fn を作成します、これは単に coeffs と intercept を引数として取り入力特徴をピンで止めてその既知の値として rv 結果を出力します。

import tensorflow_probability as tfp

# Set up training data.
features = tf.random_normal([100, 55])
outcomes = tf.random_uniform([100], minval=0, maxval=2, dtype=tf.int32)

# Pass target log-probability function to MCMC transition kernel.
log_joint = ed.make_log_joint_fn(logistic_regression)

def target_log_prob_fn(coeffs, intercept):
  """Target log-probability as a function of states."""
  return log_joint(features,
                   coeffs=coeffs,
                   intercept=intercept,
                   outcomes=outcomes)

hmc_kernel = tfp.mcmc.HamiltonianMonteCarlo(
    target_log_prob_fn=target_log_prob_fn,
    step_size=0.1,
    num_leapfrog_steps=5)
states, kernel_results = tfp.mcmc.sample_chain(
    num_results=1000,
    current_state=[tf.random_normal([55]), tf.random_normal([])],
    kernel=hmc_kernel,
    num_burnin_steps=500)

with tf.Session() as sess:
  states_, results_ = sess.run([states, kernels_results])

返された states_[0] と states_[1] はそれぞれ coeffs と intercept のための 1,000 事後サンプルを含みます。これは例えば、新しいデータ上でモデルの事後予測を評価するために使用されるかもしれません。

 

Examples

end-to-end サンプルとチュートリアルに興味がありますか?tensorflow_probability/examples/ を見てください。

 

以上






TensorFlow Probability : Tutorials : ガウス過程 潜在変数 (状態空間) モデル

TensorFlow Probability : Tutorials : ガウス過程 潜在変数 (状態空間) モデル (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 10/30/2018

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです:

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : ガウス過程 潜在変数 (状態空間) モデル

潜在変数モデルは高次元データの隠れ構造を捕捉することを企てます。サンプルは主成分分析 (PCA) と因子分析を含みます。ガウス過程は「ノンパラメトリック」モデルでこれは局所相関構造と不確かさを柔軟に捕捉できます。ガウス過程潜在変数モデル (Lawrence, 2004) はこれらの概念を結合します。

 

背景: ガウス過程

ガウス過程は確率変数の任意のコレクションで任意の有限個の部分集合に渡る周辺分布が多変量正規分布であるようなものです。回帰のコンテキストにおける GP の詳細な見え方は、TensorFlow Probability のガウス過程回帰 をチェックしてください。

GP が構成されるコレクションの確率変数の各々をラベル付けするいわゆるインデックス集合を使用します。有限のインデックス集合の場合には、単に多変量正規分布を得ます。けれども無限コレクションを考慮するとき、GP は非常に興味深いです。$\mathbb{R}^D$ のようなインデックス集合の場合、そこでは $D$ 次元空間の総ての点について確率変数を持ち、GP はランダム関数に渡る分布として考えることができます。そのような GP からの単一のドローは、それが現実化するならば、$\mathbb{R}^D$ の総ての点に (結合的に正規分布に従う) 値を割り当てるでしょう。この colab では、ある $\mathbb{R}^D$ に渡る GP に注目します。

正規分布はそれらの最初と 2 番目の統計値で完全に決定されます — 実際に、正規分布を定義する一つの方法はその高次キュムラントが総てゼロであるものとしてです。これは GP の場合にもまた、当てはまります: 平均と共分散* を記述することで GP を完全に指定します。有限次元多変量正規分布については、平均はベクトルで共分散は正方 (= square) で、対称な正定値行列であることを思い出してください。無限次元 GP では、これらの構造はインデックス集合の各点で定義される平均関数, $m : \mathbb{R}^D \to \mathbb{R}$、そして共分散「カーネル」関数, $k : \mathbb{R}^D \times \mathbb{R}^D \to \mathbb{R}$ に一般化されます。カーネル関数は 正定値 であることが要求されます、これは本質的には、有限集合の点に制限されれば、それは正定値行列を生成することを言っています。

GP の構造の殆どはその共分散カーネル関数に由来します — この関数はサンプリングされた関数の値が近くの (or それほど近くはない) 点に渡りどのように変化するかを記述します。異なる共分散関数は滑らかさの異なる度合いに働きかけます。一つの一般的に使用されるカーネル関数は「指数二次 (= exponentiated quadratic)」 (a.k.a. ガウシアン、二乗指数 (= squared exponential)、動径基底関数 (= radial basis function)), $k(x, x’) = \sigma^2 e^{(x – x^2) / \lambda^2}$ です。他の例は伝統的なテキスト Gaussian Processes for Machine Learning 内に加えて、David Duvenaud の kernel cookbook ページ で概説されます。

 

GP の適用: 回帰と潜在変数モデル

GP を使用できる一つの方法は回帰のためです: 入力 $\{x_i\}_{i=1}^N$ (インデックス集合の要素)と観測 $\{y_i\}_{i=1}^N$ の形式で多くの観測されたデータが与えられたとき、点 $\{x_j^*\}_{j=1}^M$ の新しい集合における事後予測分布を形成するためにこれらを使用できます。分布は総て正規分布ですので、これは何某かの straightforward な線形代数に要約されます (しかし注意してください: 必要な計算はデータポイントの数で 3 次の実行時間を持ちそしてデータポイントの数で 2 次の空間を必要とします — これは GP の使用における主要な制限的要因で現在の研究の大半は正確な事後推論への計算的に実行可能な代替案に焦点を合わせています)。TFP における GP 回帰 で GP 回帰を詳しくカバーします。

GP を使用できるもう一つの方法は潜在変数モデルとしてです: 高次元観測 (e.g., 画像) のコレクションが与えられたとき、ある低次元潜在構造を推測することができます。潜在構造を条件とするとき、巨大な数の出力 (画像のピクセル) は互いに独立であることを仮定します。このモデルの訓練は以下から成ります :

  1. モデル・パラメータを最適化します (カーネル関数パラメータと、e.g. 観測ノイズ分散)、そして
  2. 各訓練観測 (画像) に対して、インデックス集合の対応する点位置を見つけます。最適化の総てはデータの周辺対数尤度を最大化することにより成されます。

 

インポート

import numpy as np
import tensorflow as tf
from tensorflow_probability import distributions as tfd
from tensorflow_probability import positive_semidefinite_kernels as tfk

numpy と matplotlib から対話的な名前空間を取り込みます。

def reset_session():
  """Creates a new global, interactive session in Graph-mode."""
  global sess
  try:
    tf.reset_default_graph()
    sess.close()
  except:
    pass
  sess = tf.InteractiveSession()

reset_session()

 

MNIST データをロードする

# Load the MNIST data set and isolate a subset of it.
(x_train, y_train), (_, _) = tf.keras.datasets.mnist.load_data()
N = 1000
small_x_train = x_train[:N, ...].astype(np.float64) / 256.
small_y_train = y_train[:N]

 

訓練可能な変数を準備する

3 モデル・パラメータと潜在入力を結合的に訓練していきます。

# Create some trainable model parameters. We also constrain them to be strictly
# positive
amplitude = np.finfo(np.float64).eps + tf.nn.softplus(
    tf.get_variable(name='amplitude',
                    dtype=tf.float64,
                    initializer=np.float64(1.)))

length_scale = np.finfo(np.float64).eps + tf.nn.softplus(
    tf.get_variable(name='length_scale',
                    dtype=tf.float64, 
                    initializer=np.float64(1.)))

observation_noise_variance = np.finfo(np.float64).eps + tf.nn.softplus(
    tf.get_variable(name='observation_noise_variance',
                    dtype=tf.float64, 
                    initializer=np.float64(1.)))
# We need to flatten the images and, somewhat unintuitively, transpose from
# shape [100, 784] to [784, 100]. This is because the 784 pixels will be
# treated as *independent* conditioned on the latent inputs, meaning we really
# have a batch of 784 GP's with 100 index_points.
observations_ = small_x_train.reshape(N, -1).transpose()

# Create a collection of N 2-dimensional index points that will represent our
# latent embeddings of the data. (Lawrence, 2004) prescribes initializing these
# with PCA, but a random initialization actually gives not-too-bad results, so
# we use this for simplicity. For a fun exercise, try doing the
# PCA-initialization yourself!
init_ = np.random.normal(size=(N, 2))
latent_index_points = tf.get_variable(
    name='latent_index_points',
    dtype=tf.float64,
    initializer=init_)

 

モデルと訓練 ops を構築する

# Create our kernel and GP distribution
kernel = tfk.ExponentiatedQuadratic(amplitude, length_scale)

gp = tfd.GaussianProcess(
    kernel=kernel,
    index_points=latent_index_points,
    observation_noise_variance=observation_noise_variance)

# Compute the log probability density of each of the observations. This is what
# we will maximize below.
log_probs = gp.log_prob(observations_, name='log_prob')

# Our training loss is the negative log likelihood of the observations. Since
# the observations are IID, the log likelihoods are summed. We normalize the sum
# (take the average) to reduce the risk of extreme values subject to floating
# point issues.
loss = -tf.reduce_mean(log_probs)
optimizer = tf.train.AdamOptimizer(learning_rate=.1)
train_op = optimizer.minimize(loss)

 

予測モデルとサンプリング ops を構築する

# We'll draw samples at evenly spaced points on a 10x10 grid in the latent
# input space. 
sample_grid_points = 10
grid_ = np.linspace(-4, 4, sample_grid_points).astype(np.float64)
# Create a 10x10 grid of 2-vectors, for a total shape [10, 10, 2]
grid_ = np.stack(np.meshgrid(grid_, grid_), axis=-1)

# This part's a bit subtle! What we defined above was a batch of 784 (=28x28)
# independent GP distributions over the input space. Each one corresponds to a
# single pixel of an MNIST image. Now what we'd like to do is draw 100 (=10x10)
# *independent* samples, each one separately conditioned on all the observations
# as well as the learned latent input locations above.
#
# The GP regression model below will define a batch of 784 independent
# posteriors. We'd like to get 100 independent samples each at a different
# latent index point. We could loop over the points in the grid, but that might
# be a bit slow. Instead, we can vectorize the computation by tacking on *even
# more* batch dimensions to our GaussianProcessRegressionModel distribution.
# In the below grid_ shape, we have concatentaed
#   1. batch shape: [sample_grid_points, sample_grid_points, 1]
#   2. number of examples: [1]
#   3. number of latent input dimensions: [2]
# The `1` in the batch shape will broadcast with 784. The final result will be
# samples of shape [10, 10, 784, 1]. The `1` comes from the "number of examples"
# and we can just `np.squeeze` it off.
grid_ = grid_.reshape(sample_grid_points, sample_grid_points, 1, 1, 2)

# Create the GPRegressionModel instance which represents the posterior
# predictive at the grid of new points.
gprm = tfd.GaussianProcessRegressionModel(
    kernel=kernel,
    index_points=grid_,
    observation_index_points=latent_index_points,
    observations=observations_)

# Create a sample operation in the graph.
samples = gprm.sample()

 

訓練と結果の潜在埋め込みのプロット

# Initialize variables and train!
sess.run(tf.global_variables_initializer())
num_iters = 100
log_interval = 20
lips_ = np.zeros((num_iters, N, 2), np.float64)
for i in range(num_iters):
  _, loss_, lips_[i] = sess.run([train_op, loss, latent_index_points])
  if i % log_interval == 0 or i + 1 == num_iters:
    print("Loss at step %d: %f" % (i, loss_))
Loss at step 0: 869611.863774
Loss at step 20: 330985.762144
Loss at step 40: -90885.405676
Loss at step 60: -182680.904868
Loss at step 80: -207781.362585
Loss at step 99: -212932.604949

 

結果をプロットする

# Plot the latent locations before and after training
plt.figure(figsize=(7, 7))
plt.title("Before training")
plt.grid('off')
plt.scatter(x=init_[:, 0], y=init_[:, 1],
           c=y_train[:N], cmap=plt.get_cmap('Paired'), s=50)
plt.show()

plt.figure(figsize=(7, 7))
plt.title("After training")
plt.grid('off')
plt.scatter(x=lips_[-1, :, 0], y=lips_[-1, :, 1],
           c=y_train[:N], cmap=plt.get_cmap('Paired'), s=50)
plt.show()

 

データ上で条件付けられたサンプルのドローと潜在埋め込み

潜在空間の 2-d グリッド上の 100 点でサンプリングします。

samples_ = sess.run(samples)

# Plot the grid of samples at new points. We do a bit of tweaking of the samples
# first, squeezing off extra 1-shapes and normalizing the values.
samples_ = np.squeeze(samples_)
samples_ = ((samples_ -
             samples_.min(-1, keepdims=True)) /
            (samples_.max(-1, keepdims=True) -
             samples_.min(-1, keepdims=True)))
samples_ = samples_.reshape(sample_grid_points, sample_grid_points, 28, 28)
samples_ = samples_.transpose([0, 2, 1, 3])
samples_ = samples_.reshape(28 * sample_grid_points, 28 * sample_grid_points)
plt.figure(figsize=(7, 7))
ax = plt.subplot()
ax.grid('off')
ax.imshow(-samples_, interpolation='none', cmap='Greys')
plt.show()

 

以上






TensorFlow Probability : Tutorials : ガウス過程回帰

TensorFlow Probability : Tutorials : ガウス過程回帰 (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 10/30/2018

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです:

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : ガウス過程回帰

この colab では、TensorFlow と TensorFlow Probability を使用してガウス過程回帰を探究します。ある知られている関数(s) からのノイズを持つ観測を生成して GP モデルをそれらのデータに適合させます。それから GP 事後分布からサンプリングされた関数値をそれらの領域でグリッドに渡りプロットします。

 

背景

$\mathcal{X}$ を任意の集合とします。ガウス過程 (GP) は $\mathcal{X}$ によりインデックス付けられた確率変数のコレクションで、$\{X_1, \ldots, X_n\} \subset \mathcal{X}$ を任意の有限な部分集合とする場合に、周辺密度関数 $p(X_1 = x_1, \ldots, X_n = x_n)$ が多変量ガウス分布であるようなものです。任意のガウス分布はその最初と 2 番目の中心モーメント (平均と共分散) で完全に指定され、そして GP のものも例外ではありません。GP はその平均関数 $\mu : \mathcal{X} \to \mathbb{R}$ と共分散関数 $k : \mathcal{X} \times \mathcal{X} \to \mathbb{R}$ によって完全に指定できます。GP の表現力の殆どは共分散関数の選択にカプセル化されています。様々な理由で、共分散関数はまたカーネル関数として参照されます。対照的 (= symmetirc) で正定値 (= positive-definite) であることだけが要求されます (Rasmussen & Williams の Ch. 4 参照)。下では指数二次 (= ExponentiatedQuadratic) 共分散カーネルを活用します。その形式は :
\[
k(x, x’) := \sigma^2 \exp \left( \frac{\|x – x’\|^2}{\lambda^2} \right)
\]

$\sigma^2$ は「振幅 (= amplitude)」そして $\lambda$ は長さスケールと呼ばれます。カーネル・パラメータは最大尤度最適化手続きを通して選択可能です。

GP からのフル・サンプルは空間全体 $\mathcal{X}$ に渡る実数値関数から成りそして実際には認識することは実行不可能です; しばしば (そこでサンプルを観測するために) 点の集合を選択してそしてそれらの点で関数値を描画します。これは適切な (有限次元の) 多変量ガウス分布からサンプリングすることで獲得されます。

上の定義に従えば、任意の有限次元多変量ガウス分布もまたガウス過程であることに注意してください。通常は、GP に言及するとき、インデックス集合はある $\mathbb{R}^n$ であることが暗黙的でここでも実際にこの仮定をします。

機械学習におけるガウス過程の一般的な応用はガウス過程回帰です。考え方は、有限個の点 $\{x_1, \ldots x_N\}$ における関数のノイズを伴う観測 $\{y_1, \ldots, y_N\}$ が与えられたとき未知の関数を推定することを望むことです。生成過程を推測します :
\[
\begin{align}
f \sim \: & \textsf{GaussianProcess}\left(
\text{mean_fn}=\mu(x),
\text{covariance_fn}=k(x, x’)\right) \\
y_i \sim \: & \textsf{Normal}\left(
\text{loc}=f(x_i),
\text{scale}=\sigma\right), i = 1, \ldots, N
\end{align}
\]

上で記されたように、サンプリングされた関数は計算することは不可能です、何故ならば無限の数の点におけるその値を必要とするからです。代わりに、多変量ガウス分布からの有限個のサンプルを考えます。

$$
\begin{gather}
\begin{bmatrix}
f(x_1) \\
\vdots \\
f(x_N)
\end{bmatrix}
\sim
\textsf{MultivariateNormal} \left( \:
\text{loc}=
\begin{bmatrix}
\mu(x_1) \\
\vdots \\
\mu(x_N)
\end{bmatrix} \:,\:
\text{scale}=
\begin{bmatrix}
k(x_1, x_1) & \cdots & k(x_1, x_N) \\
\vdots & \ddots & \vdots \\
k(x_N, x_1) & \cdots & k(x_N, x_N) \\
\end{bmatrix}^{1/2}
\: \right)
\end{gather} \\
y_i \sim \textsf{Normal} \left(
\text{loc}=f(x_i),
\text{scale}=\sigma
\right)
$$

共分散行列上の指数 $\frac{1}{2}$ に注意してください: これはコレスキー分解を表しています。コレスキーを計算することは必要です、何故ならば MVN は location-scale family 分布だからです。不幸にも、コレスキー分解は計算的に高価で、$O(N^3)$ 時間と $O(N^2)$ 空間がかかります。大半の GP 文献はこの表面上は無害で小さな指数を処理することに焦点を置いています。

事前分布平均関数を定数、しばしばゼロに取ることは一般的です。また、幾つかの表記法の慣習は便利です。サンプリングされた関数値の有限ベクトルのために $\mathbf{f}$ としばしば書きます。多くの興味深い表記法が、入力ペアへの $k$ の適用に由来する共分散行列のために使用されます。(Quiñonero-Candela, 2005) に従って、行列の成分は特定の入力点における関数値の共分散であることに注意します。こうして共分散行列を $K_{AB}$ として示すことができます、そこでは $A$ と $B$ は与えられた行列次元に沿った関数値のコレクションの幾つかの indicator です。

例えば、含蓄される潜在関数値 $\mathbf{f}$ を持つ観測されたデータ $(\mathbf{x}, \mathbf{y})$ が与えられたとき、次のように書けます :

$$
K_{\mathbf{f},\mathbf{f}} = \begin{bmatrix}
k(x_1, x_1) & \cdots & k(x_1, x_N) \\
\vdots & \ddots & \vdots \\
k(x_N, x_1) & \cdots & k(x_N, x_N) \\
\end{bmatrix}
$$

同様に、以下の内のように、入力の集合を混合できます :

$$
K_{\mathbf{f},*} = \begin{bmatrix}
k(x_1, x^*_1) & \cdots & k(x_1, x^*_T) \\
\vdots & \ddots & \vdots \\
k(x_N, x^*_1) & \cdots & k(x_N, x^*_T) \\
\end{bmatrix}
$$

ここで $N$ 訓練入力、そして $T$ テスト入力があると仮定します。そして上の生成過程は次のように簡潔に書くことができます :

$$
\begin{align}
\mathbf{f} \sim \: & \textsf{MultivariateNormal} \left(
\text{loc}=\mathbf{0},
\text{scale}=K_{\mathbf{f},\mathbf{f}}^{1/2}
\right) \\
y_i \sim \: & \textsf{Normal} \left(
\text{loc}=f_i,
\text{scale}=\sigma \right), i = 1, \ldots, N
\end{align}
$$

最初の行のサンプリング演算は (一つの) 多変量正規分布からの $N$ 関数値の有限集合を生成します — 上の GP ドロー記法内のような関数全体 (= an entire function) ではありません。2 番目の行は固定された観測ノイズ $\sigma^2$ を持ち、様々な関数値を中心とする単変量正規分布 (= univariate Gaussians) からの $N$ ドローのコレクションを表します。

上の適所にある生成モデルで、事後分布問題を考えることに進めます。これは上の仮定からの観測されたノイズを持つデータ上で条件付けられた、テスト点の新しい集合における関数値に渡る事後分布を生成します。

上の適所にある表記で、相当する入力と訓練データ上で条件付けられた未来 (ノイズを持つ) 観測に渡る事後予測分布を次のように簡潔に書くことができます (より詳細については、Rasmussen & Williams の 2.2 参照)。
$$
\mathbf{y}^* \mid \mathbf{x}^*, \mathbf{x}, \mathbf{y} \sim \textsf{Normal} \left(
\text{loc}=\mathbf{\mu}^*,
\text{scale}=(\Sigma^*)^{1/2}
\right)
$$

ここで
$$
\mathbf{\mu}^* = K_{*,\mathbf{f}}\left(K_{\mathbf{f},\mathbf{f}} + \sigma^2 I \right)^{-1} \mathbf{y}
$$

そして
$$
\Sigma^* = K_{*,*} – K_{*,\mathbf{f}}
\left(K_{\mathbf{f},\mathbf{f}} + \sigma^2 I \right)^{-1} K_{\mathbf{f},*}
$$

 

インポート

import numpy as np
import tensorflow as tf
from mpl_toolkits.mplot3d import Axes3D
from tensorflow_probability import distributions as tfd
from tensorflow_probability import positive_semidefinite_kernels as tfk

%pylab inline
# Configure plot defaults
plt.rcParams['axes.facecolor'] = 'white'
plt.rcParams['grid.color'] = '#666666'
%config InlineBackend.figure_format = 'png'

numpy と matplotlib から対話的な名前空間を取り込んでいます。

 

例: ノイズを持つ正弦波データ上の正確な GP 回帰

ここではノイズを持つ正弦波から訓練データを生成し、それから GP 回帰モデルの事後分布から多くの曲線をサンプリングします。カーネル・ハイパーパラメータを最適化するために Adam を使用します (事前分布のもとのデータの負の対数尤度を最小化します)。訓練している曲線に続いて、正解関数と事後サンプルをプロットします。

def reset_session():
  """Creates a new global, interactive session in Graph-mode."""
  global sess
  try:
    tf.reset_default_graph()
    sess.close()
  except:
    pass
  sess = tf.InteractiveSession()

reset_session()
def sinusoid(x):
  return np.sin(3 * np.pi * x[..., 0])

def generate_1d_data(num_training_points, observation_noise_variance):
  """Generate noisy sinusoidal observations at a random set of points.

  Returns:
     observation_index_points, observations
  """
  index_points_ = np.random.uniform(-1., 1., (num_training_points, 1))
  index_points_ = index_points_.astype(np.float64)
  # y = f(x) + noise
  observations_ = (sinusoid(index_points_) +
                   np.random.normal(loc=0,
                                    scale=np.sqrt(observation_noise_variance),
                                    size=(num_training_points)))
  return index_points_, observations_
# Generate training data with a known noise level (we'll later try to recover
# this value from the data).
NUM_TRAINING_POINTS = 100
observation_index_points_, observations_ = generate_1d_data(
    num_training_points=NUM_TRAINING_POINTS,
    observation_noise_variance=.1)
# Create the trainable model parameters, which we'll subsequently optimize.
# Note that we constrain them to be strictly positive.
amplitude = (np.finfo(np.float64).tiny +
             tf.nn.softplus(tf.Variable(initial_value=1.,
                                        name='amplitude',
                                        dtype=np.float64)))
length_scale = (np.finfo(np.float64).tiny +
                tf.nn.softplus(tf.Variable(initial_value=1.,
                                           name='length_scale',
                                           dtype=np.float64)))

observation_noise_variance = (
    np.finfo(np.float64).tiny +
    tf.nn.softplus(tf.Variable(initial_value=1e-6,
                               name='observation_noise_variance',
                               dtype=np.float64)))
# Create the covariance kernel, which will be shared between the prior (which we
# use for maximum likelihood training) and the posterior (which we use for
# posterior predictive sampling)
kernel = tfk.ExponentiatedQuadratic(amplitude, length_scale)
# Create the GP prior distribution, which we will use to train the model
# parameters.
gp = tfd.GaussianProcess(
    kernel=kernel,
    index_points=observation_index_points_,
    observation_noise_variance=observation_noise_variance)

# This lets us compute the log likelihood of the observed data. Then we can
# maximize this quantity to find optimal model parameters.
log_likelihood = gp.log_prob(observations_)
# Define the optimization ops for maximizing likelihood (minimizing neg
# log-likelihood!)
optimizer = tf.train.AdamOptimizer(learning_rate=.01)
train_op = optimizer.minimize(-log_likelihood)
# Now we optimize the model parameters.
num_iters = 1000
# Store the likelihood values during training, so we can plot the progress
lls_ = np.zeros(num_iters, np.float64)
sess.run(tf.global_variables_initializer())
for i in range(num_iters):
  _, lls_[i] = sess.run([train_op, log_likelihood])

[amplitude_,
 length_scale_,
 observation_noise_variance_] = sess.run([
    amplitude,
    length_scale,
    observation_noise_variance])
print('Trained parameters:'.format(amplitude_))
print('amplitude: {}'.format(amplitude_))
print('length_scale: {}'.format(length_scale_))
print('observation_noise_variance: {}'.format(observation_noise_variance_))
Trained parameters:
amplitude: 0.737212035424
length_scale: 0.120468090569
observation_noise_variance: 0.0978894813992
# Plot the loss evolution
plt.figure(figsize=(12, 4))
plt.plot(lls_)
plt.xlabel("Training iteration")
plt.ylabel("Log marginal likelihood")
plt.show()

 

# Having trained the model, we'd like to sample from the posterior conditioned
# on observations. We'd like the samples to be at points other than the training
# inputs.
predictive_index_points_ = np.linspace(-1.2, 1.2, 200, dtype=np.float64)
# Reshape to [200, 1] -- 1 is the dimensionality of the feature space.
predictive_index_points_ = predictive_index_points_[..., np.newaxis]

gprm = tfd.GaussianProcessRegressionModel(
    kernel=kernel,  # Reuse the same kernel instance, with the same params
    index_points=predictive_index_points_,
    observation_index_points=observation_index_points_,
    observations=observations_,
    observation_noise_variance=observation_noise_variance,
    predictive_noise_variance=0.)

# Create op to draw  50 independent samples, each of which is a *joint* draw
# from the posterior at the predictive_index_points_. Since we have 200 input
# locations as defined above, this posterior distribution over corresponding
# function values is a 200-dimensional multivariate Gaussian distribution!
num_samples = 50
samples = gprm.sample(num_samples)
# Draw samples and visualize.
samples_ = sess.run(samples)

# Plot the true function, observations, and posterior samples.
plt.figure(figsize=(12, 4))
plt.plot(predictive_index_points_, sinusoid(predictive_index_points_),
         label='True fn')
plt.scatter(observation_index_points_[:, 0], observations_,
            label='Observations')
for i in range(num_samples):
  plt.plot(predictive_index_points_, samples_[i, :], c='r', alpha=.1,
           label='Posterior Sample' if i == 0 else None)
leg = plt.legend()
for lh in leg.legendHandles: 
    lh.set_alpha(1)
plt.xlabel(r"Index points ($\mathbb{R}^1)")
plt.ylabel("Observation space")
plt.show()

 
Note: 上のコードを何回か実行する場合、時にはそれは素晴らしく見えて他の回にはそれは酷く見えます!パラメータを訓練する最大尤度は非常に敏感で時々貧弱なモデルに収束します。最善のアプローチはモデル・ハイパーパラメータを周辺化するために MCMC を使用することです。

 

以上






TensorFlow Probability : Tutorials : TensorFlow 分布 shape の理解

TensorFlow Probability : Tutorials : TensorFlow 分布 shape の理解 (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 10/25/2018

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです:

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : TensorFlow 分布 shape を理解する

import collections

import tensorflow as tf
import tensorflow_probability as tfp
tfd = tfp.distributions
tfb = tfp.bijectors

tfe = tf.contrib.eager
tfe.enable_eager_execution()

from __future__ import print_function

 

基本

TensorFlow 分布 shape に関連する 3 つの重要な概念があります :

  • 事象 (= event) shape は分布からの単一のドローの shape を記述します; それは次元に渡り従属的かもしれません。スカラー分布については事象 shape は [] です。5-次元 MultivariateNormal については、事象 shape は [5] です。
  • バッチ shape は独立で、同一分布ではないドローを記述します, aka, 分布の「バッチ」です。
  • サンプル shape は分布族からのバッチの独立同分布なドローを記述します。

サンプル shape が sample か log_prob への特定の呼び出しに関係する一方で、事象 shape とバッチ shape は Distribution オブジェクトの特性です。

この notebook の目的は例を通してこれらの概念を示すことですので、これが直ちに明らかでなくても、don’t worry!

 

TensorFlow Eager 上のノート

notebook 全体は TensorFlow Eager を使用して書かれています。提示される概念のどれも Eager に依拠していません、けれども Eager では、分布のバッチと事象 shape は Python で Distribution オブジェクトが作成されたときに評価されます (そして従って知られます)。その一方でグラフ (non-Eager モード) では、グラフが実行されるまで事象とバッチ shape が未確定の分布を定義することが可能です。

 

スカラー分布

上で言及したように、Distribution オブジェクトは事象とバッチ shape を定義しました。分布を記述するユティリティから始めましょう :

def describe_distributions(distributions):
  print('\n'.join([str(d) for d in distributions]))

このセクションではスカラー分布を探究します: [] の事象 shape を持つ分布です。典型的な例は Poisson 分布で、rate により指定されます:

poisson_distributions = [
    tfd.Poisson(rate=1., name='One Poisson Scalar Batch'),
    tfd.Poisson(rate=[1., 10., 100.], name='Three Poissons'),
    tfd.Poisson(rate=[[1., 10., 100.,], [2., 20., 200.]],
                name='Two-by-Three Poissons'),
    tfd.Poisson(rate=[1.], name='One Poisson Vector Batch'),
    tfd.Poisson(rate=[[1.]], name='One Poisson Expanded Batch')
]

describe_distributions(poisson_distributions)
tfp.distributions.Poisson("One Poisson Scalar Batch/", batch_shape=(), event_shape=(), dtype=float32)
tfp.distributions.Poisson("Three Poissons/", batch_shape=(3,), event_shape=(), dtype=float32)
tfp.distributions.Poisson("Two-by-Three Poissons/", batch_shape=(2, 3), event_shape=(), dtype=float32)
tfp.distributions.Poisson("One Poisson Vector Batch/", batch_shape=(1,), event_shape=(), dtype=float32)
tfp.distributions.Poisson("One Poisson Expanded Batch/", batch_shape=(1, 1), event_shape=(), dtype=float32)

Poisson 分布はスカラー分布ですので、その事象 shape は常に [] です。より多くの rate を指定する場合、これらはバッチ shape に出現します。例の最後のペアは興味深いです: 単一の rate しかありませんが、その rate は non-empty shape で numpy 配列に埋め込まれていますので、その shape はバッチ shape になります。

標準正規分布もまたスカラーです。その事象 shape はちょうど Poisson のように [] ですが、ブロードキャスティングの最初の例を見るためにそれで遊びます。Normal は loc と sacle パラメータで指定されます:

normal_distributions = [
    tfd.Normal(loc=0., scale=1., name='Standard'),
    tfd.Normal(loc=[0.], scale=1., name='Standard Vector Batch'),
    tfd.Normal(loc=[0., 1., 2., 3.], scale=1., name='Different Locs'),
    tfd.Normal(loc=[0., 1., 2., 3.], scale=[[1.], [5.]],
               name='Broadcasting Scale')
]

describe_distributions(normal_distributions)
tfp.distributions.Normal("Standard/", batch_shape=(), event_shape=(), dtype=float32)
tfp.distributions.Normal("Standard Vector Batch/", batch_shape=(1,), event_shape=(), dtype=float32)
tfp.distributions.Normal("Different Locs/", batch_shape=(4,), event_shape=(), dtype=float32)
tfp.distributions.Normal("Broadcasting Scale/", batch_shape=(2, 4), event_shape=(), dtype=float32)

上の興味深い例は Broadcasting Scale 分布です。loc パラメータは shape [4] を持ち、scale パラメータは shape [2, 1] を持ちます。Numpy ブロードキャスティング・ルール を使用すれば、バッチ shape は [2, 4] です。”Broadcasting Scale” 分布を定義するための同値の (しかし洗練されておらず非推奨な) 方法は :

describe_distributions(
    [tfd.Normal(loc=[[0., 1., 2., 3], [0., 1., 2., 3.]],
                scale=[[1., 1., 1., 1.], [5., 5., 5., 5.]])])
tfp.distributions.Normal("Normal/", batch_shape=(2, 4), event_shape=(), dtype=float32)

ブロードキャスティング記法が何故有用かを見て取れます、それはまた頭痛とバグの源ですが。

 

スカラー分布のサンプリング

分布でできる 2 つの主要なことがあります: それらからサンプリングできて log_prob を計算できます。最初にサンプリングを探究しましょう。基本的なルールは分布からサンプリングするとき、結果の Tensor は shape [sample_shape, batch_shape, event_shape] を持ちます、そこでは batch_shape と event_shape は Distribution オブジェクトから提供されて、sample_shape は sample への呼び出しから提供されます。スカラー分布については、event_shape = [] ですので、sample から返される Tensor は shape [sample_shape, batch_shape] を持ちます。試してみましょう :

def describe_sample_tensor_shape(sample_shape, distribution):
    print('Sample shape:', sample_shape)
    print('Returned sample tensor shape:',
          distribution.sample(sample_shape).shape)

def describe_sample_tensor_shapes(distributions, sample_shapes):
    started = False
    for distribution in distributions:
      print(distribution)
      for sample_shape in sample_shapes:
        describe_sample_tensor_shape(sample_shape, distribution)
      print()

sample_shapes = [1, 2, [1, 5], [3, 4, 5]]
describe_sample_tensor_shapes(poisson_distributions, sample_shapes)
tfp.distributions.Poisson("One Poisson Scalar Batch/", batch_shape=(), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1,)
Sample shape: 2
Returned sample tensor shape: (2,)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5)

tfp.distributions.Poisson("Three Poissons/", batch_shape=(3,), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 3)
Sample shape: 2
Returned sample tensor shape: (2, 3)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 3)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 3)

tfp.distributions.Poisson("Two-by-Three Poissons/", batch_shape=(2, 3), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 2, 3)
Sample shape: 2
Returned sample tensor shape: (2, 2, 3)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 2, 3)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 2, 3)

tfp.distributions.Poisson("One Poisson Vector Batch/", batch_shape=(1,), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 1)
Sample shape: 2
Returned sample tensor shape: (2, 1)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 1)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 1)

tfp.distributions.Poisson("One Poisson Expanded Batch/", batch_shape=(1, 1), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 1, 1)
Sample shape: 2
Returned sample tensor shape: (2, 1, 1)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 1, 1)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 1, 1)
describe_sample_tensor_shapes(normal_distributions, sample_shapes)
tfp.distributions.Normal("Standard/", batch_shape=(), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1,)
Sample shape: 2
Returned sample tensor shape: (2,)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5)

tfp.distributions.Normal("Standard Vector Batch/", batch_shape=(1,), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 1)
Sample shape: 2
Returned sample tensor shape: (2, 1)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 1)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 1)

tfp.distributions.Normal("Different Locs/", batch_shape=(4,), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 4)
Sample shape: 2
Returned sample tensor shape: (2, 4)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 4)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 4)

tfp.distributions.Normal("Broadcasting Scale/", batch_shape=(2, 4), event_shape=(), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 2, 4)
Sample shape: 2
Returned sample tensor shape: (2, 2, 4)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 2, 4)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 2, 4)

sample について言うべきことはそれがおおよそ総てです: 返されるサンプル tensor は shape [sample_shape, batch_shape, event_shape] を持ちます。

 

スカラー分布のために log_prob を計算する

さて log_prob を見てみましょう、これは幾分扱いにくいです。log_prob は入力として (そこで分布のために log_prob を計算する) 位置を表わす (non-empty) tensor を取ります。最も straightforward なケースでは、この tensor は形式 [sample_shape, batch_shape, event_shape] の shape を持ち、そこでは batch_shape と event_shape は分布のバッチと事象 shape に一致します。スカラー分布については、event_shape = [] であることをもう一度思い出してください、従って入力 tensor は shape [sample_shape, batch_shape] を持ちます。この場合、shape [sample_shape, batch_shape] の tensor を取り戻します :

three_poissons = tfd.Poisson(rate=[1., 10., 100.], name='Three Poissons')
three_poissons
<tfp.distributions.Poisson 'Three Poissons/' batch_shape=(3,) event_shape=() dtype=float32>
three_poissons.log_prob([[1., 10., 100.], [100., 10., 1]])  # sample_shape is [2].
<tf.Tensor: id=594, shape=(2, 3), dtype=float32, numpy=
array([[  -1.       ,   -2.0785608,   -3.222351 ],
       [-364.73938  ,   -2.0785608,  -95.39483  ]], dtype=float32)>
three_poissons.log_prob([[[[1., 10., 100.], [100., 10., 1.]]]])  # sample_shape is [1, 1, 2].
<tf.Tensor: id=601, shape=(1, 1, 2, 3), dtype=float32, numpy=
array([[[[  -1.       ,   -2.0785608,   -3.222351 ],
         [-364.73938  ,   -2.0785608,  -95.39483  ]]]], dtype=float32)>

最初の例でどのように入力と出力が shape [2, 3] を持つか、そして 2 番目の例でそれらが shape [1, 1, 2, 3] を持つか注意してください。

言うべきことはそれが総てでしょう、もしそれがブロードキャスティングのためでないのであれば。ここにブロードキャスティングをひとたび考慮した場合のルールがあります。それを完全に一般化して記述してそしてスカラー分布のための単純化を記します :

  1. n = len(batch_shape) + len(event_shape) と定義します。(スカラー分布については、len(event_shape)=0。)
  2. 入力 tensor t が n 次元よりも少ない場合、それが正確に n 次元を持つまでその shape を size 1 の次元を左に追加してパッドします。結果の tensor を t’ と呼びます。
  3. t’ の n 右端次元を (そのために log_prob を計算している) 分布の [batch_shape, event_shape] に対してブロードキャストします。より詳細には: t’ が既に分布に一致している次元については、何もしません、そして t’ がシングルトンを持つような次元については、そのシングルトンを適切な回数複製 (= replicate) します。他の状況は総てエラーです。(スカラー分布については、batch_shape に対してブロードキャストするだけです、何故ならば event_shape = [] だからです。)
  4. 今は最後に log_prob を計算することができます。結果の tensor は shape [sample_shape, batch_shape] を持ち、そこでは sample_shape は t か t’ の任意の次元として n-右端次元の左側に定義されます : sample_shape = shape(t)[:-n]。(原文: The resulting tensor will have shape [sample_shape, batch_shape], where sample_shape is defined to be any dimensions of t or t’ to the left of the n-rightmost dimensions: sample_shape = shape(t)[:-n].)

それが意味するところを知らないのであればこれは混乱しますので、幾つかの例を動作させましょう :

three_poissons.log_prob([10.])
<tf.Tensor: id=608, shape=(3,), dtype=float32, numpy=array([-16.104412 ,  -2.0785608, -69.052704 ], dtype=float32)>

tensor [10.] (with shape [1]) は 3 の batch_shape に渡りブロードキャストされますので、3 つ総ての Poisson の対数確率を値 10 で評価します。

three_poissons.log_prob([[[1.], [10.]], [[100.], [1000.]]])
<tf.Tensor: id=615, shape=(2, 2, 3), dtype=float32, numpy=
array([[[-1.0000000e+00, -7.6974149e+00, -9.5394829e+01],
        [-1.6104412e+01, -2.0785608e+00, -6.9052704e+01]],

       [[-3.6473938e+02, -1.4348087e+02, -3.2223511e+00],
        [-5.9131279e+03, -3.6195427e+03, -1.4069575e+03]]], dtype=float32)>

上の例では、入力 tensor は shape [2, 2, 1] を持ち、一方で distributions オブジェクトは 3 のバッチ shape を持ちます。従って [2, 2] サンプル次元の各々について、提供される単一の値は 3 の Poisson の各々へのブロードキャストを得ます。

それを考える多分有用な方法は: three_poissons は batch_shape = [2, 3] を持ちますので、log_prob への呼び出しは最後の次元が 1 か 3 である Tensor を取らなければなりません; 他の総てはエラーです。(numpy ブロードキャスティング・ルールはスカラーの特別なケースを shape [1] の Tensor とまったく同値であるとして扱います。)

batch_shape = [2, 3] を持つより複雑な Poisson 分布で遊ぶことにより私達の技術をテストしましょう :

poisson_2_by_3 = tfd.Poisson(
    rate=[[1., 10., 100.,], [2., 20., 200.]],
    name='Two-by-Three Poissons')
poisson_2_by_3.log_prob(1.)
<tf.Tensor: id=623, shape=(2, 3), dtype=float32, numpy=
array([[  -1.       ,   -7.697415 ,  -95.39483  ],
       [  -1.3068528,  -17.004269 , -194.70168  ]], dtype=float32)>
poisson_2_by_3.log_prob([1.])  # Exactly equivalent to above, demonstrating the scalar special case.
<tf.Tensor: id=630, shape=(2, 3), dtype=float32, numpy=
array([[  -1.       ,   -7.697415 ,  -95.39483  ],
       [  -1.3068528,  -17.004269 , -194.70168  ]], dtype=float32)>
poisson_2_by_3.log_prob([[1., 1., 1.], [1., 1., 1.]])  # Another way to write the same thing. No broadcasting.
<tf.Tensor: id=637, shape=(2, 3), dtype=float32, numpy=
array([[  -1.       ,   -7.697415 ,  -95.39483  ],
       [  -1.3068528,  -17.004269 , -194.70168  ]], dtype=float32)>
poisson_2_by_3.log_prob([[1., 10., 100.]])  # Input is [1, 3] broadcast to [2, 3].
<tf.Tensor: id=644, shape=(2, 3), dtype=float32, numpy=
array([[ -1.       ,  -2.0785608,  -3.222351 ],
       [ -1.3068528,  -5.14709  , -33.907654 ]], dtype=float32)>
poisson_2_by_3.log_prob([[1., 10., 100.], [1., 10., 100.]])  # Equivalent to above. No broadcasting.
<tf.Tensor: id=651, shape=(2, 3), dtype=float32, numpy=
array([[ -1.       ,  -2.0785608,  -3.222351 ],
       [ -1.3068528,  -5.14709  , -33.907654 ]], dtype=float32)>
poisson_2_by_3.log_prob([[1., 1., 1.], [2., 2., 2.]])  # No broadcasting.
<tf.Tensor: id=658, shape=(2, 3), dtype=float32, numpy=
array([[  -1.       ,   -7.697415 ,  -95.39483  ],
       [  -1.3068528,  -14.701683 , -190.09651  ]], dtype=float32)>
poisson_2_by_3.log_prob([[1.], [2.]])  # Equivalent to above. Input shape [2, 1] broadcast to [2, 3].
<tf.Tensor: id=665, shape=(2, 3), dtype=float32, numpy=
array([[  -1.       ,   -7.697415 ,  -95.39483  ],
       [  -1.3068528,  -14.701683 , -190.09651  ]], dtype=float32)>

上の例はバッチに渡るブロードキャスティングを伴いましたが、サンプル shape は empty でした。値のコレクションがあり、バッチの各ポイントでの各値の対数確率を得たいとします。それを手動で行えます :

poisson_2_by_3.log_prob([[[1., 1., 1.], [1., 1., 1.]], [[2., 2., 2.], [2., 2., 2.]]])  # Input shape [2, 2, 3].
<tf.Tensor: id=672, shape=(2, 2, 3), dtype=float32, numpy=
array([[[  -1.       ,   -7.697415 ,  -95.39483  ],
        [  -1.3068528,  -17.004269 , -194.70168  ]],

       [[  -1.6931472,   -6.087977 ,  -91.4828   ],
        [  -1.3068528,  -14.701683 , -190.09651  ]]], dtype=float32)>

あるいはブロードキャスティングに最後のバッチ次元を処理させることもできます :

poisson_2_by_3.log_prob([[[1.], [1.]], [[2.], [2.]]])  # Input shape [2, 2, 1].
<tf.Tensor: id=679, shape=(2, 2, 3), dtype=float32, numpy=
array([[[  -1.       ,   -7.697415 ,  -95.39483  ],
        [  -1.3068528,  -17.004269 , -194.70168  ]],

       [[  -1.6931472,   -6.087977 ,  -91.4828   ],
        [  -1.3068528,  -14.701683 , -190.09651  ]]], dtype=float32)>

ブロードキャスティングにちょうど最初のバッチ次元を処理させることもまたできます (多分幾分自然ではないですが) :

poisson_2_by_3.log_prob([[[1., 1., 1.]], [[2., 2., 2.]]])  # Input shape [2, 1, 3].
<tf.Tensor: id=686, shape=(2, 2, 3), dtype=float32, numpy=
array([[[  -1.       ,   -7.697415 ,  -95.39483  ],
        [  -1.3068528,  -17.004269 , -194.70168  ]],

       [[  -1.6931472,   -6.087977 ,  -91.4828   ],
        [  -1.3068528,  -14.701683 , -190.09651  ]]], dtype=float32)>

あるいはブロードキャスティングに両者のバッチ次元を処理させられるでしょう :

poisson_2_by_3.log_prob([[[1.]], [[2.]]])  # Input shape [2, 1, 1].
<tf.Tensor: id=693, shape=(2, 2, 3), dtype=float32, numpy=
array([[[  -1.       ,   -7.697415 ,  -95.39483  ],
        [  -1.3068528,  -17.004269 , -194.70168  ]],

       [[  -1.6931472,   -6.087977 ,  -91.4828   ],
        [  -1.3068528,  -14.701683 , -190.09651  ]]], dtype=float32)>

望む 2 つの値だけを持ったとき上は素晴らしく動作しましたが、総てのバッチポイントで評価することを望む値の長いリストを持つと仮定します。そのためには、shape の右側に size 1 の特別な次元を追加する、次の記法が非常に有用です :

poisson_2_by_3.log_prob(tf.constant([1., 2.])[..., tf.newaxis, tf.newaxis])
<tf.Tensor: id=704, shape=(2, 2, 3), dtype=float32, numpy=
array([[[  -1.       ,   -7.697415 ,  -95.39483  ],
        [  -1.3068528,  -17.004269 , -194.70168  ]],

       [[  -1.6931472,   -6.087977 ,  -91.4828   ],
        [  -1.3068528,  -14.701683 , -190.09651  ]]], dtype=float32)>

これは (知るに値する) strided slice 記法 のインスタンスです。

完全性のために three_poissons に戻ります、同じ例は次のように見えます :

three_poissons.log_prob([[1.], [10.], [50.], [100.]])
<tf.Tensor: id=711, shape=(4, 3), dtype=float32, numpy=
array([[  -1.       ,   -7.697415 ,  -95.39483  ],
       [ -16.104412 ,   -2.0785608,  -69.052704 ],
       [-149.47777  ,  -43.34851  ,  -18.219254 ],
       [-364.73938  , -143.48087  ,   -3.222351 ]], dtype=float32)>
three_poissons.log_prob(tf.constant([1., 10., 50., 100.])[..., tf.newaxis])  # Equivalent to above.
<tf.Tensor: id=722, shape=(4, 3), dtype=float32, numpy=
array([[  -1.       ,   -7.697415 ,  -95.39483  ],
       [ -16.104412 ,   -2.0785608,  -69.052704 ],
       [-149.47777  ,  -43.34851  ,  -18.219254 ],
       [-364.73938  , -143.48087  ,   -3.222351 ]], dtype=float32)>

 

多変量分布

さて多変量分布に向かいます、これは non-empty 事象 shape を持ちます。多変量分布を見てみましょう。

multinomial_distributions = [
    # Multinomial is a vector-valued distribution: if we have k classes,
    # an individual sample from the distribution has k values in it, so the
    # event_shape is `[k]`.
    tfd.Multinomial(total_count=100., probs=[.5, .4, .1],
                    name='One Multinomial'),
    tfd.Multinomial(total_count=[100., 1000.], probs=[.5, .4, .1],
                    name='Two Multinomials Same Probs'),
    tfd.Multinomial(total_count=100., probs=[[.5, .4, .1], [.1, .2, .7]],
                    name='Two Multinomials Same Counts'),
    tfd.Multinomial(total_count=[100., 1000.],
                    probs=[[.5, .4, .1], [.1, .2, .7]],
                    name='Two Multinomials Different Everything')

]

describe_distributions(multinomial_distributions)
tfp.distributions.Multinomial("One Multinomial/", batch_shape=(), event_shape=(3,), dtype=float32)
tfp.distributions.Multinomial("Two Multinomials Same Probs/", batch_shape=(2,), event_shape=(3,), dtype=float32)
tfp.distributions.Multinomial("Two Multinomials Same Counts/", batch_shape=(2,), event_shape=(3,), dtype=float32)
tfp.distributions.Multinomial("Two Multinomials Different Everything/", batch_shape=(2,), event_shape=(3,), dtype=float32)

最後の 3 つの例で、batch_shape は常に [2] ですが、どちらかが共有された total_count か共有された probs を持つ (またはどちらも持たない) ようにブロードキャスティングをどのように使用できるかに注意してください、何故ならばそれらは同じ shape を持つように内部でブロードキャストしています。

既に私達が知っていることを考えれば、サンプリングは straightforward です :

describe_sample_tensor_shapes(multinomial_distributions, sample_shapes)
tfp.distributions.Multinomial("One Multinomial/", batch_shape=(), event_shape=(3,), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 3)
Sample shape: 2
Returned sample tensor shape: (2, 3)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 3)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 3)

tfp.distributions.Multinomial("Two Multinomials Same Probs/", batch_shape=(2,), event_shape=(3,), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 2, 3)
Sample shape: 2
Returned sample tensor shape: (2, 2, 3)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 2, 3)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 2, 3)

tfp.distributions.Multinomial("Two Multinomials Same Counts/", batch_shape=(2,), event_shape=(3,), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 2, 3)
Sample shape: 2
Returned sample tensor shape: (2, 2, 3)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 2, 3)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 2, 3)

tfp.distributions.Multinomial("Two Multinomials Different Everything/", batch_shape=(2,), event_shape=(3,), dtype=float32)
Sample shape: 1
Returned sample tensor shape: (1, 2, 3)
Sample shape: 2
Returned sample tensor shape: (2, 2, 3)
Sample shape: [1, 5]
Returned sample tensor shape: (1, 5, 2, 3)
Sample shape: [3, 4, 5]
Returned sample tensor shape: (3, 4, 5, 2, 3)

対数確率の計算も等しく straightforward です。対角多変量正規分布の例を動かしてみましょう。(多変量はそれほど都合よくはブロードキャストされません、何故ならば counts と probabilities 上の制約はブロードキャスティングがしばしば承認しがたい値を生成することを意味するからです。) 同じ mean しかし異なる scale (標準偏差) の 2 つの 3-次元分布のバッチを使用します :

two_multivariate_normals = tfd.MultivariateNormalDiag(loc=[1., 2., 3.], scale_identity_multiplier=[1., 2.])
two_multivariate_normals
<tfp.distributions.MultivariateNormalDiag 'MultivariateNormalDiag/' batch_shape=(2,) event_shape=(3,) dtype=float32>

(scale が multiples of the identity (単位行列の倍数) であるような分布を使用しましたが、これは制限ではありません; scale_identity_multiplier の代わりに scale を渡せるでしょう。)

さて各バッチポイントの対数確率をその mean と shifted mean で評価しましょう :

two_multivariate_normals.log_prob([[[1., 2., 3.]], [[3., 4., 5.]]])  # Input has shape [2,1,3].
<tf.Tensor: id=2394, shape=(2, 2), dtype=float32, numpy=
array([[-2.7568154, -4.836257 ],
       [-8.756816 , -6.336257 ]], dtype=float32)>

正確に等価に、constant の中 (= middle) に特別な shape=1 次元を挿入するために https://www.tensorflow.org/api_docs/cc/class/tensorflow/ops/strided-slice を使用できます :

two_multivariate_normals.log_prob(
    tf.constant([[1., 2., 3.], [3., 4., 5.]])[:, tf.newaxis, :])  # Equivalent to above.
<tf.Tensor: id=2475, shape=(2, 2), dtype=float32, numpy=
array([[-2.7568154, -4.836257 ],
       [-8.756816 , -6.336257 ]], dtype=float32)>

他方、extra 次元を挿入しないのであれば、[1., 2., 3.] を最初のバッチポイントに渡して [3., 4., 5.] を 2 番目に渡します :

two_multivariate_normals.log_prob(tf.constant([[1., 2., 3.], [3., 4., 5.]]))
<tf.Tensor: id=2552, shape=(2,), dtype=float32, numpy=array([-2.7568154, -6.336257 ], dtype=float32)>

 

Shape 操作テクニック

Reshape Bijector

分布の event_shape を reshape するために Reshape bijector が使用できます。例を見てみましょう :

six_way_multinomial = tfd.Multinomial(total_count=1000., probs=[.3, .25, .2, .15, .08, .02])
six_way_multinomial
<tfp.distributions.Multinomial 'Multinomial/' batch_shape=() event_shape=(6,) dtype=float32>

[6] の事象 shape で多変量分布を作成しました。Reshape Bijector はこれを [2, 3] の事象 shape を持つ分布として扱うことを可能にします。

Bijector は ${\mathbb R}^n$ の開部分集合 (= open subset) の微分可能な、one-to-one 関数を表します。Bijector は TransformedDistribution とともに使用されます、これは 基底分布 $p(x)$ と ($Y = g(X)$ を表わす) Bijector によって分布 $p(y)$ をモデル化します。実際に見てみましょう :

transformed_multinomial = tfd.TransformedDistribution(
    distribution=six_way_multinomial,
    bijector=tfb.Reshape(event_shape_out=[2, 3]))
transformed_multinomial
<tfp.distributions.TransformedDistribution 'reshapeMultinomial/' batch_shape=() event_shape=(2, 3) dtype=float32>
six_way_multinomial.log_prob([500., 100., 100., 150., 100., 50.])
<tf.Tensor: id=2628, shape=(), dtype=float32, numpy=-178.21973>
transformed_multinomial.log_prob([[500., 100., 100.], [150., 100., 50.]])
<tf.Tensor: id=2703, shape=(), dtype=float32, numpy=-178.21973>

これが Reshape bijector のできる唯一のことです: それは事象次元をバッチ次元に変更することはできません or vice-versa。

 

Independent 分布

Independent 分布は独立で、必ずしも同一でない分布 (aka 分布のバッチ)を単一の分布として扱うために使用されます。より簡潔には、batch_shape の次元を event_shape の次元に変換することを可能にします。例で示します :

two_by_five_bernoulli = tfd.Bernoulli(
    probs=[[.05, .1, .15, .2, .25], [.3, .35, .4, .45, .5]],
    name="Two By Five Bernoulli")
two_by_five_bernoulli
<tfp.distributions.Bernoulli 'Two By Five Bernoulli/' batch_shape=(2, 5) event_shape=() dtype=int32>

これを表の関連する確率を持つコインの 2 by 5 配列として考えることができます。1 と 0 の特定の、任意のセットの確率を評価しましょう :

pattern = [[1., 0., 0., 1., 0.], [0., 0., 1., 1., 1.]]
two_by_five_bernoulli.log_prob(pattern)
<tf.Tensor: id=2722, shape=(2, 5), dtype=float32, numpy=
array([[-2.9957323 , -0.10536051, -0.16251893, -1.609438  , -0.28768206],
       [-0.35667494, -0.43078294, -0.9162907 , -0.7985077 , -0.6931472 ]],
      dtype=float32)>

これを 2 つの異なる「5 つの Bernoulli のセット」に変えるために Independent を使用することができます、与えられたパターンで現れるコイン投げの「行 (= row)」を単一の結果として考えたい場合にこれは有用です :

two_sets_of_five = tfd.Independent(
    distribution=two_by_five_bernoulli,
    reinterpreted_batch_ndims=1,
    name="Two Sets Of Five")
two_sets_of_five
<tfp.distributions.Independent 'Two Sets Of Five/' batch_shape=(2,) event_shape=(5,) dtype=int32>

数学的には、5 つの各「セット」の対数確率をセットの 5 つの「独立した」コイン投げの対数確率を合計することで計算しています、それは分布がその名前を得るところです :

two_sets_of_five.log_prob(pattern)
<tf.Tensor: id=2739, shape=(2,), dtype=float32, numpy=array([-5.160732 , -3.1954036], dtype=float32)>

さらに進んでそして (そこでは個々の事象は 2-by-5 Bernoulli 分布のセットであるような) 分布を作成するために Independent を使用できます :

one_set_of_two_by_five = tfd.Independent(
    distribution=two_by_five_bernoulli, reinterpreted_batch_ndims=2,
    name="One Set Of Two By Five")
one_set_of_two_by_five.log_prob(pattern)
<tf.Tensor: id=2756, shape=(), dtype=float32, numpy=-8.356134>

注目すべき点はサンプルの観点から、Independent の使用は何も変更しないことです :

describe_sample_tensor_shapes(
    [two_by_five_bernoulli,
     two_sets_of_five,
     one_set_of_two_by_five],
    [[3, 5]])
tfp.distributions.Bernoulli("Two By Five Bernoulli/", batch_shape=(2, 5), event_shape=(), dtype=int32)
Sample shape: [3, 5]
Returned sample tensor shape: (3, 5, 2, 5)

tfp.distributions.Independent("Two Sets Of Five/", batch_shape=(2,), event_shape=(5,), dtype=int32)
Sample shape: [3, 5]
Returned sample tensor shape: (3, 5, 2, 5)

tfp.distributions.Independent("One Set Of Two By Five/", batch_shape=(), event_shape=(2, 5), dtype=int32)
Sample shape: [3, 5]
Returned sample tensor shape: (3, 5, 2, 5)

 
課題: As a parting exercise for the reader, we suggest considering the differences and similarities between a vector batch of Normal distributions and a MultivariateNormalDiag distribution from a sampling and log probability perspective. How can we use Independent to construct a MultivariateNormalDiag from a batch of Normals? (Note that MultivariateNormalDiag is not actually implemented this way.)

 

以上






TensorFlow Probability : Tutorials : 分布イントロダクション

TensorFlow Probability : Tutorials : 分布イントロダクション (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 10/24/2018

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです:

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : 分布イントロダクション

この notebook では、TensorFlow 分布 (TFD for short) を調査します。この notebook の目的は tensor shape の TFD の処理を理解することを含めて、優しく貴方に学習カーブを上げて頂くことです。この notebook は抽象的な概念の前に例を提示しようとします。最初に物事を進めるための標準的で簡単な方法を提示して、最も一般的な抽象的な眺めは最後まで取っておきます。もし貴方がより抽象的でリファレンス・スタイルのチュートリアルを好みタイプであれば、Understanding TensorFlow Distributions Shapes をチェックしてください。ここでのマテリアルについてどのような疑問でもあれば、TensorFlow Probability mailing list にコンタクト (or 参加) することを躊躇しないでください。We’re happy to help.

始める前に、適切なライブラリをインポートする必要があります。私達の全体的なライブラリは tensorflow_probability です。慣習により、一般的に distributions ライブラリを tfd として参照します。

Tensorflow Eager は TensorFlow のための命令型実行環境です。TensorFlow eager では、総ての TF 演算は直ちに評価されて結果を生成します。これは TensorFlow の標準的な「グラフ」モードとは対称的です、そこでは TF 演算はグラフに (後で実行される) ノードを追加します。ここで表わされる概念のどれもそれに依存はしてませんが、この notebook 全体は TF Eager を使用して書かれます、そして TFP はグラフモードでも利用できます。

!pip install -q tensorflow-probability
import collections

import tensorflow as tf
import tensorflow_probability as tfp
tfd = tfp.distributions

tfe = tf.contrib.eager
try:
  tfe.enable_eager_execution()
except ValueError:
  pass

import matplotlib.pyplot as plt

from __future__ import print_function

 

基本的な単変量分布

Let’s dive right in そして正規分布を作成しましょう :

n = tfd.Normal(loc=0., scale=1.)
n
<tfp.distributions.Normal 'Normal/' batch_shape=() event_shape=() dtype=float32>

それからサンプルをドローできます :

n.sample()
<tf.Tensor: id=23, shape=(), dtype=float32, numpy=-1.2775375>

複数のサンプルをドローできます :

n.sample(3)
<tf.Tensor: id=44, shape=(3,), dtype=float32, numpy=array([ 0.49041903,  1.3672128 , -0.39134106], dtype=float32)>

log prob を評価できます :

n.log_prob(0.)
<tf.Tensor: id=55, shape=(), dtype=float32, numpy=-0.9189385>

複数の対数確率を評価できます :

n.log_prob([0., 2., 4.])
<tf.Tensor: id=64, shape=(3,), dtype=float32, numpy=array([-0.9189385, -2.9189386, -8.918939 ], dtype=float32)>

広範囲の分布を持ちます。Bernoulli を試してみましょう:

b = tfd.Bernoulli(probs=0.7)
b
<tfp.distributions.Bernoulli 'Bernoulli/' batch_shape=() event_shape=() dtype=int32>
b.sample()
<tf.Tensor: id=92, shape=(), dtype=int32, numpy=1>
b.sample(8)
<tf.Tensor: id=114, shape=(8,), dtype=int32, numpy=array([0, 1, 1, 1, 1, 1, 0, 1], dtype=int32)>
b.log_prob(1)
<tf.Tensor: id=129, shape=(), dtype=float32, numpy=-0.35667494>
b.log_prob([1, 0, 1, 0])
<tf.Tensor: id=149, shape=(4,), dtype=float32, numpy=array([-0.35667494, -1.2039728 , -0.35667494, -1.2039728 ], dtype=float32)>

 

多変量分布

対角共分散 (= diagonal covariance) で多変量正規分布を作成します:

nd = tfd.MultivariateNormalDiag(loc=[0., 10.], scale_diag=[1., 4.])
nd
<tfp.distributions.MultivariateNormalDiag 'MultivariateNormalDiag/' batch_shape=() event_shape=(2,) dtype=float32>

これを先に作成した単変量正規分布と比較します、何が違うのでしょう?

tfd.Normal(loc=0., scale=1.)
<tfp.distributions.Normal 'Normal/' batch_shape=() event_shape=() dtype=float32>

単変量正規分布は () の event_shape を持つことを見ます、これはそれがスカラー分布であることを示します。多変量正規分布は 2 の event_shape を持ち、これはこの分布の基本的な事象空間 (= event space) は 2-次元であることを示しています。

サンプリングはちょうど前のように動作します :

nd.sample()
<tf.Tensor: id=317, shape=(2,), dtype=float32, numpy=array([-0.43942595,  9.620465  ], dtype=float32)>
nd.sample(5)
<tf.Tensor: id=425, shape=(5, 2), dtype=float32, numpy=
array([[-0.10185293,  9.465504  ],
       [ 1.3469841 , 11.134725  ],
       [ 3.049023  , -3.3670492 ],
       [ 0.26367915, 18.198818  ],
       [-0.8832267 , 11.224263  ]], dtype=float32)>
nd.log_prob([0., 10])
<tf.Tensor: id=526, shape=(), dtype=float32, numpy=-3.2241714>

多変量正規分布は一般には対角共分散を持ちません。TFD は多変量戦記分布を作成するための複数の方法を提供し、フル共分散仕様を含みます。それはここで利用します。

nd = tfd.MultivariateNormalFullCovariance(
    loc = [0., 5], covariance_matrix = [[1., .7], [.7, 1.]])
data = nd.sample(200)
plt.scatter(data[:, 0], data[:, 1], color='blue', alpha=0.4)
plt.axis([-5, 5, 0, 10])
plt.title("Data set")
plt.show()

 

複数の分布

最初の Bernoulli 分布は一枚の公正なコインのフリップを表しました。独立な Bernoulli 分布のバッチもまた作成できます、単一の Distribution オブジェクトで、それぞれはそれら自身のパラメータを伴います:

b3 = tfd.Bernoulli(probs=[.3, .5, .7])
b3
<tfp.distributions.Bernoulli 'Bernoulli/' batch_shape=(3,) event_shape=() dtype=int32>

これが何を意味するかを明らかにすることは重要です。上の呼び出しは 3 つの独立な Bernoulli 分布を定義し、これらは偶々同じ Python Distribution オブジェクトに含まれています。3 つの分布は個々に操作はできません。batch_shape がどのように (3,) であるかに注意してください、これは 3 つの分布のバッチを示します、そして event_shape は () です、これは個々の分布は単変量事象空間を持つことを示します。

sample を呼びだせば、3 つ総てからサンプルを得ます :

b3.sample()
<tf.Tensor: id=744, shape=(3,), dtype=int32, numpy=array([1, 0, 1], dtype=int32)>
b3.sample(6)
<tf.Tensor: id=766, shape=(6, 3), dtype=int32, numpy=
array([[0, 0, 1],
       [0, 1, 1],
       [1, 0, 1],
       [0, 0, 1],
       [0, 1, 1],
       [0, 0, 1]], dtype=int32)>

prob を呼び出す場合、(これは log_prob と同じ shape セマンティクスを持ちます; 明確にするためにこれらの小さな Bernoulli サンプルで prob を使用します、log_prob は通常はアプリケーションで好まれますが) それにベクトルを渡してその値を生成する各コインの確率を評価できます :

b3.prob([1, 1, 0])
<tf.Tensor: id=782, shape=(3,), dtype=float32, numpy=array([0.29999998, 0.5       , 0.29999998], dtype=float32)>

何故 API はバッチ shape を含むのでしょう?意味的には、分布のリストを作成して for ループでそれらに渡り反復することで同じ計算を遂行できるでしょう (少なくとも Eager モードで、TF graph モードでは tf.while ループを使用必要があるでしょう)。けれども、同値にパラメータ化された分布の (潜在的に巨大な) セットを持つことは極めて一般的で、可能なときはいつでもベクトル化された計算の使用はハードウェア・アクセラレータを使用して高速な計算を遂行することを可能にする点において主要な構成要素です。

 

バッチを事象に集めるために Independent を使用する

前のセクションでは、b3 を作成しました、3 つのコイン投げを表わす単一の Distribution オブジェクトです。ベクトル $v$ で b3.prob を呼び出した場合、$i$ 番目のエントリは $i$ 番目のコインが値 $v[i]$ を取る確率でした。

代わりに、同じ基礎的な族 (= family) から独立な確率変数に渡る「同時」分布を指定したいと仮定します。これは数学的には異なるオブジェクトで、そこではこの新しい分布に対して、ベクトル $v$ 上の prob は (コインの全体セットがベクトル $v$ に一致する) 確率を表わす 単一の値 を返します。

これをどのように成し遂げるのでしょう?Independent と呼ばれる「高階」分布を使用します、これは分布を取り (batch shape が event shape に移された) 新しい分布を生成します :

b3_joint = tfd.Independent(b3, reinterpreted_batch_ndims=1)
b3_joint
<tfp.distributions.Independent 'IndependentBernoulli/' batch_shape=() event_shape=(3,) dtype=int32>

shape を元の b3 のそれと比較します :

b3 
<tfp.distributions.Bernoulli 'Bernoulli/' batch_shape=(3,) event_shape=() dtype=int32>

約束したように、その Independent は batch shape を event shape に移したことを見ます : b3_joint は 3-次元事象空間 (event_shape = (3,)) に渡る単一の分布 (batch_shape = ()) です

セマンティクスをチェックしてみましょう :

b3_joint.prob([1, 1, 0])
<tf.Tensor: id=801, shape=(), dtype=float32, numpy=0.044999998>

同じ結果を得るための代わりの方法は b3 を使用して確率を計算して乗算により手動で reduction を行なうことです (あるいは、対数確率が使用されるより通常のケースでは、合計します) :

tf.reduce_prod(b3.prob([1, 1, 0]))
<tf.Tensor: id=818, shape=(), dtype=float32, numpy=0.044999994>

Indpendent は望まれる概念をより明示的に表わすことをユーザに可能にします。これを極めて有用であるとみなします、それは厳密には必要ではありませんが。

Fun facts:

  • b3.sample と b3_joint.sample は異なる概念の実装を持ちますが、見分けがつかない出力を持ちます: 独立な分布のバッチとIndependent を使用したバッチから作成された単一の分布の間の違いはサンプリングのときではなく、確率を計算するときに現出します。
  • MultivariateNormalDiag はスカラー Nomal と Independent 分布を使用して自明に実装されるでしょう (それは実際にはこのように実装されていませんが、可能でしょう)。

 

多変量分布のバッチ

3 つの full-共分散 2-次元多変量正規分布のバッチを作成します :

nd_batch = tfd.MultivariateNormalFullCovariance(
    loc = [[0., 0.], [1., 1.], [2., 2.]],
    covariance_matrix = [[[1., .1], [.1, 1.]], 
                         [[1., .3], [.3, 1.]],
                         [[1., .5], [.5, 1.]]])
nd_batch
<tfp.distributions.MultivariateNormalFullCovariance 'MultivariateNormalFullCovariance/' batch_shape=(3,) event_shape=(2,) dtype=float32>

batch_shape = (3,) を見ますので、3 つの独立な多変量正規分布があり、そして event_shape = (2,) ですので、各多変量正規分布は 2-次元です。この例では、個々の分布は独立な要素を持ちません。

サンプリングが動作します :

nd_batch.sample(4)
<tf.Tensor: id=992, shape=(4, 3, 2), dtype=float32, numpy=
array([[[ 0.8355589 ,  1.0287906 ],
        [ 0.6796141 ,  1.7522928 ],
        [ 0.603493  ,  0.13745153]],

       [[-0.49636343, -1.3271412 ],
        [ 0.89452106,  0.6159041 ],
        [ 0.7503762 ,  1.1700599 ]],

       [[-0.3203717 ,  0.99095964],
        [ 0.22984672, -1.4632237 ],
        [ 2.1868286 ,  2.476635  ]],

       [[ 1.0310459 , -0.35648566],
        [ 1.6214066 ,  0.5076788 ],
        [ 1.9965721 ,  2.3196375 ]]], dtype=float32)>

batch_shape = (3,) and event_shape = (2,) ですから、log_prob に shape (3, 2) の tensor を渡します :

nd_batch.log_prob([[0., 0.], [1., 1.], [2., 2.]])
<tf.Tensor: id=1090, shape=(3,), dtype=float32, numpy=array([-1.8328519, -1.7907217, -1.694036 ], dtype=float32)>

 

ブロードキャスティング, aka 何故これはそれほど混乱させるのでしょう?

ここまで行なったことから抽出すると、総ての分布は batch shape B と event shape E を持ちます。BE を event shapes の連結とします (訳注: 原文ママ) :

  • 単変量スカラー分布 n と b について、BE = ()..
  • 2-次元多変量正規分布 nd について、BE = (2)。
  • b3 と b3_joint の両者について、BE = (3)。
  • 多変量正規分布 ndb のバッチについて、BE = (3, 2)。

ここまで使用した「評価ルール」は :

  • 引数なしの sample は shape BE の tensor を返します; スカラー n を伴うサンプリングは “n by BE” tensor を返します。
  • prob と log_prob は shape BE の tensor を取り shape B の結果を返します。

prob と log_prob のための実際の「評価ルール」はより複雑で、潜在的なパワーとスピードを、しかし複雑さと挑戦も提供する方法でです。実際のルールは (本質的に) log_prob への引数は BE に対して ブロードキャスト可能 です; 任意の「特別な (= extra)」次元は出力で保持されます

含蓄を調べてみましょう。単変量正規分布 n, BE = () について、log_prob はスカラーを想定します。log_prob に non-empty shape を持つ tensor を渡す場合、それらは出力のバッチ次元として出現します:

n = tfd.Normal(loc=0., scale=1.)
n
<tfp.distributions.Normal 'Normal/' batch_shape=() event_shape=() dtype=float32>
n.log_prob(0.)
<tf.Tensor: id=1099, shape=(), dtype=float32, numpy=-0.9189385>
n.log_prob([0.])
<tf.Tensor: id=1108, shape=(1,), dtype=float32, numpy=array([-0.9189385], dtype=float32)>
n.log_prob([[0., 1.], [-1., 2.]])
<tf.Tensor: id=1117, shape=(2, 2), dtype=float32, numpy=
array([[-0.9189385, -1.4189385],
       [-1.4189385, -2.9189386]], dtype=float32)>

2-次元多変量正規分布 nd へ向かいましょう (例示目的でパラメータは変更されています) :

nd = tfd.MultivariateNormalDiag(loc=[0., 1.], scale_diag=[1., 1.])
nd
<tfp.distributions.MultivariateNormalDiag 'MultivariateNormalDiag/' batch_shape=() event_shape=(2,) dtype=float32>

log_prob は shape (2,) を持つ引数を「期待」しますが、それはこの shape に対してブロードキャストする任意の引数を受け取るでしょう :

nd.log_prob([0., 0.])
<tf.Tensor: id=1272, shape=(), dtype=float32, numpy=-2.337877>

しかし「より多くの」サンプルを渡して、そしてそれらの log_prob 総てを一度に評価することもできます :

nd.log_prob([[0., 0.],
             [1., 1.],
             [2., 2.]])
<tf.Tensor: id=1352, shape=(3,), dtype=float32, numpy=array([-2.337877 , -2.337877 , -4.3378773], dtype=float32)>

多分魅力は少ないですが、事象次元に渡りブロードキャストすることもできます :

nd.log_prob([0.])
<tf.Tensor: id=1432, shape=(), dtype=float32, numpy=-2.337877>
nd.log_prob([[0.], [1.], [2.]])
<tf.Tensor: id=1512, shape=(3,), dtype=float32, numpy=array([-2.337877 , -2.337877 , -4.3378773], dtype=float32)>

このようにブロードキャストすることは私達の “enable broadcasting whenever possible” デザインの結果です; この使用方法は幾分議論の的になりそして TFP の将来バージョンで潜在的に除去されるかもしれません。

さて 3 つのコインの例を再度見てみましょう :

b3 = tfd.Bernoulli(probs=[.3, .5, .7])

ここで、各コインが表が出る確率を表わすためのブロードキャスティングの使用は非常に直感的です :

b3.prob([1])
<tf.Tensor: id=1539, shape=(3,), dtype=float32, numpy=array([0.29999998, 0.5       , 0.7       ], dtype=float32)>

(これを b3.prob([1., 1., 1.]) と比較してください、これは b3 が導入されたところで使用されたでしょう。)

さて各コインについて、コインが表になる確率と裏になる確率を知ることを望むとします。次を試すことを想像するかもしれません :

b3.log_prob([0, 1])

不幸なことに、これは 長くて not-very-readable スタックトレースを伴うエラーを生成します。b3 は BE = (3) を持ちますので、b3.prob には (3,) に対してブロードキャスト可能な何かを渡さなければなりません。[0, 1] は shape (2) を持ちますので、それはブロードキャストしないでエラーを作成します。代わりに、次のように言わなければなりません :

b3.prob([[0], [1]])
<tf.Tensor: id=1560, shape=(2, 3), dtype=float32, numpy=
array([[0.7       , 0.5       , 0.29999998],
       [0.29999998, 0.5       , 0.7       ]], dtype=float32)>

Why? [[0], [1]] は shape (2, 1) を持ちますので、(2, 3) のブロードキャスト shape を作成するために shape (3) に対してブロードキャストします。

ブロードキャスティングは非常にパワフルです: それが使用されるメモリの総量において 1 桁削減を可能にするケースがあるでしょう、そしてそれはしばしばユーザ・コードをより短くします。けれども、それでプログラムすることは挑戦的である可能性があります。log_prob を呼び出してエラーを得る場合、ブロードキャストの失敗は殆どいつも問題かもしれません。

 

Going Farther

このチュートリアルでは、単純な紹介を (願わくば) 提供しました。更に先に進むための幾つかの指針は :

  • event_shape, batch_shape と sample_shape は任意のランクであり得ます (このチュートリアルではそれらは常にスカラーかランク 1 です)。これはパワーを増大しますが再度プログラミングの挑戦に繋がるかもしれません、特にブロードキャスティングを伴うときには。shape 操作を更にに深く掘り下げて調べるためには、Understanding TensorFlow Distributions Shapes を見てください。
  • TFP は Bijectors として知られるパワフルな抽象を含みます、これは TransformedDistribution と連携して (既存の分布の可逆な (= invertible) 変換である) 新しい分布を簡単に作成するための柔軟で、合成的な方法を与えます。私達はこれについてのチュートリアルを間もなく書こうとしますが、当面は、ドキュメント をチェックしてください。
 

以上






TensorFlow Probability : 概要

TensorFlow Probability : 概要 (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 10/23/2018

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Overview を翻訳した上で
適宜、補足説明したものです:


* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability 概要

TensorFlow Probability は TensorFlow における確率的推論と統計解析のためのライブラリです。TensorFlow エコシステムの一部として、TensorFlow Probability は自動微分を使用した深層ネットワーク・勾配ベースの推論と確率的方法の統合、巨大なデータセットへのスケーラビリティ、ハードウェア・アクセラレーション (GPU) を伴うモデルそして分散計算を提供します。

 

コンポーネント

私達の確率的機械学習ツールは次のように構造化されています :

層 0: TensorFlow

数値演算 — 特に、LinearOperator クラス — は効率的な計算のための特定の構造 (対角, 低ランク, etc.) を開発できる行列-free な実装を可能にします。それは TensorFlow Probability チームにより構築されて維持されてコア TensorFlow の tf.linalg の一部です。

 

層 1: 統計的ビルディング・ブロック

  • Distributions (tfp.distributions, tf.distributions): 確率分布とバッチとブロードキャスティング・セマンティクスを持つ関連する統計の巨大なコレクション。
  • Bijectors (tfp.bijectors): 確率変数のリバーシブルで構成可能な変換。Bijectors は対数正規分布のような古典的な例から masked autoregressive flows のような洗練された深層学習モデルまで、変換された分布のリッチなクラスを提供します、

 

層 2: モデル構築

  • Edward2 (tfp.edward2): 柔軟な確率的モデルをプログラムとして指定するための確率的プログラミング言語。
  • 確率的層 (= Probabilistic layers) (tfp.layers): それらが表わす関数に渡る不確かさを持つニューラルネットワーク層で、TensorFlow 層を拡張しています。
  • 訓練可能な分布 (tfp.trainable_distributions): 単一の tensor でパラメータ化された確率分布で、確率分布を出力するニューラルネットを構築することを容易にします。

 

層 3: 確率的推論

  • マルコフ連鎖モンテカルロ (tfp.mcmc): サンプリングを通して積分を近似するためのアルゴリズム。ハミルトニアンモンテカルロ、ランダム・ウォーク・メトロポリス・ヘイスティングス、そしてカスタム推移カーネルを構築するための機能を含みます。
  • 変分推論 (tfp.vi): 最適化を通して積分を近似するためのアルゴリズム。
  • オプティマイザ (tfp.optimizer): TensorFlow Optimizers を拡張した確率的最適化メソッド。確率的勾配 Langevin Dynamics を含みます。
  • モンテカルロ (tfp.monte_carlo): モンテカルロ期待値を計算するためのツール。

 

サンプル

ナビゲーションでリストされている Python notebook チュートリアルに加えて、利用可能な幾つかのサンプル・スクリプトがあります :

 

以上






AI導入支援 #2 ウェビナー

スモールスタートを可能としたAI導入支援   Vol.2
[無料 WEB セミナー] [詳細]
「画像認識 AI PoC スターターパック」の紹介
既に AI 技術を実ビジネスで活用し、成果を上げている日本企業も多く存在しており、競争優位なビジネスを展開しております。
しかしながら AI を導入したくとも PoC (概念実証) だけでも高額な費用がかかり取組めていない企業も少なくないようです。A I導入時には欠かせない PoC を手軽にしかも短期間で認知度を確認可能とするサービの紹介と共に、AI 技術の特性と具体的な導入プロセスに加え運用時のポイントについても解説いたします。
日時:2021年10月13日(水)
会場:WEBセミナー
共催:クラスキャット、日本FLOW(株)
後援:働き方改革推進コンソーシアム
参加費: 無料 (事前登録制)
人工知能開発支援
◆ クラスキャットは 人工知能研究開発支援 サービスを提供しています :
  • テクニカルコンサルティングサービス
  • 実証実験 (プロトタイプ構築)
  • アプリケーションへの実装
  • 人工知能研修サービス
◆ お問合せ先 ◆
(株)クラスキャット
セールス・インフォメーション
E-Mail:sales-info@classcat.com