ホーム » TensorFlow Probability

TensorFlow Probability」カテゴリーアーカイブ

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] を持つか注意してください。

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

  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]。

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

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!

これらのクールなブログ投稿とサンプルを調べてください :

GitHub 上のより多くのサンプルとノートブックは ここ です!

 

以上






TensorFlow Probability 0.6.0 リリースノート

TensorFlow Probability 0.6.0 リリースノート (翻訳)

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

* 本ページは、github TensorFlow Probability の releases の TensorFlow Probability 0.6.0 を翻訳したものです:

 

本文

これは TensorFlow Probability の 0.6 リリースです。それは TensorFlow 1.13.1 に対してテストされてステーブルです。

 

変更ノート

  • tfp.positive_semidefinite_kernels.RationalQuadratic を追加します。
  • tfpl.MultivariateNormalTriL で float64 をサポートします。
  • IndependentLogistic と IndependentPoisson 分布層を追加します。
  • Edward2 確率変数の値を設定するために make_value_setter interceptor を追加します。
  • LinearGaussianStateSpaceModel のメンバー関数として、Kalman Smoother の実装。
  • Bijector キャッシングは eager モードで実行するとき一つの方向でのみ有効にされます。forward(x) か inverse(y) を同じ x か y 値で繰り返し計算する場合には eager モードで何某かのパフォーマンス後退を引き起こすかもしれません。
  • tfpl.Independent{Bernoulli,Normal} で rank-0/empty event_shape を扱います。
  • eager モードで追加のテストを実行します。
  • quantiles(x, n, …) が tfp.stats に追加されました。
  • tensorflow_probability を Tensorflow 2.0 TensorShape インデキシングと互換にします。
  • Chi, Chi2 のための KL ダイバージェンスをテストするとき scipy.special 関数を使用します。
  • STS モデルから予測 (= forecast) を作成するためのメソッドを追加します。
  • MixtureSameFamily 分布を追加します。
  • Chi 分布を追加します。
  • doc typo tfp.Distribution -> tfd.Distribution を修正します。
  • Gumbel-Gumbel KL ダイバージェンスを追加します。
  • HalfNormal-HalfNormal KL ダイバージェンスを追加します。
  • Chi2-Chi2 KL ダイバージェンスのユニットテストを追加します。
  • Exponential-Exponential KL ダイバージェンスのユニットテストを追加します。
  • Normal-Normal KL ダイバージェンスのためのサンプリングテストを追加します。
  • IndependentNormal 分布層を追加します。
  • HiddenMarkovModel に posterior_marginals を追加しました。
  • Pareto-Pareto KL ダイバージェンスを追加します。
  • 構造時系列モデルのための LinearRegression コンポーネントを追加します。
  • python Dataset オブジェクト作成の間に (それを Iterator 作成時に行なう代わりに) dataset ops をグラフを追加します (あるいは Eager execution でカーネルを作成します)。
  • テキスト・メッセージ HMC ベンチマーク。
  • switching ポアソン過程を多重変化点検出のための HMM としてエンコードするサンプル・ノートブックを追加します。
  • make_simple_step_size_update_policy に num_adaptation_steps 引数を必要とします。
  • printed ベンチマーク文字列における s/eight_hmc_schools/eight_schools_hmc/。
  • Keras モデルを通して tfd.Distribution インスタンスを plumb することを可能にするために tfp.layers.DistributionLambda を追加します。
  • tfp.math.batch_interp_regular_1d_grid を追加しています。
  • in-depth サンプルを含めるために fill_triangular の記述を更新します。
  • bijector/distribution 合成, eg, tfb.Exp(tfd.Normal(0,1)) を有効にします。
  • tfp.stats.percentile に linear と midpoint interpolation が追加されました。
  • distributions をそれらが利用する bijectors のみを含むようにします。
  • tfp.math.interp_regular_1d_grid が追加されました。
  • tfp.stats.correlation が追加されました (ピアソン相関)。
  • 最近追加された Distributions を含めるために edward2 RVs のリストを更新します。
  • 連続一様分布の密度は upper エンドポイントを含みます。
  • tfp.glm.fit_sparse でバッチ化入力のためのサポートを追加します。
  • interp_regular_1d_grid が tfp.math に追加されました。
  • HiddenMarkovModel 分布を追加しました。
  • Student’s T 過程を追加します。
  • 観測が静的にスカラーであると知られているとき matrix ops を回避することにより LinearGaussianStateSpaceModel を最適化します。
  • stddev, cholesky が tfp.stats に追加されました。
  • 構造時系列モデルを変分推論と HMC でデータに fit するためのメソッドを追加します。
  • Expm1 bijector (Y = Exp(X) – 1) を追加します。
  • tfp.stats に新しい stats 名前空間 covariance と variance が追加されました。
  • 総ての利用可能な MCMC カーネルを TransformedTransitionKernel と互換にします。

 
以上

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 0.5.0 リリースノート

TensorFlow Probability 0.5.0 リリースノート (翻訳)

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

* 本ページは、github TensorFlow Probability の releases の TensorFlow Probability 0.5.0 を翻訳したものです:

 

本文

これは TensorFlow Probability の 0.5.0 リリースです。それは TensorFlow 1.12 に対してテストされてステーブルです。

 

パッケージング変更

このリリースの時点で、別々の GPU-固有のビルドをもはやパッケージ化しません。ユーザは彼らが使用することを望む TensorFlow のバージョン (CPU or GPU) を選択することができて、そして TensorFlow Probability は両者で動作します。

その結果、私達のパッケージ requirements で TensorFlow 依存性をもはや明示的にリストアップしません (何故ならばユーザがどちらのバージョンを望むかを知ることができないからです)。TFP が TensorFlow パッケージなしで、あるいはサポートされていない TensorFlow バージョンとともにインストールされる場合、import 時に ImportError を発行します。

 

Distributions & Bijectors

  • 総ての Distributions は tf.distributions から tfp.distributions へ再配置されました (TF のものは deprecated となり TF 2.0 で除去されます)。
  • 三角分布 (= Triangular distribution) を追加します。
  • ジップ分布 (= Zipf distribution) を追加します。
  • NormalCDF (訳注: 正規累積分布関数) Bijector を追加します。
  • 多変量スチューデントの t-分布を追加します。
  • RationalQuadratic カーネルを追加します。

 

ドキュメント & サンプル

  • 変分推論を使用して GLMM をどのように fit させるかを示すサンプルを追加します。
  • ガウス過程潜在変数モデルを colab に導入します。
  • ガウス過程回帰サンプルを colab に導入します。
  • GLM アルゴリズムを示して (それらのアルゴリズムを活用する) GLM についての幾つかの結果を導出する notebook を追加します。

 
以上

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

 

以上






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