ホーム » 「TensorFlow Probability 0.10 Guide」タグがついた投稿
タグアーカイブ: TensorFlow Probability 0.10 Guide
TensorFlow Probability 0.10 : ガイド : Optimizer
TensorFlow Probability 0.10 : ガイド : Optimizer (翻訳/解説)
翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 06/19/2020 (0.10.0)
* 本ページは、TensorFlow Probability の以下のドキュメントを翻訳した上で適宜、補足説明したものです：
* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。
ガイド : Optimizer
Abstract
この colab では TensorFlow Probability で実装された様々な optimizer をどのように使用するかを実演します。
依存性 & 要件
#@title Import { display-mode: "form" }
%matplotlib inline
import contextlib
import functools
import os
import time
import numpy as np
import pandas as pd
import scipy as sp
from six.moves import urllib
from sklearn import preprocessing
import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
import tensorflow_probability as tfp
BFGS と L-BFGS Optimizers
準ニュートン法はポピュラーな first order 最適化アルゴリズムのクラスです。これらのメソッドは探索方向を見つけるために正確なヘッセ (行列) への正定値近似を利用します。
Broyden-Fletcher-Goldfarb-Shanno アルゴリズム (BFGS) はこの一般的なアイデアの特定の実装です。それはミディアムサイズの問題のために適用可能で選択される方法です、そこでは勾配はどこでも連続です (e.g. $L_2$ ペナルティを持つ線形回帰)。
L-BFGS は BFGS の制限されたメモリバージョンで、そのヘッセ行列が合理的なコストで計算できないかスパースでないようなより大きな問題を解くために有用です。ヘッセ行列の完全密 $n \times n$ 近似をストアする代わりに、それらはこれらの近似を暗黙的に表す長さ $n$ の 2, 3 のベクトルをセーブするだけです。
#@title Helper functions
CACHE_DIR = os.path.join(os.sep, 'tmp', 'datasets')
def make_val_and_grad_fn(value_fn):
@functools.wraps(value_fn)
def val_and_grad(x):
return tfp.math.value_and_gradient(value_fn, x)
return val_and_grad
@contextlib.contextmanager
def timed_execution():
t0 = time.time()
yield
dt = time.time() - t0
print('Evaluation took: %f seconds' % dt)
def np_value(tensor):
"""Get numpy value out of possibly nested tuple of tensors."""
if isinstance(tensor, tuple):
return type(tensor)(*(np_value(t) for t in tensor))
else:
return tensor.numpy()
def run(optimizer):
"""Run an optimizer and measure it's evaluation time."""
optimizer() # Warmup.
with timed_execution():
result = optimizer()
return np_value(result)
単純な二次関数上の L-BFGS
# Fix numpy seed for reproducibility
np.random.seed(12345)
# The objective must be supplied as a function that takes a single
# (Tensor) argument and returns a tuple. The first component of the
# tuple is the value of the objective at the supplied point and the
# second value is the gradient at the supplied point. The value must
# be a scalar and the gradient must have the same shape as the
# supplied argument.
# The `make_val_and_grad_fn` decorator helps transforming a function
# returning the objective value into one that returns both the gradient
# and the value. It also works for both eager and graph mode.
dim = 10
minimum = np.ones([dim])
scales = np.exp(np.random.randn(dim))
@make_val_and_grad_fn
def quadratic(x):
return tf.reduce_sum(scales * (x - minimum) ** 2, axis=-1)
# The minimization routine also requires you to supply an initial
# starting point for the search. For this example we choose a random
# starting point.
start = np.random.randn(dim)
# Finally an optional argument called tolerance let's you choose the
# stopping point of the search. The tolerance specifies the maximum
# (supremum) norm of the gradient vector at which the algorithm terminates.
# If you don't have a specific need for higher or lower accuracy, leaving
# this parameter unspecified (and hence using the default value of 1e-8)
# should be good enough.
tolerance = 1e-10
@tf.function
def quadratic_with_lbfgs():
return tfp.optimizer.lbfgs_minimize(
quadratic,
initial_position=tf.constant(start),
tolerance=tolerance)
results = run(quadratic_with_lbfgs)
# The optimization results contain multiple pieces of information. The most
# important fields are: 'converged' and 'position'.
# Converged is a boolean scalar tensor. As the name implies, it indicates
# whether the norm of the gradient at the final point was within tolerance.
# Position is the location of the minimum found. It is important to check
# that converged is True before using the value of the position.
print('L-BFGS Results')
print('Converged:', results.converged)
print('Location of the minimum:', results.position)
print('Number of iterations:', results.num_iterations)
Evaluation took: 0.014586 seconds L-BFGS Results Converged: True Location of the minimum: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] Number of iterations: 10
BFGS で同じ問題
@tf.function
def quadratic_with_bfgs():
return tfp.optimizer.bfgs_minimize(
quadratic,
initial_position=tf.constant(start),
tolerance=tolerance)
results = run(quadratic_with_bfgs)
print('BFGS Results')
print('Converged:', results.converged)
print('Location of the minimum:', results.position)
print('Number of iterations:', results.num_iterations)
Evaluation took: 0.010468 seconds BFGS Results Converged: True Location of the minimum: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] Number of iterations: 10
L1 ペナルティを持つ線形回帰 : 前立腺がん
次の本からの例 : The Elements of Statistical Learning, Data Mining, Inference, and Prediction by Trevor Hastie, Robert Tibshirani and Jerome Friedman.
これは L1 ペナルティを持つ最適化問題であることに注意してください。
データセットを取得する
def cache_or_download_file(cache_dir, url_base, filename):
"""Read a cached file or download it."""
filepath = os.path.join(cache_dir, filename)
if tf.io.gfile.exists(filepath):
return filepath
if not tf.io.gfile.exists(cache_dir):
tf.io.gfile.makedirs(cache_dir)
url = url_base + filename
print("Downloading {url} to {filepath}.".format(url=url, filepath=filepath))
urllib.request.urlretrieve(url, filepath)
return filepath
def get_prostate_dataset(cache_dir=CACHE_DIR):
"""Download the prostate dataset and read as Pandas dataframe."""
url_base = 'http://web.stanford.edu/~hastie/ElemStatLearn/datasets/'
return pd.read_csv(
cache_or_download_file(cache_dir, url_base, 'prostate.data'),
delim_whitespace=True, index_col=0)
prostate_df = get_prostate_dataset()
Downloading http://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data to /tmp/datasets/prostate.data.
問題定義
np.random.seed(12345)
feature_names = ['lcavol', 'lweight', 'age', 'lbph', 'svi', 'lcp',
'gleason', 'pgg45']
# Normalize features
scalar = preprocessing.StandardScaler()
prostate_df[feature_names] = pd.DataFrame(
scalar.fit_transform(
prostate_df[feature_names].astype('float64')))
# select training set
prostate_df_train = prostate_df[prostate_df.train == 'T']
# Select features and labels
features = prostate_df_train[feature_names]
labels = prostate_df_train[['lpsa']]
# Create tensors
feat = tf.constant(features.values, dtype=tf.float64)
lab = tf.constant(labels.values, dtype=tf.float64)
dtype = feat.dtype
regularization = 0 # regularization parameter
dim = 8 # number of features
# We pick a random starting point for the search
start = np.random.randn(dim + 1)
def regression_loss(params):
"""Compute loss for linear regression model with L1 penalty
Args:
params: A real tensor of shape [dim + 1]. The zeroth component
is the intercept term and the rest of the components are the
beta coefficients.
Returns:
The mean square error loss including L1 penalty.
"""
params = tf.squeeze(params)
intercept, beta = params[0], params[1:]
pred = tf.matmul(feat, tf.expand_dims(beta, axis=-1)) + intercept
mse_loss = tf.reduce_sum(
tf.cast(
tf.losses.mean_squared_error(y_true=lab, y_pred=pred), tf.float64))
l1_penalty = regularization * tf.reduce_sum(tf.abs(beta))
total_loss = mse_loss + l1_penalty
return total_loss
L-BFGS で解く
L-BFGS を使用して fit します。L1 ペナルティは導関数不連続性を導入しさえしますが、実際には、L-BFGS は依然として非常に上手く動作します。
@tf.function
def l1_regression_with_lbfgs():
return tfp.optimizer.lbfgs_minimize(
make_val_and_grad_fn(regression_loss),
initial_position=tf.constant(start),
tolerance=1e-8)
results = run(l1_regression_with_lbfgs)
minimum = results.position
fitted_intercept = minimum[0]
fitted_beta = minimum[1:]
print('L-BFGS Results')
print('Converged:', results.converged)
print('Intercept: Fitted ({})'.format(fitted_intercept))
print('Beta: Fitted {}'.format(fitted_beta))
Evaluation took: 0.017987 seconds L-BFGS Results Converged: True Intercept: Fitted (2.3879985744556484) Beta: Fitted [ 0.68626215 0.28193532 -0.17030254 0.10799274 0.33634988 -0.24888523 0.11992237 0.08689026]
Nelder Mead で解く
Nelder Mead 法 は最もポピュラーな導関数が不要な最小化法の一つです。この optimizer は勾配情報を利用しません、そしてターゲット関数の微分可能性について仮定しません ; 従ってそれは非滑らかな目的関数のために適切です、例えば L1 ペナルティを持つ最適化問題です。
$n$-次元の最適化問題についてそれは非退化単体に広がる (= span) $n+1$ 解のセットを維持します。それは続いて頂点の各々における関数値を使用して移動のセット (反射 (reflection)、膨張 (expansion)、縮小 (shrinkage) そして 収縮 (contraction)) に基づいて単体を修正します。
# Nelder mead expects an initial_vertex of shape [n + 1, 1].
initial_vertex = tf.expand_dims(tf.constant(start, dtype=dtype), axis=-1)
@tf.function
def l1_regression_with_nelder_mead():
return tfp.optimizer.nelder_mead_minimize(
regression_loss,
initial_vertex=initial_vertex,
func_tolerance=1e-10,
position_tolerance=1e-10)
results = run(l1_regression_with_nelder_mead)
minimum = results.position.reshape([-1])
fitted_intercept = minimum[0]
fitted_beta = minimum[1:]
print('Nelder Mead Results')
print('Converged:', results.converged)
print('Intercept: Fitted ({})'.format(fitted_intercept))
print('Beta: Fitted {}'.format(fitted_beta))
Evaluation took: 0.325643 seconds Nelder Mead Results Converged: True Intercept: Fitted (2.387998456121595) Beta: Fitted [ 0.68626266 0.28193456 -0.17030291 0.10799375 0.33635132 -0.24888703 0.11992244 0.08689023]
L2 ペナルティを持つロジスティック回帰
この例のために、分類のための合成データセットを作成してパラメータを fit するために L-BFGS optimizer を使用します。
np.random.seed(12345)
dim = 5 # The number of features
n_obs = 10000 # The number of observations
betas = np.random.randn(dim) # The true beta
intercept = np.random.randn() # The true intercept
features = np.random.randn(n_obs, dim) # The feature matrix
probs = sp.special.expit(
np.matmul(features, np.expand_dims(betas, -1)) + intercept)
labels = sp.stats.bernoulli.rvs(probs) # The true labels
regularization = 0.8
feat = tf.constant(features)
lab = tf.constant(labels, dtype=feat.dtype)
@make_val_and_grad_fn
def negative_log_likelihood(params):
"""Negative log likelihood for logistic model with L2 penalty
Args:
params: A real tensor of shape [dim + 1]. The zeroth component
is the intercept term and the rest of the components are the
beta coefficients.
Returns:
The negative log likelihood plus the penalty term.
"""
intercept, beta = params[0], params[1:]
logit = tf.matmul(feat, tf.expand_dims(beta, -1)) + intercept
log_likelihood = tf.reduce_sum(tf.nn.sigmoid_cross_entropy_with_logits(
labels=lab, logits=logit))
l2_penalty = regularization * tf.reduce_sum(beta ** 2)
total_loss = log_likelihood + l2_penalty
return total_loss
start = np.random.randn(dim + 1)
@tf.function
def l2_regression_with_lbfgs():
return tfp.optimizer.lbfgs_minimize(
negative_log_likelihood,
initial_position=tf.constant(start),
tolerance=1e-8)
results = run(l2_regression_with_lbfgs)
minimum = results.position
fitted_intercept = minimum[0]
fitted_beta = minimum[1:]
print('Converged:', results.converged)
print('Intercept: Fitted ({}), Actual ({})'.format(fitted_intercept, intercept))
print('Beta:\n\tFitted {},\n\tActual {}'.format(fitted_beta, betas))
Evaluation took: 0.056751 seconds Converged: True Intercept: Fitted (1.4111415084244365), Actual (1.3934058329729904) Beta: Fitted [-0.18016612 0.53121578 -0.56420632 -0.5336374 2.00499675], Actual [-0.20470766 0.47894334 -0.51943872 -0.5557303 1.96578057]
バッチ処理サポート
BFGS と L-BFGS の両者はバッチ化計算をサポートします、例えば多くの異なる開始点から単一関数を ; あるいは単一ポイントから複数のパラメトリック関数を最適化します。
単一関数、マルチ開始点
Himmelblau 関数は標準的な最適化テストケースです。関数は次により与えられます :
$$f(x, y) = (x^2 + y – 11)^2 + (x + y^2 – 7)^2$$
関数は以下に位置する 4 つの最小値を持ちます :
- (3, 2),
- (-2.805118, 3.131312),
- (-3.779310, -3.283186),
- (3.584428, -1.848126).
これら総ての最小値は適切な開始点から到達されるかもしれません。
# The function to minimize must take as input a tensor of shape [..., n]. In
# this n=2 is the size of the domain of the input and [...] are batching
# dimensions. The return value must be of shape [...], i.e. a batch of scalars
# with the objective value of the function evaluated at each input point.
@make_val_and_grad_fn
def himmelblau(coord):
x, y = coord[..., 0], coord[..., 1]
return (x * x + y - 11) ** 2 + (x + y * y - 7) ** 2
starts = tf.constant([[1, 1],
[-2, 2],
[-1, -1],
[1, -2]], dtype='float64')
# The stopping_condition allows to further specify when should the search stop.
# The default, tfp.optimizer.converged_all, will proceed until all points have
# either converged or failed. There is also a tfp.optimizer.converged_any to
# stop as soon as the first point converges, or all have failed.
@tf.function
def batch_multiple_starts():
return tfp.optimizer.lbfgs_minimize(
himmelblau, initial_position=starts,
stopping_condition=tfp.optimizer.converged_all,
tolerance=1e-8)
results = run(batch_multiple_starts)
print('Converged:', results.converged)
print('Minima:', results.position)
Evaluation took: 0.019095 seconds Converged: [ True True True True] Minima: [[ 3. 2. ] [-2.80511809 3.13131252] [-3.77931025 -3.28318599] [ 3.58442834 -1.84812653]]
複数の関数
デモ目的のため、この例では多くの高次元のランダムに生成された二次 bowls を同時に最適化します。
np.random.seed(12345)
dim = 100
batches = 500
minimum = np.random.randn(batches, dim)
scales = np.exp(np.random.randn(batches, dim))
@make_val_and_grad_fn
def quadratic(x):
return tf.reduce_sum(input_tensor=scales * (x - minimum)**2, axis=-1)
# Make all starting points (1, 1, ..., 1). Note not all starting points need
# to be the same.
start = tf.ones((batches, dim), dtype='float64')
@tf.function
def batch_multiple_functions():
return tfp.optimizer.lbfgs_minimize(
quadratic, initial_position=start,
stopping_condition=tfp.optimizer.converged_all,
max_iterations=100,
tolerance=1e-8)
results = run(batch_multiple_functions)
print('All converged:', np.all(results.converged))
print('Largest error:', np.max(results.position - minimum))
Evaluation took: 1.994132 seconds All converged: True Largest error: 4.4131473142527966e-08
以上
TensorFlow Probability 0.10 : ガイド : 学習可能な分布 Zoo
TensorFlow Probability 0.10 : ガイド : 学習可能な分布 Zoo (翻訳/解説)
翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 06/18/2020 (0.10.0)
* 本ページは、TensorFlow Probability の以下のドキュメントを翻訳した上で適宜、補足説明したものです：
* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。
ガイド : 学習可能な分布 Zoo
この colab では学習可能な (「訓練可能な」) 分布を構築する様々な例を示します (分布を説明する努力はしません、それらをどのように構築するかを示すだけです)。
import numpy as np import tensorflow.compat.v2 as tf import tensorflow_probability as tfp from tensorflow_probability.python.internal import prefer_static tfb = tfp.bijectors tfd = tfp.distributions tf.enable_v2_behavior()
event_size = 4 num_components = 3
chol(Cov) のための Scaled Identity を持つ学習可能な多変量正規分布
learnable_mvn_scaled_identity = tfd.Independent(
tfd.Normal(
loc=tf.Variable(tf.zeros(event_size), name='loc'),
scale=tfp.util.TransformedVariable(
tf.ones([event_size, 1]),
bijector=tfb.Exp()),
name='scale'),
reinterpreted_batch_ndims=1,
name='learnable_mvn_scaled_identity')
print(learnable_mvn_scaled_identity)
print(learnable_mvn_scaled_identity.trainable_variables)
tfp.distributions.Independent("learnable_mvn_scaled_identity", batch_shape=[4], event_shape=[4], dtype=float32)
(<tf.Variable 'Variable:0' shape=(4,) dtype=float32, numpy=array([0., 0., 0., 0.], dtype=float32)>, <tf.Variable 'Variable:0' shape=(4, 1) dtype=float32, numpy=
array([[0.],
[0.],
[0.],
[0.]], dtype=float32)>)
chol(Cov) のための Diagonal を持つ学習可能な多変量正規分布
learnable_mvndiag = tfd.Independent(
tfd.Normal(
loc=tf.Variable(tf.zeros(event_size), name='loc'),
scale=tfp.util.TransformedVariable(
tf.ones(event_size),
bijector=tfb.Softplus()), # Use Softplus...cuz why not?
name='scale'),
reinterpreted_batch_ndims=1,
name='learnable_mvn_diag')
print(learnable_mvndiag)
print(learnable_mvndiag.trainable_variables)
tfp.distributions.Independent("learnable_mvn_diag", batch_shape=[], event_shape=[4], dtype=float32)
(<tf.Variable 'Variable:0' shape=(4,) dtype=float32, numpy=array([0., 0., 0., 0.], dtype=float32)>, <tf.Variable 'Variable:0' shape=(4,) dtype=float32, numpy=array([0., 0., 0., 0.], dtype=float32)>)
多変量正規分布の混合 (球面 (= spherical))
learnable_mix_mvn_scaled_identity = tfd.MixtureSameFamily(
mixture_distribution=tfd.Categorical(
logits=tf.Variable(
# Changing the `1.` intializes with a geometric decay.
-tf.math.log(1.) * tf.range(num_components, dtype=tf.float32),
name='logits')),
components_distribution=tfd.Independent(
tfd.Normal(
loc=tf.Variable(
tf.random.normal([num_components, event_size]),
name='loc'),
scale=tfp.util.TransformedVariable(
10. * tf.ones([num_components, 1]),
bijector=tfb.Softplus()), # Use Softplus...cuz why not?
name='scale'),
reinterpreted_batch_ndims=1),
name='learnable_mix_mvn_scaled_identity')
print(learnable_mix_mvn_scaled_identity)
print(learnable_mix_mvn_scaled_identity.trainable_variables)
tfp.distributions.MixtureSameFamily("learnable_mix_mvn_scaled_identity", batch_shape=[], event_shape=[4], dtype=float32)
(<tf.Variable 'Variable:0' shape=(3, 4) dtype=float32, numpy=
array([[-0.11800266, -1.3127382 , -0.1813932 , 0.24975719],
[ 0.97209346, -0.91694 , 0.84734786, -1.3201759 ],
[ 0.15194772, 0.33787215, -0.4269769 , -0.26286215]],
dtype=float32)>, <tf.Variable 'Variable:0' shape=(3, 1) dtype=float32, numpy=
array([[-4.600166],
[-4.600166],
[-4.600166]], dtype=float32)>, <tf.Variable 'logits:0' shape=(3,) dtype=float32, numpy=array([-0., -0., -0.], dtype=float32)>)
多変量正規分布の混合 (球面) with first mix weight unlearnable
learnable_mix_mvndiag_first_fixed = tfd.MixtureSameFamily(
mixture_distribution=tfd.Categorical(
logits=tfp.util.TransformedVariable(
# Initialize logits as geometric decay.
-tf.math.log(1.5) * tf.range(num_components, dtype=tf.float32),
tfb.Pad(paddings=[[1, 0]], constant_values=0)),
name='logits'),
components_distribution=tfd.Independent(
tfd.Normal(
loc=tf.Variable(
# Use Rademacher...cuz why not?
tfp.math.random_rademacher([num_components, event_size]),
name='loc'),
scale=tfp.util.TransformedVariable(
10. * tf.ones([num_components, 1]),
bijector=tfb.Softplus()), # Use Softplus...cuz why not?
name='scale'),
reinterpreted_batch_ndims=1),
name='learnable_mix_mvndiag_first_fixed')
print(learnable_mix_mvndiag_first_fixed)
print(learnable_mix_mvndiag_first_fixed.trainable_variables)
tfp.distributions.MixtureSameFamily("learnable_mix_mvndiag_first_fixed", batch_shape=[], event_shape=[4], dtype=float32)
(<tf.Variable 'Variable:0' shape=(3, 4) dtype=float32, numpy=
array([[-1., -1., 1., -1.],
[-1., -1., -1., -1.],
[ 1., 1., 1., -1.]], dtype=float32)>, <tf.Variable 'Variable:0' shape=(3, 4) dtype=float32, numpy=
array([[-4.600166, -4.600166, -4.600166, -4.600166],
[-4.600166, -4.600166, -4.600166, -4.600166],
[-4.600166, -4.600166, -4.600166, -4.600166]], dtype=float32)>, <tf.Variable 'Variable:0' shape=(2,) dtype=float32, numpy=array([-0.4054651, -0.8109302], dtype=float32)>)
多変量正規分布の混合 (full Cov)
earnable_mix_mvntril = tfd.MixtureSameFamily(
mixture_distribution=tfd.Categorical(
logits=tf.Variable(
# Changing the `1.` intializes with a geometric decay.
-tf.math.log(1.) * tf.range(num_components, dtype=tf.float32),
name='logits')),
components_distribution=tfd.MultivariateNormalTriL(
loc=tf.Variable(tf.zeros([num_components, event_size]), name='loc'),
scale_tril=tfp.util.TransformedVariable(
10. * tf.eye(event_size, batch_shape=[num_components]),
bijector=tfb.FillScaleTriL()),
name='scale_tril'),
name='learnable_mix_mvntril')
print(learnable_mix_mvntril)
print(learnable_mix_mvntril.trainable_variables)
tfp.distributions.MixtureSameFamily("learnable_mix_mvntril", batch_shape=[], event_shape=[4], dtype=float32)
(<tf.Variable 'Variable:0' shape=(3, 4) dtype=float32, numpy=
array([[ 1.4036556 , -0.22486973, 0.8365339 , 1.1744921 ],
[-0.14385273, 1.5095806 , 0.78327304, -0.64133334],
[-0.22640549, 1.908316 , 1.1216396 , -1.2109828 ]],
dtype=float32)>, <tf.Variable 'Variable:0' shape=(3, 10) dtype=float32, numpy=
array([[-4.6011715, 0. , 0. , 0. , -4.6011715,
-4.6011715, 0. , 0. , 0. , -4.6011715],
[-4.6011715, 0. , 0. , 0. , -4.6011715,
-4.6011715, 0. , 0. , 0. , -4.6011715],
[-4.6011715, 0. , 0. , 0. , -4.6011715,
-4.6011715, 0. , 0. , 0. , -4.6011715]],
dtype=float32)>, <tf.Variable 'logits:0' shape=(3,) dtype=float32, numpy=array([-0., -0., -0.], dtype=float32)>)
多変量正規分布の混合 (full Cov) with unlearnable first mix & first component
# Make a bijector which pads an eye to what otherwise fills a tril.
num_tril_nonzero = lambda num_rows: num_rows * (num_rows + 1) // 2
num_tril_rows = lambda nnz: prefer_static.cast(
prefer_static.sqrt(0.25 + 2. * prefer_static.cast(nnz, tf.float32)) - 0.5,
tf.int32)
# TFP doesn't have a concat bijector, so we roll out our own.
class PadEye(tfb.Bijector):
def __init__(self, tril_fn=None):
if tril_fn is None:
tril_fn = tfb.FillScaleTriL()
self._tril_fn = getattr(tril_fn, 'inverse', tril_fn)
super(PadEye, self).__init__(
forward_min_event_ndims=2,
inverse_min_event_ndims=2,
is_constant_jacobian=True,
name='PadEye')
def _forward(self, x):
num_rows = int(num_tril_rows(tf.compat.dimension_value(x.shape[-1])))
eye = tf.eye(num_rows, batch_shape=prefer_static.shape(x)[:-2])
return tf.concat([self._tril_fn(eye)[..., tf.newaxis, :], x],
axis=prefer_static.rank(x) - 2)
def _inverse(self, y):
return y[..., 1:, :]
def _forward_log_det_jacobian(self, x):
return tf.zeros([], dtype=x.dtype)
def _inverse_log_det_jacobian(self, y):
return tf.zeros([], dtype=y.dtype)
def _forward_event_shape(self, in_shape):
n = prefer_static.size(in_shape)
return in_shape + prefer_static.one_hot(n - 2, depth=n, dtype=tf.int32)
def _inverse_event_shape(self, out_shape):
n = prefer_static.size(out_shape)
return out_shape - prefer_static.one_hot(n - 2, depth=n, dtype=tf.int32)
tril_bijector = tfb.FillScaleTriL(diag_bijector=tfb.Softplus())
learnable_mix_mvntril_fixed_first = tfd.MixtureSameFamily(
mixture_distribution=tfd.Categorical(
logits=tfp.util.TransformedVariable(
# Changing the `1.` intializes with a geometric decay.
-tf.math.log(1.) * tf.range(num_components, dtype=tf.float32),
bijector=tfb.Pad(paddings=[(1, 0)]),
name='logits')),
components_distribution=tfd.MultivariateNormalTriL(
loc=tfp.util.TransformedVariable(
tf.zeros([num_components, event_size]),
bijector=tfb.Pad(paddings=[(1, 0)], axis=-2),
name='loc'),
scale_tril=tfp.util.TransformedVariable(
10. * tf.eye(event_size, batch_shape=[num_components]),
bijector=tfb.Chain([tril_bijector, PadEye(tril_bijector)]),
name='scale_tril')),
name='learnable_mix_mvntril_fixed_first')
print(learnable_mix_mvntril_fixed_first)
print(learnable_mix_mvntril_fixed_first.trainable_variables)
tfp.distributions.MixtureSameFamily("learnable_mix_mvntril_fixed_first", batch_shape=[], event_shape=[4], dtype=float32)
(<tf.Variable 'loc:0' shape=(2, 4) dtype=float32, numpy=
array([[ 0.53900903, -0.17989647, -1.196744 , -1.0601326 ],
[ 0.46199334, 1.2968503 , 0.20908853, -0.36455044]],
dtype=float32)>, <tf.Variable 'Variable:0' shape=(2, 10) dtype=float32, numpy=
array([[-4.6011715, 0. , 0. , 0. , -4.6011715,
-4.6011715, 0. , 0. , 0. , -4.6011715],
[-4.6011715, 0. , 0. , 0. , -4.6011715,
-4.6011715, 0. , 0. , 0. , -4.6011715]],
dtype=float32)>, <tf.Variable 'logits:0' shape=(2,) dtype=float32, numpy=array([-0., -0.], dtype=float32)>)
以上
TensorFlow Probability 0.10 : ガイド : TensorFlow Distributions : 優しいイントロダクション
TensorFlow Probability 0.10 : ガイド : TensorFlow Distributions : 優しいイントロダクション (翻訳/解説)
翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 06/18/2020 (0.10.0)
* 本ページは、TensorFlow Probability の以下のドキュメントを翻訳した上で適宜、補足説明したものです：
* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。
ガイド : TensorFlow Distributions : 優しいイントロダクション
このノートブックでは、TensorFlow 分布 (TFD for short) を探求します。このノートブックの目的は tensor shape の TFD の処理を理解することを含めて、貴方に学習カーブを優しく上げて頂くことです。このノートブックは抽象的な概念の前に例を提示しようとします。最初に物事を進めるための標準的で容易な方法を提示して、最も一般的な抽象的な眺めは最後まで取っておきます。もし貴方がより抽象的でリファレンス・スタイルのチュートリアルを好みタイプであれば、TensorFlow 分布 shape の理解 を調べてください。ここでのマテリアルについてどのような質問でもあれば、TensorFlow Probability メーリングリスト にコンタクト (or 参加) することを躊躇しないでください。We’re happy to help.
始める前に、適切なライブラリをインポートする必要があります。私達の全体的なライブラリは tensorflow_probability です。慣習により、一般的に distributions ライブラリを tfd として参照します。
Tensorflow Eager は TensorFlow のための命令型実行環境です。TensorFlow eager では、総ての TF 演算は直ちに評価されて結果を生成します。これは TensorFlow の標準的な「グラフ」モードとは対称的です、そこでは TF 演算はグラフに (後で実行される) ノードを追加します。ここで表わされる概念のどれもそれに依存はしてませんが、このノートブック全体は TF Eager を使用して書かれます、そして TFP はグラフモードでも利用できます。
import collections import tensorflow as tf import tensorflow_probability as tfp tfd = tfp.distributions try: tf.compat.v1.enable_eager_execution() except ValueError: pass import matplotlib.pyplot as plt
基本的な単変量分布
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: shape=(), dtype=float32, numpy=0.25322816>
複数のサンプルをドローできます :
n.sample(3)
<tf.Tensor: shape=(3,), dtype=float32, numpy=array([-1.4658079, -0.5653636, 0.9314412], dtype=float32)>
log prob を評価できます :
n.log_prob(0.)
<tf.Tensor: shape=(), dtype=float32, numpy=-0.9189385>
複数の対数確率を評価できます :
n.log_prob([0., 2., 4.])
<tf.Tensor: 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: shape=(), dtype=int32, numpy=1>
b.sample(8)
<tf.Tensor: shape=(8,), dtype=int32, numpy=array([1, 0, 0, 0, 1, 0, 1, 0], dtype=int32)>
b.log_prob(1)
<tf.Tensor: shape=(), dtype=float32, numpy=-0.35667497>
b.log_prob([1, 0, 1, 0])
<tf.Tensor: shape=(4,), dtype=float32, numpy=array([-0.35667497, -1.2039728 , -0.35667497, -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: shape=(2,), dtype=float32, numpy=array([-1.2489667, 15.025171 ], dtype=float32)>
nd.sample(5)
<tf.Tensor: shape=(5, 2), dtype=float32, numpy=
array([[-1.5439653 , 8.9968405 ],
[-0.38730723, 12.448896 ],
[-0.8697963 , 9.330035 ],
[-1.2541095 , 10.268944 ],
[ 2.3475595 , 13.184147 ]], dtype=float32)>
nd.log_prob([0., 10])
<tf.Tensor: 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: shape=(3,), dtype=int32, numpy=array([0, 1, 1], dtype=int32)>
b3.sample(6)
<tf.Tensor: shape=(6, 3), dtype=int32, numpy=
array([[1, 0, 1],
[0, 1, 1],
[0, 0, 1],
[0, 0, 1],
[0, 0, 1],
[0, 1, 0]], dtype=int32)>
prob を呼び出す場合、(これは log_prob と同じ shape セマンティクスを持ちます ; 明確にするためにこれらの小さな Bernoulli サンプルで prob を使用します、アプリケーションでは通常は log_prob が好まれますが) それにベクトルを渡してその値を生成する各コインの確率を評価できます :
b3.prob([1, 1, 0])
<tf.Tensor: shape=(3,), dtype=float32, numpy=array([0.29999998, 0.5 , 0.29999998], dtype=float32)>
何故 API はバッチ shape を含むのでしょう？意味的には、分布のリストを作成して for ループでそれらに渡り反復することで同じ計算を遂行できるでしょう (少なくとも Eager モードで、TF graph モードでは tf.while ループが必要でしょう)。けれども、同一にパラメータ化された分布の (潜在的に巨大な) セットを持つことは極めて一般的で、そして可能なときはいつでもベクトル化された計算の使用はハードウェア・アクセラレータを使用して高速な計算を遂行することを可能にする点において主要な構成要素です。
Independent を使用してバッチを事象に集積する (= aggregate)
前のセクションでは、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 はバッチ shape を事象 shape に移したことを見ます : b3_joint は 3-次元事象空間 (event_shape = (3,)) に渡る単一の分布 (batch_shape = ()) です。
セマンティクスをチェックしましょう :
b3_joint.prob([1, 1, 0])
<tf.Tensor: shape=(), dtype=float32, numpy=0.044999998>
同じ結果を得るための代替の方法は b3 を使用して確率を計算して乗算により手動で reduction を行なうことです (あるいは、対数確率が使用されるより通常のケースでは、合計します) :
tf.reduce_prod(b3.prob([1, 1, 0]))
<tf.Tensor: 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: shape=(4, 3, 2), dtype=float32, numpy=
array([[[ 0.7367498 , 2.730996 ],
[-0.74080074, -0.36466932],
[ 0.6516018 , 0.9391426 ]],
[[ 1.038303 , 0.12231752],
[-0.94788766, -1.204232 ],
[ 4.059758 , 3.035752 ]],
[[ 0.56903946, -0.06875849],
[-0.35127294, 0.5311631 ],
[ 3.4635801 , 4.565582 ]],
[[-0.15989424, -0.25715637],
[ 0.87479895, 0.97391707],
[ 0.5211419 , 2.32108 ]]], dtype=float32)>
batch_shape = (3,) そして event_shape = (2,) ですから、log_prob に shape (3, 2) の tensor を渡します :
nd_batch.log_prob([[0., 0.], [1., 1.], [2., 2.]])
<tf.Tensor: shape=(3,), dtype=float32, numpy=array([-1.8328519, -1.7907217, -1.694036 ], dtype=float32)>
ブロードキャスティング, aka 何故これはそれほど混乱させるのでしょう？
ここまで行なったことから要約すると、総ての分布はバッチ shape B と事象 shape E を持ちます。BE を 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: shape=(), dtype=float32, numpy=-0.9189385>
n.log_prob([0.])
<tf.Tensor: shape=(1,), dtype=float32, numpy=array([-0.9189385], dtype=float32)>
n.log_prob([[0., 1.], [-1., 2.]])
<tf.Tensor: 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: shape=(), dtype=float32, numpy=-2.337877>
しかし「より多くの」サンプルを渡して、そしてそれらの log_prob 総てを一度に評価することもできます :
nd.log_prob([[0., 0.],
[1., 1.],
[2., 2.]])
<tf.Tensor: shape=(3,), dtype=float32, numpy=array([-2.337877 , -2.337877 , -4.3378773], dtype=float32)>
多分魅力は少ないですが、事象次元に渡りブロードキャストすることもできます :
nd.log_prob([0.])
<tf.Tensor: shape=(), dtype=float32, numpy=-2.337877>
nd.log_prob([[0.], [1.], [2.]])
<tf.Tensor: 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: 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: shape=(2, 3), dtype=float32, numpy=
array([[0.7, 0.5, 0.3],
[0.3, 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 操作を更に深く掘り下げて調べるためには、TensorFlow 分布 shape の理解 を見てください。
- TFP は Bijector として知られるパワフルな抽象を含みます、これは TransformedDistribution と連携して (既存の分布の可逆な (= invertible) 変換である) 新しい分布を容易に作成するための柔軟で、合成的な方法を与えます。私達はこれについてのチュートリアルを間もなく書こうとしますが、当面は、この ドキュメント をチェックしてください。
以上
TensorFlow Probability 0.10 : ガイド : TensorFlow 分布 shape の理解
TensorFlow Probability 0.10 : ガイド : TensorFlow 分布 shape の理解 (翻訳/解説)
翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 06/18/2020 (0.10.0)
* 本ページは、TensorFlow Probability の以下のドキュメントを翻訳した上で適宜、補足説明したものです：
* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。
ガイド : TensorFlow 分布 shape を理解する
import collections import tensorflow as tf tf.compat.v2.enable_v2_behavior() import tensorflow_probability as tfp tfd = tfp.distributions tfb = tfp.bijectors
基本
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 オブジェクトが作成されたときに評価されます (そして従って知られます)、その一方でグラフ (非 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 を計算する) 位置を表わす (非-empty) tensor を取ります。最も簡単なケースでは、この 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: shape=(2, 3), dtype=float32, numpy=
array([[ -1. , -2.0785608, -3.2223587],
[-364.73938 , -2.0785608, -95.39484 ]], dtype=float32)>
three_poissons.log_prob([[[[1., 10., 100.], [100., 10., 1.]]]]) # sample_shape is [1, 1, 2].
<tf.Tensor: shape=(1, 1, 2, 3), dtype=float32, numpy=
array([[[[ -1. , -2.0785608, -3.2223587],
[-364.73938 , -2.0785608, -95.39484 ]]]], dtype=float32)>
最初の例でどのように入力と出力が shape [2, 3] を持つか、そして 2 番目の例でそれらが shape [1, 1, 2, 3] を持つか注意してください。
言うべきことはそれが総てでしょう、もしそれがブロードキャスティングのためでないのであれば。ここにブロードキャスティングをひとたび考慮した場合のルールがあります。それを完全に一般化して記述してそしてスカラー分布のための単純化を記します :
- n = len(batch_shape) + len(event_shape) と定義します。(スカラー分布については、len(event_shape)=0。)
- 入力 tensor t が n 次元よりも少ない場合、それが正確に n 次元を持つまでその shape を size 1 の次元を左に追加してパッドします。結果の tensor を t’ と呼びます。
- t’ の n 右端次元を (そのために log_prob を計算している) 分布の [batch_shape, event_shape] に対してブロードキャストします。より詳細には : t’ が既に分布に一致している次元については、何もしません、そして t’ がシングルトンを持つような次元については、そのシングルトンを適切な回数複製 (= replicate) します。他の状況は総てエラーです。(スカラー分布については、batch_shape に対してブロードキャストするだけです、何故ならば event_shape = [] だからです。)
- 今は最後に log_prob を計算することができます。結果の tensor は shape [sample_shape, batch_shape] を持ち、そこでは sample_shape は t か t’ の任意の次元として n-右端次元の左側に定義されます : sample_shape = shape(t)[:-n]。
それが何を意味するか知らないのであればこれは混乱しますので、幾つかの例を動作させましょう :
three_poissons.log_prob([10.])
<tf.Tensor: shape=(3,), dtype=float32, numpy=array([-16.104412 , -2.0785608, -69.05272 ], dtype=float32)>
tensor [10.] (with shape [1]) は 3 の batch_shape に渡りブロードキャストされますので、3 つ総ての Poisson の対数確率を値 10 で評価します。
three_poissons.log_prob([[[1.], [10.]], [[100.], [1000.]]])
<tf.Tensor: shape=(2, 2, 3), dtype=float32, numpy=
array([[[-1.0000000e+00, -7.6974149e+00, -9.5394836e+01],
[-1.6104412e+01, -2.0785608e+00, -6.9052719e+01]],
[[-3.6473938e+02, -1.4348087e+02, -3.2223587e+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: shape=(2, 3), dtype=float32, numpy=
array([[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -17.004269 , -194.70169 ]], dtype=float32)>
poisson_2_by_3.log_prob([1.]) # Exactly equivalent to above, demonstrating the scalar special case.
<tf.Tensor: shape=(2, 3), dtype=float32, numpy=
array([[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -17.004269 , -194.70169 ]], 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: shape=(2, 3), dtype=float32, numpy=
array([[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -17.004269 , -194.70169 ]], dtype=float32)>
poisson_2_by_3.log_prob([[1., 10., 100.]]) # Input is [1, 3] broadcast to [2, 3].
<tf.Tensor: shape=(2, 3), dtype=float32, numpy=
array([[ -1. , -2.0785608, -3.2223587],
[ -1.3068528, -5.14709 , -33.90767 ]], dtype=float32)>
poisson_2_by_3.log_prob([[1., 10., 100.], [1., 10., 100.]]) # Equivalent to above. No broadcasting.
<tf.Tensor: shape=(2, 3), dtype=float32, numpy=
array([[ -1. , -2.0785608, -3.2223587],
[ -1.3068528, -5.14709 , -33.90767 ]], dtype=float32)>
poisson_2_by_3.log_prob([[1., 1., 1.], [2., 2., 2.]]) # No broadcasting.
<tf.Tensor: shape=(2, 3), dtype=float32, numpy=
array([[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -14.701683 , -190.09653 ]], dtype=float32)>
poisson_2_by_3.log_prob([[1.], [2.]]) # Equivalent to above. Input shape [2, 1] broadcast to [2, 3].
<tf.Tensor: shape=(2, 3), dtype=float32, numpy=
array([[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -14.701683 , -190.09653 ]], 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: shape=(2, 2, 3), dtype=float32, numpy=
array([[[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -17.004269 , -194.70169 ]],
[[ -1.6931472, -6.087977 , -91.48282 ],
[ -1.3068528, -14.701683 , -190.09653 ]]], dtype=float32)>
あるいはブロードキャスティングに最後のバッチ次元を処理させることもできます :
poisson_2_by_3.log_prob([[[1.], [1.]], [[2.], [2.]]]) # Input shape [2, 2, 1].
<tf.Tensor: shape=(2, 2, 3), dtype=float32, numpy=
array([[[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -17.004269 , -194.70169 ]],
[[ -1.6931472, -6.087977 , -91.48282 ],
[ -1.3068528, -14.701683 , -190.09653 ]]], dtype=float32)>
ブロードキャスティングにちょうど最初のバッチ次元を処理させることもできます (多分幾分自然ではないですが) :
poisson_2_by_3.log_prob([[[1., 1., 1.]], [[2., 2., 2.]]]) # Input shape [2, 1, 3].
<tf.Tensor: shape=(2, 2, 3), dtype=float32, numpy=
array([[[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -17.004269 , -194.70169 ]],
[[ -1.6931472, -6.087977 , -91.48282 ],
[ -1.3068528, -14.701683 , -190.09653 ]]], dtype=float32)>
あるいはブロードキャスティングに両者のバッチ次元を処理させられるでしょう :
poisson_2_by_3.log_prob([[[1.]], [[2.]]]) # Input shape [2, 1, 1].
<tf.Tensor: shape=(2, 2, 3), dtype=float32, numpy=
array([[[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -17.004269 , -194.70169 ]],
[[ -1.6931472, -6.087977 , -91.48282 ],
[ -1.3068528, -14.701683 , -190.09653 ]]], dtype=float32)>
望んだ 2 つの値だけを持ったとき上は素晴らしく動作しましたが、総てのバッチポイントで評価することを望む値の長いリストを持つと仮定します。そのためには、shape の右側に size 1 の特別な次元を追加する、次の記法が非常に有用です :
poisson_2_by_3.log_prob(tf.constant([1., 2.])[..., tf.newaxis, tf.newaxis])
<tf.Tensor: shape=(2, 2, 3), dtype=float32, numpy=
array([[[ -1. , -7.697415 , -95.39484 ],
[ -1.3068528, -17.004269 , -194.70169 ]],
[[ -1.6931472, -6.087977 , -91.48282 ],
[ -1.3068528, -14.701683 , -190.09653 ]]], dtype=float32)>
これは (知るに値する) strided slice 記法 (訳注: リンク切れ) のインスタンスです。
完全性のために three_poissons に戻ります、同じ例は次のように見えます :
three_poissons.log_prob([[1.], [10.], [50.], [100.]])
<tf.Tensor: shape=(4, 3), dtype=float32, numpy=
array([[ -1. , -7.697415 , -95.39484 ],
[ -16.104412 , -2.0785608, -69.05272 ],
[-149.47777 , -43.34851 , -18.219261 ],
[-364.73938 , -143.48087 , -3.2223587]], dtype=float32)>
three_poissons.log_prob(tf.constant([1., 10., 50., 100.])[..., tf.newaxis]) # Equivalent to above.
<tf.Tensor: shape=(4, 3), dtype=float32, numpy=
array([[ -1. , -7.697415 , -95.39484 ],
[ -16.104412 , -2.0785608, -69.05272 ],
[-149.47777 , -43.34851 , -18.219261 ],
[-364.73938 , -143.48087 , -3.2223587]], dtype=float32)>
多変量分布
さて多変量分布に向かいます、これは非-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 を持つように内部でブロードキャストしています。
既に私達が知っていることを考えれば、サンプリングは簡単です :
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)
対数確率の計算も等しく簡単です。対角多変量正規分布の例を動かしてみましょう。(多変量はそれほど都合よくはブロードキャストされません、何故ならば 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: 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: 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: 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: shape=(), dtype=float32, numpy=-178.22021>
transformed_multinomial.log_prob([[500., 100., 100.], [150., 100., 50.]])
<tf.Tensor: shape=(), dtype=float32, numpy=-178.22021>
これが Reshape bijector のできる唯一のことです : それは事象次元をバッチ次元に変更することはできません or vice-versa。
独立分布
独立 (= Independent) 分布は独立で、必ずしも同一でない分布 (aka 分布のバッチ) を単一の分布として扱うために使用されます。より簡潔には、Independent は 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: shape=(2, 5), dtype=float32, numpy=
array([[-2.9957323 , -0.10536052, -0.16251892, -1.609438 , -0.2876821 ],
[-0.35667497, -0.4307829 , -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: 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: 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 0.10 : ガイド : 概要 – TensorFlow Probability のツアー
TensorFlow Probability 0.10 : ガイド : 概要 – TensorFlow Probability のツアー (翻訳/解説)
翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 06/17/2020 (0.10.0)
* 本ページは、TensorFlow Probability の以下のドキュメントを翻訳した上で適宜、補足説明したものです：
* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。
ガイド : 概要 – TensorFlow Probability のツアー
依存性 & 要件
!pip install -q --upgrade seaborn
#@title Import { display-mode: "form" }
from pprint import pprint
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
import tensorflow_probability as tfp
sns.reset_defaults()
sns.set_context(context='talk',font_scale=0.7)
plt.rcParams['image.cmap'] = 'viridis'
%matplotlib inline
tfd = tfp.distributions
tfb = tfp.bijectors
#@title Utils { display-mode: "form" }
def print_subclasses_from_module(module, base_class, maxwidth=80):
import functools, inspect, sys
subclasses = [name for name, obj in inspect.getmembers(module)
if inspect.isclass(obj) and issubclass(obj, base_class)]
def red(acc, x):
if not acc or len(acc[-1]) + len(x) + 2 > maxwidth:
acc.append(x)
else:
acc[-1] += ", " + x
return acc
print('\n'.join(functools.reduce(red, subclasses, [])))
アウトライン
- TensorFlow
- TensorFlow Probability
- Distributions
- Bijectors
- MCMC
- …and more!
序: TensorFlow
TensorFlow は科学計算ライブラリです。
それは以下をサポートします :
- 多くの数学的演算
- 効率的なベクトル化計算
- 容易なハードウェア・アクセラレーション
- 自動微分
ベクトル化
- ベクトル化は物事を高速にします！
- それはまた shape について多く考えることも意味します。
mats = tf.random.uniform(shape=[1000, 10, 10])
vecs = tf.random.uniform(shape=[1000, 10, 1])
def for_loop_solve():
return np.array(
[tf.linalg.solve(mats[i, ...], vecs[i, ...]) for i in range(1000)])
def vectorized_solve():
return tf.linalg.solve(mats, vecs)
# Vectorization for the win!
%timeit for_loop_solve()
%timeit vectorized_solve()
1 loops, best of 3: 2.49 s per loop 1000 loops, best of 3: 836 µs per loop
ハードウェア・アクセラレーション
# Code can run seamlessly on a GPU, just change Colab runtime type
# in the 'Runtime' menu.
if tf.test.gpu_device_name() == '/device:GPU:0':
print("Using a GPU")
else:
print("Using a CPU")
Using a CPU
自動微分
a = tf.constant(np.pi) b = tf.constant(np.e) with tf.GradientTape() as tape: tape.watch([a, b]) c = .5 * (a**2 + b**2) grads = tape.gradient(c, [a, b]) print(grads[0]) print(grads[1])
tf.Tensor(3.1415927, shape=(), dtype=float32) tf.Tensor(2.7182817, shape=(), dtype=float32)
TensorFlow Probability
TensorFlow Probability は TensorFlow で確率的推論と統計的分析のためのライブラリです。
低位モジュールコンポーネントの構成を通してモデル化、推論と評価をサポートします。
低位ビルディングブロック
- Distributions
- Bijectors
高位構成概念
- マルコフ連鎖モンテカルロ法
- 確率的層
- 構造時系列
- 一般化線形モデル
- Edward2
- Optimizers
Distributions
tfp.distributions.Distribution は 2 つの中心的なメソッド: sample と log_prob を持つクラスです。
TFP has a lot of distributions!
print_subclasses_from_module(tfp.distributions, tfp.distributions.Distribution)
Autoregressive, BatchReshape, Bernoulli, Beta, BetaBinomial, Binomial, Blockwise Categorical, Cauchy, Chi, Chi2, CholeskyLKJ, ContinuousBernoulli, Deterministic Dirichlet, DirichletMultinomial, Distribution, DoublesidedMaxwell, Empirical ExpRelaxedOneHotCategorical, Exponential, FiniteDiscrete, Gamma, GammaGamma GaussianProcess, GaussianProcessRegressionModel, GeneralizedPareto, Geometric Gumbel, HalfCauchy, HalfNormal, HalfStudentT, HiddenMarkovModel, Horseshoe Independent, InverseGamma, InverseGaussian, JointDistribution JointDistributionCoroutine, JointDistributionCoroutineAutoBatched JointDistributionNamed, JointDistributionNamedAutoBatched JointDistributionSequential, JointDistributionSequentialAutoBatched, Kumaraswamy LKJ, Laplace, LinearGaussianStateSpaceModel, LogNormal, Logistic, LogitNormal Mixture, MixtureSameFamily, Multinomial, MultivariateNormalDiag MultivariateNormalDiagPlusLowRank, MultivariateNormalFullCovariance MultivariateNormalLinearOperator, MultivariateNormalTriL MultivariateStudentTLinearOperator, NegativeBinomial, Normal, OneHotCategorical OrderedLogistic, PERT, Pareto, PixelCNN, PlackettLuce, Poisson PoissonLogNormalQuadratureCompound, ProbitBernoulli, QuantizedDistribution RelaxedBernoulli, RelaxedOneHotCategorical, Sample, SinhArcsinh, StudentT StudentTProcess, TransformedDistribution, Triangular, TruncatedNormal, Uniform VariationalGaussianProcess, VectorDeterministic, VectorDiffeomixture VectorExponentialDiag, VonMises, VonMisesFisher, WishartLinearOperator WishartTriL, Zipf
単純なスカラー変量分布
# A standard normal normal = tfd.Normal(loc=0., scale=1.) print(normal)
tfp.distributions.Normal("Normal", batch_shape=[], event_shape=[], dtype=float32)
# Plot 1000 samples from a standard normal
samples = normal.sample(1000)
sns.distplot(samples)
plt.title("Samples from a standard Normal")
plt.show()
# Compute the log_prob of a point in the event space of `normal` normal.log_prob(0.)
<tf.Tensor: shape=(), dtype=float32, numpy=-0.9189385>
# Compute the log_prob of a few points normal.log_prob([-1., 0., 1.])
<tf.Tensor: shape=(3,), dtype=float32, numpy=array([-1.4189385, -0.9189385, -1.4189385], dtype=float32)>
Distributions と Shapes
Numpy ndarray と TensorFlow Tensor は shape を持ちます。
TensorFlow Probability Distributions は shape セマンティクスを持ちます — shapes を意味的に別個のピースに分割します、全体総てのためにメモリの同じチャンク (Tensor/ndarray) が利用されてさえも。
- バッチ shape は制限パラメータを持つ Distributions のコレクション を意味します。
- 事象 (= Event) shape は Distribution からの サンプリングの shape を意味します。
常にバッチ shape は「左側」にそして事象 shape は「右側」に置きます。
スカラー変量分布のバッチ
バッチは「ベクトル化」分布のようなものです : 独立インスタンスでその計算は並列に発生します。
# Create a batch of 3 normals, and plot 1000 samples from each
normals = tfd.Normal([-2.5, 0., 2.5], 1.) # The scale parameter broadacasts!
print("Batch shape:", normals.batch_shape)
print("Event shape:", normals.event_shape)
Batch shape: (3,) Event shape: ()
# Samples' shapes go on the left!
samples = normals.sample(1000)
print("Shape of samples:", samples.shape)
Shape of samples: (1000, 3)
# Sample shapes can themselves be more complicated
print("Shape of samples:", normals.sample([10, 10, 10]).shape)
Shape of samples: (10, 10, 10, 3)
# A batch of normals gives a batch of log_probs. print(normals.log_prob([-2.5, 0., 2.5]))
tf.Tensor([-0.9189385 -0.9189385 -0.9189385], shape=(3,), dtype=float32)
# The computation broadcasts, so a batch of normals applied to a scalar # also gives a batch of log_probs. print(normals.log_prob(0.))
tf.Tensor([-4.0439386 -0.9189385 -4.0439386], shape=(3,), dtype=float32)
# Normal numpy-like broadcasting rules apply!
xs = np.linspace(-6, 6, 200)
try:
normals.log_prob(xs)
except Exception as e:
print("TFP error:", e.message)
TFP error: Incompatible shapes: [200] vs. [3] [Op:SquaredDifference]
# That fails for the same reason this does:
try:
np.zeros(200) + np.zeros(3)
except Exception as e:
print("Numpy error:", e)
Numpy error: operands could not be broadcast together with shapes (200,) (3,)
# But this would work:
a = np.zeros([200, 1]) + np.zeros(3)
print("Broadcast shape:", a.shape)
Broadcast shape: (200, 3)
# And so will this!
xs = np.linspace(-6, 6, 200)[..., np.newaxis]
# => shape = [200, 1]
lps = normals.log_prob(xs)
print("Broadcast log_prob shape:", lps.shape)
Broadcast log_prob shape: (200, 3)
# Summarizing visually
for i in range(3):
sns.distplot(samples[:, i], kde=False, norm_hist=True)
plt.plot(np.tile(xs, 3), normals.prob(xs), c='k', alpha=.5)
plt.title("Samples from 3 Normals, and their PDF's")
plt.show()
ベクトル変量分布
mvn = tfd.MultivariateNormalDiag(loc=[0., 0.], scale_diag = [1., 1.])
print("Batch shape:", mvn.batch_shape)
print("Event shape:", mvn.event_shape)
Batch shape: () Event shape: (2,)
samples = mvn.sample(1000)
print("Samples shape:", samples.shape)
Samples shape: (1000, 2)
行列変量分布
lkj = tfd.LKJ(dimension=10, concentration=[1.5, 3.0])
print("Batch shape: ", lkj.batch_shape)
print("Event shape: ", lkj.event_shape)
Batch shape: (2,) Event shape: (10, 10)
samples = lkj.sample()
print("Samples shape: ", samples.shape)
Samples shape: (2, 10, 10)
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(6, 3)) sns.heatmap(samples[0, ...], ax=axes[0], cbar=False) sns.heatmap(samples[1, ...], ax=axes[1], cbar=False) fig.tight_layout() plt.show()
ガウス過程
kernel = tfp.math.psd_kernels.ExponentiatedQuadratic()
xs = np.linspace(-5., 5., 200).reshape([-1, 1])
gp = tfd.GaussianProcess(kernel, index_points=xs)
print("Batch shape:", gp.batch_shape)
print("Event shape:", gp.event_shape)
Batch shape: () Event shape: (200,)
upper, lower = gp.mean() + [2 * gp.stddev(), -2 * gp.stddev()] plt.plot(xs, gp.mean()) plt.fill_between(xs[..., 0], upper, lower, color='k', alpha=.1) for _ in range(5): plt.plot(xs, gp.sample(), c='r', alpha=.3) plt.title(r"GP prior mean, $2\sigma$ intervals, and samples") plt.show() # *** Bonus question *** # Why do so many of these functions lie outside the 95% intervals?
GP 回帰
# Suppose we have some observed data obs_x = [[-3.], [0.], [2.]] # Shape 3x1 (3 1-D vectors) obs_y = [3., -2., 2.] # Shape 3 (3 scalars) gprm = tfd.GaussianProcessRegressionModel(kernel, xs, obs_x, obs_y)
upper, lower = gprm.mean() + [2 * gprm.stddev(), -2 * gprm.stddev()] plt.plot(xs, gprm.mean()) plt.fill_between(xs[..., 0], upper, lower, color='k', alpha=.1) for _ in range(5): plt.plot(xs, gprm.sample(), c='r', alpha=.3) plt.scatter(obs_x, obs_y, c='k', zorder=3) plt.title(r"GP posterior mean, $2\sigma$ intervals, and samples") plt.show()
Bijector
Bijector は (殆どの場合) 可逆 (= invertible)、滑らかな関数を表します。それらは分布を変換するために使用できて、サンプリングして log_probs を計算する能力を維持します。それらは tfp.bijectors モジュールに在り得ます。
各 bijector は少なくとも 3 つのメソッドを実装します :
- forward,
- inverse そして
- (少なくとも) forward_log_det_jacobian と inverse_log_det_jacobian の一つ。
これらの構成要素で、分布を変換して結果から依然としてサンプルと log prob を得ることができます！
数学では、幾分雑です :
- $X$ は pdf $p(x)$ を持つ確率変数です。
- $g$ は $X$ の空間上の滑らかで、可逆な関数です。
- $Y = g(X)$ は新しい、変換された確率変数です。
- $p(Y=y) = p(X=g^{-1}(y)) \cdot |\nabla g^{-1}(y)|$
キャッシング
Bijector はまた forward と inverse 計算、そして log-det-Jacobians をキャッシュもします、これは潜在的に非常に高価な演算を繰り返すことを節約します！
print_subclasses_from_module(tfp.bijectors, tfp.bijectors.Bijector)
AbsoluteValue, Affine, AffineLinearOperator, AffineScalar, BatchNormalization Bijector, Blockwise, Chain, CholeskyOuterProduct, CholeskyToInvCholesky CorrelationCholesky, Cumsum, DiscreteCosineTransform, Exp, Expm1, FFJORD FillScaleTriL, FillTriangular, FrechetCDF, GeneralizedPareto, GumbelCDF Identity, Inline, Invert, IteratedSigmoidCentered, KumaraswamyCDF, LambertWTail Log, Log1p, MaskedAutoregressiveFlow, MatrixInverseTriL, MatvecLU, NormalCDF Ordered, Pad, Permute, PowerTransform, RationalQuadraticSpline, RealNVP Reciprocal, Reshape, Scale, ScaleMatvecDiag, ScaleMatvecLU ScaleMatvecLinearOperator, ScaleMatvecTriL, ScaleTriL, Shift, Sigmoid SinhArcsinh, SoftClip, Softfloor, SoftmaxCentered, Softplus, Softsign, Square Tanh, TransformDiagonal, Transpose, WeibullCDF
単純な Bijector
normal_cdf = tfp.bijectors.NormalCDF() xs = np.linspace(-4., 4., 200) plt.plot(xs, normal_cdf.forward(xs)) plt.show()
plt.plot(xs, normal_cdf.forward_log_det_jacobian(xs, event_ndims=0)) plt.show()
分布を変換する Bijector
exp_bijector = tfp.bijectors.Exp() log_normal = exp_bijector(tfd.Normal(0., .5)) samples = log_normal.sample(1000) xs = np.linspace(1e-10, np.max(samples), 200) sns.distplot(samples, norm_hist=True, kde=False) plt.plot(xs, log_normal.prob(xs), c='k', alpha=.75) plt.show()
Bijector をバッチ化する
# Create a batch of bijectors of shape [3,]
softplus = tfp.bijectors.Softplus(
hinge_softness=[1., .5, .1])
print("Hinge softness shape:", softplus.hinge_softness.shape)
Hinge softness shape: (3,)
# For broadcasting, we want this to be shape [200, 1]
xs = np.linspace(-4., 4., 200)[..., np.newaxis]
ys = softplus.forward(xs)
print("Forward shape:", ys.shape)
Forward shape: (200, 3)
# Visualization
lines = plt.plot(np.tile(xs, 3), ys)
for line, hs in zip(lines, softplus.hinge_softness):
line.set_label("Softness: %1.1f" % hs)
plt.legend()
plt.show()
キャッシング
# This bijector represents a matrix outer product on the forward pass, # and a cholesky decomposition on the inverse pass. The latter costs O(N^3)! bij = tfb.CholeskyOuterProduct() size = 2500 # Make a big, lower-triangular matrix big_lower_triangular = tf.eye(size) # Squaring it gives us a positive-definite matrix big_positive_definite = bij.forward(big_lower_triangular) # Caching for the win! %timeit bij.inverse(big_positive_definite) %timeit tf.linalg.cholesky(big_positive_definite)
10000 loops, best of 3: 112 µs per loop 1 loops, best of 3: 437 ms per loop
MCMC
TFP はハミルトニアン・モンテカルロを含む、幾つかの標準的なマルコフ連鎖モンテカルロ・アルゴリズムのための組込みサポートを持ちます。
データセットを生成する
# Generate some data def f(x, w): # Pad x with 1's so we can add bias via matmul x = tf.pad(x, [[1, 0], [0, 0]], constant_values=1) linop = tf.linalg.LinearOperatorFullMatrix(w[..., np.newaxis]) result = linop.matmul(x, adjoint=True) return result[..., 0, :] num_features = 2 num_examples = 50 noise_scale = .5 true_w = np.array([-1., 2., 3.]) xs = np.random.uniform(-1., 1., [num_features, num_examples]) ys = f(xs, true_w) + np.random.normal(0., noise_scale, size=num_examples)
# Visualize the data set plt.scatter(*xs, c=ys, s=100, linewidths=0) grid = np.meshgrid(*([np.linspace(-1, 1, 100)] * 2)) xs_grid = np.stack(grid, axis=0) fs_grid = f(xs_grid.reshape([num_features, -1]), true_w) fs_grid = np.reshape(fs_grid, [100, 100]) plt.colorbar() plt.contour(xs_grid[0, ...], xs_grid[1, ...], fs_grid, 20, linewidths=1) plt.show()
結合 log-prob 関数を定義する
非正規化事後 (分布) は結合 log prob の部分適用を形成するためデータに渡り閉じた結果です。
# Define the joint_log_prob function, and our unnormalized posterior.
def joint_log_prob(w, x, y):
# Our model in maths is
# w ~ MVN([0, 0, 0], diag([1, 1, 1]))
# y_i ~ Normal(w @ x_i, noise_scale), i=1..N
rv_w = tfd.MultivariateNormalDiag(
loc=np.zeros(num_features + 1),
scale_diag=np.ones(num_features + 1))
rv_y = tfd.Normal(f(x, w), noise_scale)
return (rv_w.log_prob(w) +
tf.reduce_sum(rv_y.log_prob(y), axis=-1))
# Create our unnormalized target density by currying x and y from the joint. def unnormalized_posterior(w): return joint_log_prob(w, xs, ys)
HMC TransitionKernel を構築して sample_chain を呼び出す
# Create an HMC TransitionKernel hmc_kernel = tfp.mcmc.HamiltonianMonteCarlo( target_log_prob_fn=unnormalized_posterior, step_size=np.float64(.1), num_leapfrog_steps=2)
# We wrap sample_chain in tf.function, telling TF to precompile a reusable
# computation graph, which will dramatically improve performance.
@tf.function
def run_chain(initial_state, num_results=1000, num_burnin_steps=500):
return tfp.mcmc.sample_chain(
num_results=num_results,
num_burnin_steps=num_burnin_steps,
current_state=initial_state,
kernel=hmc_kernel,
trace_fn=lambda current_state, kernel_results: kernel_results)
initial_state = np.zeros(num_features + 1)
samples, kernel_results = run_chain(initial_state)
print("Acceptance rate:", kernel_results.is_accepted.numpy().mean())
Acceptance rate: 0.913
That’s not great! We’d like an acceptance rate closer to .65.
(“Optimal Scaling for Various Metropolis-Hastings Algorithms“, Roberts & Rosenthal, 2001 参照)
適応可能な (= adaptive) ステップサイズ
HMC TransitionKernel を SimpleStepSizeAdaptation “meta-kernel” にラップすることができます、これは burnin の間に HMC ステップサイズを適応させるために幾つかの (むしろ単純なヒューリスティックな) ロジックを適用します。ステップサイズを適応させるために burnin の 80% を割り当ててから残りの 20% を単に混合 (= mixing) に進ませます。
# Apply a simple step size adaptation during burnin
@tf.function
def run_chain(initial_state, num_results=1000, num_burnin_steps=500):
adaptive_kernel = tfp.mcmc.SimpleStepSizeAdaptation(
hmc_kernel,
num_adaptation_steps=int(.8 * num_burnin_steps),
target_accept_prob=np.float64(.65))
return tfp.mcmc.sample_chain(
num_results=num_results,
num_burnin_steps=num_burnin_steps,
current_state=initial_state,
kernel=adaptive_kernel,
trace_fn=lambda cs, kr: kr)
samples, kernel_results = run_chain(
initial_state=np.zeros(num_features+1))
print("Acceptance rate:", kernel_results.inner_results.is_accepted.numpy().mean())
Acceptance rate: 0.431
# Trace plots
colors = ['b', 'g', 'r']
for i in range(3):
plt.plot(samples[:, i], c=colors[i], alpha=.3)
plt.hlines(true_w[i], 0, 1000, zorder=4, color=colors[i], label="$w_{}$".format(i))
plt.legend(loc='upper right')
plt.show()
# Histogram of samples
for i in range(3):
sns.distplot(samples[:, i], color=colors[i])
ymax = plt.ylim()[1]
for i in range(3):
plt.vlines(true_w[i], 0, ymax, color=colors[i])
plt.ylim(0, ymax)
plt.show()
診断
追跡プロットは素晴らしいですが、診断 (= diagnostics) は更に素晴らしいです！
最初に多重連鎖を実行する必要があります。これは initial_state tensor のバッチを与えるくらいに単純です。
# Instead of a single set of initial w's, we create a batch of 8.
num_chains = 8
initial_state = np.zeros([num_chains, num_features + 1])
chains, kernel_results = run_chain(initial_state)
r_hat = tfp.mcmc.potential_scale_reduction(chains)
print("Acceptance rate:", kernel_results.inner_results.is_accepted.numpy().mean())
print("R-hat diagnostic (per latent variable):", r_hat.numpy())
Acceptance rate: 0.6545 R-hat diagnostic (per latent variable): [1.00242245 1.00223462 1.00079678]
ノイズスケールのサンプリング
# Define the joint_log_prob function, and our unnormalized posterior.
def joint_log_prob(w, sigma, x, y):
# Our model in maths is
# w ~ MVN([0, 0, 0], diag([1, 1, 1]))
# y_i ~ Normal(w @ x_i, noise_scale), i=1..N
rv_w = tfd.MultivariateNormalDiag(
loc=np.zeros(num_features + 1),
scale_diag=np.ones(num_features + 1))
rv_sigma = tfd.LogNormal(np.float64(1.), np.float64(5.))
rv_y = tfd.Normal(f(x, w), sigma[..., np.newaxis])
return (rv_w.log_prob(w) +
rv_sigma.log_prob(sigma) +
tf.reduce_sum(rv_y.log_prob(y), axis=-1))
# Create our unnormalized target density by currying x and y from the joint.
def unnormalized_posterior(w, sigma):
return joint_log_prob(w, sigma, xs, ys)
# Create an HMC TransitionKernel
hmc_kernel = tfp.mcmc.HamiltonianMonteCarlo(
target_log_prob_fn=unnormalized_posterior,
step_size=np.float64(.1),
num_leapfrog_steps=4)
# Create a TransformedTransitionKernl
transformed_kernel = tfp.mcmc.TransformedTransitionKernel(
inner_kernel=hmc_kernel,
bijector=[tfb.Identity(), # w
tfb.Invert(tfb.Softplus())]) # sigma
# Apply a simple step size adaptation during burnin
@tf.function
def run_chain(initial_state, num_results=1000, num_burnin_steps=500):
adaptive_kernel = tfp.mcmc.SimpleStepSizeAdaptation(
transformed_kernel,
num_adaptation_steps=int(.8 * num_burnin_steps),
target_accept_prob=np.float64(.75))
return tfp.mcmc.sample_chain(
num_results=num_results,
num_burnin_steps=num_burnin_steps,
current_state=initial_state,
kernel=adaptive_kernel,
trace_fn=lambda cs, kr: kr)
# Instead of a single set of initial w's, we create a batch of 8.
num_chains = 8
initial_state = [np.zeros([num_chains, num_features + 1]),
.54 * np.ones([num_chains], dtype=np.float64)]
chains, kernel_results = run_chain(initial_state)
r_hat = tfp.mcmc.potential_scale_reduction(chains)
print("Acceptance rate:", kernel_results.inner_results.inner_results.is_accepted.numpy().mean())
print("R-hat diagnostic (per w variable):", r_hat[0].numpy())
print("R-hat diagnostic (sigma):", r_hat[1].numpy())
Acceptance rate: 0.75925 R-hat diagnostic (per w variable): [1.0009577 1.00296745 1.00102734] R-hat diagnostic (sigma): 1.0042382564283778
w_chains, sigma_chains = chains
# Trace plots of w (one of 8 chains)
colors = ['b', 'g', 'r', 'teal']
fig, axes = plt.subplots(4, num_chains, figsize=(4 * num_chains, 8))
for j in range(num_chains):
for i in range(3):
ax = axes[i][j]
ax.plot(w_chains[:, j, i], c=colors[i], alpha=.3)
ax.hlines(true_w[i], 0, 1000, zorder=4, color=colors[i], label="$w_{}$".format(i))
ax.legend(loc='upper right')
ax = axes[3][j]
ax.plot(sigma_chains[:, j], alpha=.3, c=colors[3])
ax.hlines(noise_scale, 0, 1000, zorder=4, color=colors[3], label=r"$\sigma$".format(i))
ax.legend(loc='upper right')
fig.tight_layout()
plt.show()
# Histogram of samples of w
fig, axes = plt.subplots(4, num_chains, figsize=(4 * num_chains, 8))
for j in range(num_chains):
for i in range(3):
ax = axes[i][j]
sns.distplot(w_chains[:, j, i], color=colors[i], norm_hist=True, ax=ax, hist_kws={'alpha': .3})
for i in range(3):
ax = axes[i][j]
ymax = ax.get_ylim()[1]
ax.vlines(true_w[i], 0, ymax, color=colors[i], label="$w_{}$".format(i), linewidth=3)
ax.set_ylim(0, ymax)
ax.legend(loc='upper right')
ax = axes[3][j]
sns.distplot(sigma_chains[:, j], color=colors[3], norm_hist=True, ax=ax, hist_kws={'alpha': .3})
ymax = ax.get_ylim()[1]
ax.vlines(noise_scale, 0, ymax, color=colors[3], label=r"$\sigma$".format(i), linewidth=3)
ax.set_ylim(0, ymax)
ax.legend(loc='upper right')
fig.tight_layout()
plt.show()
There’s a lot more!
これらのクールなブログ投稿とサンプルを調べてください :
- 構造時系列サポート ブログ colab
- 確率的 Keras 層 (入力: Tensor, 出力: Distribution!) ブログ colab
- ガウス過程回帰 colab と潜在変数モデリング colab
GitHub 上のより多くのサンプルとノートブックは ここ です！
以上
人工知能開発支援
◆ クラスキャットは 人工知能研究開発支援 サービスを提供しています :
(株)クラスキャット
セールス・インフォメーション
E-Mail：sales-info@classcat.com
- テクニカルコンサルティングサービス
- 実証実験 (プロトタイプ構築)
- アプリケーションへの実装
- 人工知能研修サービス
(株)クラスキャット
セールス・インフォメーション
E-Mail：sales-info@classcat.com