PyMC5を扱う上で知っておくと便利なPyTensorの機能

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は以下のように読み込みます。

``` python
import pytensor.tensor as pt
```

# PyTensorの基礎的な操作

基本的にはNumPyと同じような感覚で使えます。

## 配列の生成

NumPyの`np.ones()`のようにPyTensorでも要素が全て1の配列を作ることができます。

``` python
array1 = pt.ones(shape=(2,2))
```
これで、2×2の`TensorVariable`が作れます。

これはNumPyの配列に似ていますが、PyTensorの場合、以下を実行しても配列の中身を見ることはできません。

``` python
print(array1)
```

出力

```
Alloc.0
```

配列の中身を確認したい場合は`eval()`メソッドを使います。

``` python
print(array1.eval())
```

出力

```
[[1. 1.]
 [1. 1.]]
```

データ型などは`type`で確認できます。

``` python
print(array1.type)
```

出力

```
Matrix(float64, shape=(2, 2))
```

## reshape

これもNumPyの`reshape()`と使い方は一緒です。

``` python
array2 = array1.reshape((1,-1))
print(array2.eval())
```

出力

```
[[1. 1. 1. 1.]]
```

## NumPyとの合わせ技

NumPyと組み合わせることでより自由に配列を作ることが出来ます。

まず、通常のNumPyアレイを作ります。

``` python
import numpy as np
np_arr = np.array([1,2,3])
print(np_arr)
```

出力

```
[1 2 3]
```

これを次のようにするとPyTensorの配列(ここではベクトル)に変換できます

``` python
array3 = pt.as_tensor_variable(np_arr)
print(array3.eval())
print(array3.type)
```

出力

```
[1 2 3]
Vector(int64, shape=(3,))
```

また、NumPyアレイの形式で要素番号を指定した配列生成し、これを使用するとかなり自由にPytensorの配列を変形させることが出来ます。

``` python
# 要素番号の配列を生成
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でモデリングするとき次のように記述するとします。

``` python
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`はモデルのデバッグで非常に役立ちます。

``` python
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))
```

ちなみに

``` python
print(mod.param1_mu)
```

とすると、`param1_mu`とだけ返されます。

次に、シェープサイズが異なる2つのハイパーパラメータを持つ`param1`の配列がどうなっているのか確認してみます。

``` python
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の配列になっていました。