PyMC5はバックエンドでPyTensorが使われています。単純なモデルを記述する分にはPyTensorを意識しなくても扱えますが、モデルが複雑化していくとPyTensorに触れる必要が生じてきます。PyMC5を扱う上で知っておくと便利なPyTensorの基礎的な機能について紹介します。
PyMCの公式でもPyTensorの解説はされているのでこちらも合わせて確認してください。
https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/pymc_pytensor.html
PyTensorのチュートリアルも用意されています。
https://pytensor.readthedocs.io/en/latest/tutorial/index.html
PyTensorとは
PyTensorは、多次元配列を含む数学的な式を定義し、最適化し、効率的に評価することができるPythonライブラリです。
もう少し、簡単に説明すると、NumPyのように多次元配列の数学的演算を行う機能があります。また、グラフフレームワークにより数式を最適化して効率的に計算することができるようです。
インストール
インストールは他のPythonライブラリと同じです。
pip install pytensor
PyMC5をインストールしている場合は一緒にインストールされています。
PyTensorは以下のように読み込みます。
import pytensor.tensor as pt
PyTensorの基礎的な操作
基本的にはNumPyと同じような感覚で使えます。
配列の生成
NumPyのnp.ones()のようにPyTensorでも要素が全て1の配列を作ることができます。
array1 = pt.ones(shape=(2,2))
これで、2×2のTensorVariableが作れます。
これはNumPyの配列に似ていますが、PyTensorの場合、以下を実行しても配列の中身を見ることはできません。
print(array1)
出力
Alloc.0
配列の中身を確認したい場合はeval()メソッドを使います。
print(array1.eval())
出力
[[1. 1.]
[1. 1.]]
データ型などはtypeで確認できます。
print(array1.type)
出力
Matrix(float64, shape=(2, 2))
reshape
これもNumPyのreshape()と使い方は一緒です。
array2 = array1.reshape((1,-1))
print(array2.eval())
出力
[[1. 1. 1. 1.]]
NumPyとの合わせ技
NumPyと組み合わせることでより自由に配列を作ることが出来ます。
まず、通常のNumPyアレイを作ります。
import numpy as np
np_arr = np.array([1,2,3])
print(np_arr)
出力
[1 2 3]
これを次のようにするとPyTensorの配列(ここではベクトル)に変換できます
array3 = pt.as_tensor_variable(np_arr)
print(array3.eval())
print(array3.type)
出力
[1 2 3]
Vector(int64, shape=(3,))
また、NumPyアレイの形式で要素番号を指定した配列生成し、これを使用するとかなり自由にPytensorの配列を変形させることが出来ます。
# 要素番号の配列を生成
idx_np = np.array([[0,2,1],[2,0,1],[2,0,1]])
# PyTensorを変形
array4 = array3[idx_np]
print(array4.eval())
print(array4.type)
出力
[[1 3 2]
[3 1 2]
[3 1 2]]
Matrix(int64, shape=(?, ?))
PyMC5での活用
これがPyMC5でモデリングするときに非常に役立ちます。PyMC5で設定するパラメータはPyTensorの配列となっているので上で紹介したメソッドなどを使用することが出来ます。
例えばPyMC5でモデリングするとき次のように記述するとします。
import pymc as pm
# モデリング
class MyModel(pm.Model):
def __init__(self, *args, **kwargs) -> None:
super().__init__()
self.param1_mu = pm.Normal('param1_mu', mu = 0, sigma = 10, shape = (2,2))
self.param1_sd = pm.Uniform('param1_sd', lower = 0, upper = 100)
self.param1 = pm.Normal('param1', mu = self.param1_mu, sigma = self.param1_sd)
mod = MyModel()
ここでハイパーパラメータparam1_muがshapeで設定した通りの配列になっているかeval()とtypeを使って確認してみます。
eval()とtypeはモデルのデバッグで非常に役立ちます。
print('正規分布から発生した乱数の配列(param1_mu)')
print(mod.param1_mu.eval())
print('\nparam1_muのデータ型の確認')
print(mod.param1_mu.type)
設定どおりの配列になっています。
正規分布から発生した乱数の配列(param1_mu)
[[ 4.8983854 13.04799909]
[ -6.75311399 -14.84662657]]
param1_muのデータ型の確認
Matrix(float64, shape=(2, 2))
ちなみに
print(mod.param1_mu)
とすると、param1_muとだけ返されます。
次に、シェープサイズが異なる2つのハイパーパラメータを持つparam1の配列がどうなっているのか確認してみます。
print('正規分布から発生した乱数の配列(param1)')
print(mod.param1.eval())
print('\nparam1のデータ型の確認')
print(mod.param1.type)
出力
正規分布から発生した乱数の配列(param1)
[[ 18.33583847 -45.5451446 ]
[ 41.19690527 33.46987694]]
param1のデータ型の確認
Matrix(float64, shape=(2, 2))
2×2の配列になっていました。