からっぽのしょこ

読んだら書く!書いたら読む!読書読読書読書♪同じ事は二度調べ(たく)ない

【Python】3.5:線形回帰の例【緑ベイズ入門のノート】

はじめに

 『ベイズ推論による機械学習入門』の学習時のノートです。基本的な内容は「数式の行間を読んでみた」とそれを「Rで組んでみた」になります。「数式」と「プログラム」から理解するのが目標です。

 この記事は3.5節の内容です。ベイズ推論を用いて線形回帰モデルのパラメータの学習と未観測データの予測をPythonで実装します。

 省略してある内容等ありますので、本と併せて読んでください。初学者な自分が理解できるレベルまで落として書き下していますので、分かる人にはかなりくどくなっています。同じような立場の人のお役に立てれば幸いです。

【数式読解編】

www.anarchive-beta.com

【前節の内容】

www.anarchive-beta.com

【他の節一覧】

www.anarchive-beta.com

【この節の内容】

・Pythonでやってみよう

 人工的に生成したデータを用いて、パラメータを推定してみましょう。

 利用するライブラリを読み込みます。

# 利用ライブラリ
import numpy as np
import matplotlib.pyplot as plt


・観測モデルの設定

 この例では、入力値$x_n$に対して$\mathbf{x}_n = (1, x_n, x_n^2, x_n^3)$としたベクトルを用います。このとき次元数$M = 4$であり、観測モデル$p(\mathbf{y} | \mathbf{X} ,\mathbf{w})$のパラメータを$\mathbf{w} = (w_1, w_2, \cdots, w_4)$とします。また平均パラメータ$0$、分散パラメータ$\sigma^2 = \lambda^{-1}$の1次元のガウス分布に従うノイズ成分$\epsilon_n$を加えた

$$ \begin{align} y_n &= \mathbf{w}^{\top} \mathbf{x}_n + \epsilon_n \\ &= w_1 + w_2 x_n + w_3 x_n^2 + w_4 x_n^3 + \epsilon_n \tag{3.141} \end{align} $$

を出力値$y_n$とします。ここで$\lambda$を精度パラメータと呼びます。

 まずは、$x_n$から$\mathbf{x}_n$を作成する関数を定義しておきます。

# xのベクトル作成関数を定義
def x_vector(x_n, M):
    x_mn = np.zeros(shape=(M, len(x_n)))
    for m in range(M):
        x_mn[m, :] = np.power(x_n, m)
    return x_mn

# 確認
x_vector(np.array([2.0]), 5)
array([[ 1.],
       [ 2.],
       [ 4.],
       [ 8.],
       [16.]])

 $\mathbf{x}_n$は縦ベクトルであるため、$M$行1列の配列を返します。また$x^0 = 1$であることを利用しています。

 次に、観測モデルのパラメータを設定します。

# 観測モデルのパラメータを指定
M_truth = 4
sigma = 1.0
lmd = 1.0 / sigma**2
w_m = np.random.choice(
    np.arange(-1.0, 1.0, step=0.1), size=M_truth, replace=True
).reshape([M_truth, -1])

# 確認
print(w_m)
[[-0.7]
 [-0.4]
 [-0.2]
 [ 0.1]]

 w_mも$M$行1列の配列で縦ベクトルを表現します。

 観測モデルを設定できたので、グラフを作成してモデルを確認しましょう。

# 作図用のx軸の値
x_line = np.arange(-3.0, 3.0, step=0.01)
y_line = np.dot(w_m.T, x_vector(x_line, M_truth)).flatten()

# ノイズを含まない観測モデルを作図
fig = plt.figure(figsize=(12, 8))
plt.plot(x_line, y_line)
plt.suptitle('Observation Model', fontsize=20)
plt.title('M=' + str(M_truth), loc='left', fontsize=20)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

f:id:anemptyarchive:20201114180209p:plain
観測モデル

 設定した観測モデルに従って乱数(入力データ)を生成します。

# データをサンプリング
N = 50
smp_x_n = np.random.choice(
    np.arange(min(x_line), max(x_line), step=0.01), size=N, replace=True
)
y_1n = np.dot(w_m.T, x_vector(smp_x_n, M_truth))
y_1n += np.random.normal(loc=0.0, scale=1 / lmd, size=N) # ノイズ成分

# 確認
print(y_1n.shape)
(1, 50)

 データの生成に対して分布の仮定をおいていないため、np.random.choice()で一様な確率に従いランダムに値を生成します。ただし$x_n$がとり得る範囲は、グラフの描画範囲にしておきます。サンプリングしたデータsmp_x_nは、計算に用いないためnp.random.choice()の返り値のままの1次元配列です。

 観測モデル(3.141)に従い、出力値$\mathbf{y} = \{y_1, \cdots, y_N\}$を計算します。計算結果は、1行$N$列の配列(横ベクトル)になります。

 観測データの散布図を先ほどの観測モデルのグラフと重ねて確認しておきます。

# 観測データの散布図を作成
fig = plt.figure(figsize=(12, 8))
plt.scatter(smp_x_n, y_1n) # 観測データ
plt.plot(x_line, y_line) # 観測モデル
plt.suptitle('Observation Model', fontsize=20)
plt.title('M=' + str(M_truth) + ', N=' + str(N), loc='left', fontsize=20)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

f:id:anemptyarchive:20201114180239p:plain
観測データ

 以上が観測モデルに関する設定です。

 続いて、事前分布$p(\mathbf{w})$のパラメータを設定します。

# 事前分布のパラメータを指定
M = 4
m_m = np.zeros((M, 1))
sigma_mm = np.identity(M)
lambda_mm = np.linalg.inv(sigma_mm**2)

# xのベクトル作成
x_mn = x_vector(smp_x_n, M)
x_mline = x_vector(x_line, M)

# 確認
print(x_mn.shape)
(4, 50)

 事前分布は多次元ガウス分布を仮定するのでした。次元数をMに指定します。この値はパラメータ$\mathbf{w}$の要素数に対応します。

 この値に合うように、平均パラメータ$\mathbf{m} = (m_1, \cdots, m_M)$、精度行列パラメータ$\boldsymbol{\Lambda}^{-1} = \boldsymbol{\Sigma} = (\sigma_{11}^2, \cdots, \sigma_{MM}^2)$を指定する必要があります。
 この例では、平均m_mを全ての要素が0の$M$行1列の配列(縦ベクトル)、標準偏差と相関係数の行列$(\sigma_{11}, \cdots, \sigma_{MM})$を$M$行$M$列の単位行列としてその2乗の逆行列を精度行列lambda_mmとします。$M$行$M$列の単位行列はnp.identity(M)で作成し、逆行列はnp.linalg.inv()で計算します。この場合lambda_mmも$M$行$M$列の単位行列になるので、直接指定しても問題ありません。

 またサンプリングしたデータ$\{x_1, \cdots, x_N\}$をもとに、入力値$\mathbf{X} = \{\mathbf{x}_1, \cdots, \mathbf{x}_N\}$に変換します。

 設定した事前分布から$\mathbf{w}$をサンプリングしたモデルを確認しましょう。

# サンプリング数を指定
smp_size = 5

# 事前分布からwをサンプリング
smp_model_arr = np.empty((smp_size, len(x_line)))
for i in range(smp_size):
    # wをサンプリング
    smp_w_m = np.random.multivariate_normal(
        mean=m_m.flatten(), cov=np.linalg.inv(lambda_mm), size=1
    ).reshape(M, -1)
    
    # 出力値を計算
    tmp_y_line = np.dot(smp_w_m.T, x_mline).flatten()
    smp_model_arr[i] = tmp_y_line.copy()

# 事前分布からサンプリングしたモデルを作図
fig = plt.figure(figsize=(12, 8))
for i in range(smp_size):
    plt.plot(x_line, smp_model_arr[i], label=str(i+1)) # 事前分布からサンプリングしたwによるモデル
plt.plot(x_line, y_line, color='blue', linestyle='dotted', label='model') # 観測モデル
plt.title('Sampling from Piror Distribution', fontsize=20)
plt.xlabel('x')
plt.ylabel('y')
plt.ylim(min(min(y_1n.flatten()), min(y_line)), max(max(y_1n.flatten()), max(y_line)))
plt.grid()
plt.show()

f:id:anemptyarchive:20201114235541p:plain
事前分布からサンプリングしたパラメータによるモデル

 この時点では当然観測モデルとは程遠いグラフになります。

・事後分布の更新

 観測データy_nx_mnを用いて、事後分布$p(\mathbf{w} | \mathbf{y}, \mathbf{X})$(のパラメータ)を求めます。

# 事後分布のパラメータを計算
lambda_hat_mm = lmd * np.dot(x_mn, x_mn.T) + lambda_mm
tmp_m_m = lmd * np.dot(x_mn, y_1n.T)
tmp_m_m += np.dot(lambda_mm, m_m)
m_hat_m = np.dot(np.linalg.inv(lambda_hat_mm), tmp_m_m)

 事後分布(多次元ガウス分布)の精度行列パラメータ$\hat{\boldsymbol{\Lambda}}^{-1} = \hat{\boldsymbol{\Sigma}} = (\hat{\sigma}_{11}^2, \cdots, \hat{\sigma}_{MM}^2)$、平均パラメータ$\hat{\mathbf{m}} = (\hat{m}_1, \cdots, \hat{m}_M)$を、次の式で計算します。

$$ \begin{align} \hat{\boldsymbol{\Lambda}} &= \lambda \sum_{n=1}^N \mathbf{x}_n \mathbf{x}_n^{\top} + \boldsymbol{\Lambda} \\ \hat{\mathbf{m}} &= \hat{\boldsymbol{\Lambda}}^{-1} \left( \lambda \sum_{n=1}^N y_n \mathbf{x}_n + \boldsymbol{\Lambda} \mathbf{m} \right) \tag{3.148} \end{align} $$

 ただし$\sum_{n=1}^N y_n \mathbf{x}_n$の計算を効率よく処理するために、$\mathbf{X} \mathbf{y}^{\top}$で計算しています。計算結果をそれぞれlambda_hat_mmm_hat_mとします。

 事後分布からも$\mathbf{w}$をサンプリングしてモデルを確認しましょう。

# 事後分布からwをサンプリング
smp_model_arr = np.empty((smp_size, len(x_line)))
for i in range(smp_size):
    # wをサンプリング
    smp_w_m = np.random.multivariate_normal(
        mean=m_hat_m.flatten(), cov=np.linalg.inv(lambda_hat_mm), size=1
    ).reshape(M, -1)
    
    # 出力値を計算
    tmp_y_line = np.dot(smp_w_m.T, x_mline).flatten()
    smp_model_arr[i] = tmp_y_line.copy()

# 事後分布からサンプリングしたモデルを作図
fig = plt.figure(figsize=(12, 8))
for i in range(smp_size):
    plt.plot(x_line, smp_model_arr[i], label=str(i+1)) # 事後分布からサンプリングしたwによるモデル
plt.scatter(smp_x_n, y_1n.flatten()) # 観測データ
plt.plot(x_line, y_line, color='blue', linestyle='dotted', label='model') # 観測モデル
plt.title('Sampling from Posterior Distribution', fontsize=20)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

f:id:anemptyarchive:20201114180329p:plain
事後分布からサンプリングしたパラメータによるモデル

 観測データから学習(パラメータを更新)することで、どれも観測モデルに近い形になっていることが分かります。

・予測分布の更新

 最後に、事後分布パラメータm_mlambda_hat_mmを用いて、予測分布$(y_{*} | \mathbf{x}_{*}, \mathbf{y}, \mathbf{X})$(のパラメータ)を求めます。

# 予測分布のパラメータを計算
sigma2_star_hat_line = np.repeat(1 / lmd, len(x_line))
for i in range(len(x_line)):
    sigma2_star_hat_line[i] += x_mline[:, i].T.dot(np.linalg.inv(lambda_hat_mm)).dot(x_mline[:, i])
mu_star_hat_line = np.dot(m_hat_m.T, x_mline).flatten()

 予測分布(1次元のガウス分布)の平均パラメータ$\hat{\mu}_{*}$、精度パラメータ$\hat{\lambda}_{*}^{-1} = \hat{\sigma}^2$を、次の式で計算します。

$$ \begin{aligned} \hat{\mu}_{*} &= \hat{\mathbf{m}}^{\top} \mathbf{x}_{*} \\ \hat{\lambda}_{*}^{-1} &= \lambda^{-1} + \mathbf{x}_{*}^{\top} \hat{\boldsymbol{\Lambda}}^{-1} \mathbf{x}_{*} \end{aligned} $$

 グラフの描画範囲全体の未知の入力データ$\mathbf{X}_{*}$に対して計算する必要があるため、x_mlineを使います。計算結果をそれぞれsigma2_star_hat_linemu_star_hat_lineとします。

 平均mu_star_hat_lineと平均から標準偏差sqrt(sigma2_star_hat_line)を足し引きしたものとあわせて予測分布を作図します。

# 予測分布を作図
fig = plt.figure(figsize=(12, 8))
plt.plot(x_line, mu_star_hat_line, color='orange', label='predict') # 予測分布の期待値
plt.plot(x_line, mu_star_hat_line + np.sqrt(sigma2_star_hat_line), 
         color='#00A968', linestyle='--', label='$+\sigma$') # +sigma
plt.plot(x_line, mu_star_hat_line - np.sqrt(sigma2_star_hat_line), 
         color='#00A968', linestyle='--', label='$-\sigma$') # -sigma
plt.plot(x_line, y_line, linestyle=':', color='blue', label='model') # 観測モデル
plt.scatter(smp_x_n, y_1n.flatten(), color='chocolate') # 観測データ
plt.suptitle('Predict Distribution', fontsize=20)
plt.title('M=' + str(M) + ', N=' + str(N), loc='left', fontsize=20)
plt.xlabel('x')
plt.ylabel('y')
plt.ylim(min(min(y_1n.flatten()), min(y_line)), max(max(y_1n.flatten()), max(y_line)))
plt.grid()
plt.legend()
plt.show()

f:id:anemptyarchive:20201114180351p:plain
予測分布

 以上が推論処理です。

 いくつか値を変えて試してみましょう。

 データ数$N$を減らしてパラメータ数$M$を増やした場合は、次のようになります。

f:id:anemptyarchive:20201114180448p:plain
予測分布:過学習

 モデルへの当てはまりが悪くなり、またx軸の各値に対するy軸の値の標準偏差(緑の破線の縦方向の範囲)も広くなっているのが分かります。

 次に描画範囲を広くして、観測データにない範囲の入力値に対する予測分布を確認してみましょう。

f:id:anemptyarchive:20201114180509p:plain
予測分布:過学習

 $M$の値が真の値よりも大きいと、観測データの範囲を越えたとたんに予測できていません。

f:id:anemptyarchive:20201114180526p:plain
予測分布

 $M$が真の値と一致していると、(ある程度)予測ができています。

・おまけ

 matplotlib.animationモジュールを利用して、gif画像を作成するコードです。ただしいまいち理解せずに書いたものなので参考程度にお願いします…

 観測データの増加が、予測精度へ与える影響を確認します。

・作図用のPythonコード(クリックで展開)

# 利用ライブラリ
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# xのベクトル作成関数を定義
def x_vector(x_n, M):
    x_mn = np.zeros(shape=(M, len(x_n)))
    for m in range(M):
        x_mn[m, :] = np.power(x_n, m)
    return x_mn

# 観測モデルのパラメータを指定
M_truth = 4
sigma = 1.0
lmd = 1.0 / sigma**2
w_m = np.random.choice(
    np.arange(-1.0, 1.0, step=0.1), size=M_truth, replace=True
).reshape([M_truth, -1])

# 作図用のx軸の値
x_line = np.arange(-3.0, 3.0, step=0.01)
y_line = np.dot(w_m.T, x_vector(x_line, M_truth)).flatten()
# gif画像でデータ数の変化による学習への影響を確認

# データ数を指定:(試行回数)
N = 100

# 事前分布のパラメータを指定
M = 10
m_m = np.zeros((M, 1))
sigma_mm = np.identity(M)
lambda_mm = np.linalg.inv(sigma_mm**2)

# 作図用のx軸の値
x_mline = x_vector(x_line, M)

# 推論
smp_x_n = np.array([])
y_1n = np.array([[]])
for n in range(N):
    
    # データをサンプリング
    smp_x = np.random.choice(
        np.arange(min(x_line), max(x_line), step=0.01), size=1, replace=True
    )
    smp_x_n = np.append(smp_x_n, smp_x)
    
    # 出力値を計算
    tmp_y = np.dot(w_m.T, x_vector(smp_x_n[n:n+1], M_truth))
    tmp_y += np.random.normal(loc=0.0, scale=1 / lmd, size=1) # ノイズ成分
    y_1n = np.append(y_1n, tmp_y, axis=1)
    
    # xのベクトル作成
    x_mn = x_vector(smp_x_n[n:n+1], M)
    
    # 事後分布のパラメータを更新
    old_lambda_mm = lambda_mm.copy()
    lambda_mm += lmd * np.dot(x_mn, x_mn.T)
    tmp_m_m = lmd * np.dot(x_mn, y_1n[0:1, n:n+1].T)
    tmp_m_m += np.dot(old_lambda_mm, m_m)
    m_m = np.dot(np.linalg.inv(lambda_mm), tmp_m_m)
    
    # 予測分布のパラメータを計算
    sigma2_star_line = np.repeat(1 / lmd, len(x_line))
    for i in range(len(x_line)):
        sigma2_star_line[i] += x_mline[:, i].T.dot(np.linalg.inv(lambda_mm)).dot(x_mline[:, i])
    mu_star_line = np.dot(m_m.T, x_mline).flatten()
    
    # 推移を記録
    if n == 0: # 初回
        # 2次元配列を作成
        sigma2_star_arr = np.array([sigma2_star_line])
        mu_star_arr = np.array([mu_star_line])
    elif n > 0: # 初回以降
        # 配列に追加
        sigma2_star_arr = np.append(sigma2_star_arr, [sigma2_star_line.copy()], axis=0)
        mu_star_arr = np.append(mu_star_arr, [mu_star_line.copy()], axis=0)

# 予測分布を作図

fig = plt.figure(figsize=(12, 8))
fig.suptitle('Predict Distribution', fontsize=20)
ax = fig.add_subplot(1, 1, 1)

# 作図処理を関数として定義
def update(n):
    
    # 前フレームのグラフを初期化
    ax.cla()
    
    # nフレーム目のグラフを描画
    ax.plot(x_line, mu_star_arr[n], color='orange', label='predict') # 予測分布の期待値
    ax.plot(x_line, mu_star_arr[n] + np.sqrt(sigma2_star_arr[n]), 
            color='#00A968', linestyle='--', label='$+\sigma$') # +sigma
    ax.plot(x_line, mu_star_arr[n] - np.sqrt(sigma2_star_arr[n]), 
            color='#00A968', linestyle='--', label='$-\sigma$') # -sigma
    ax.plot(x_line, y_line, linestyle=':', color='blue', label='model') # 観測モデル
    ax.scatter(smp_x_n[0:n+1], y_1n[0, 0:n+1].flatten(), color='chocolate') # 観測データ
    
    # グラフの設定
    ax.set_title('M=' + str(M) + ', N=' + str(n+1), loc='left', fontsize=20)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_ylim([min(min(y_1n.flatten()), min(y_line)), max(max(y_1n.flatten()), max(y_line))])
    ax.grid(True)
    ax.legend(loc='upper right')
    
# gif画像を作成
ani = animation.FuncAnimation(fig, update, frames=len(mu_star_arr), interval=100)
ani.save("figname.gif")

 異なる点のみを簡単に解説します。

 各データによってどのように学習する(推定値が変化する)のかを確認するため、こちらはfor文で1データずつ処理します。よって観測データ数Nがイタレーション数になります。

 一度の処理で事後分布のパラメータを計算するのではなく、事前分布のパラメータに対して繰り返し観測データの情報を与えることで更新(上書き)していきます。$\sum_{n=1}^N$の計算は、N回繰り返し計算することで実行されます。
 ただし事後分布のパラメータの計算(3.148)の計算において、更新前(事前分布)のパラメータ$\boldsymbol{\Lambda}$を使うため、old_lambda_mmとして値を一時的に保存しておきます。

f:id:anemptyarchive:20201114182744g:plain
予測分布の推移:過学習

 (わざわざ学習してるっぽい振る舞いをさせていますが、これが良くないことは先ほど確認した通りです。)

 同じ観測データを用いたとき、パラメータの数(事前分布の次元数)の違いが予測精度へ与える影響を確認します。

・作図用のPythonコード(クリックで展開)

# gif画像でパラメータ数の変化による学習への影響を確認

# データをサンプリング
N = 10
smp_x_n = np.random.choice(
    np.arange(min(x_line), max(x_line), step=0.01), size=N, replace=True
)
y_1n = np.dot(w_m.T, x_vector(smp_x_n, M_truth))
y_1n += np.random.normal(loc=0.0, scale=1 / lmd, size=N) # ノイズ成分

# パラメータ数(次元数)を指定:(試行回数)
max_M = 25

# 推論
for m in range(max_M):
    
    # 事前分布のパラメータを設定
    m_m = np.zeros((m, 1))
    sigma_mm = np.identity(m)
    lambda_mm = np.linalg.inv(sigma_mm**2)
    
    # xのベクトル作成
    x_mn = x_vector(smp_x_n, m)
    x_mline = x_vector(x_line, m)
    
    # 事後分布のパラメータを更新
    old_lambda_mm = lambda_mm.copy()
    lambda_mm += lmd * np.dot(x_mn, x_mn.T)
    tmp_m_m = lmd * np.dot(x_mn, y_1n.T)
    tmp_m_m += np.dot(old_lambda_mm, m_m)
    m_m = np.dot(np.linalg.inv(lambda_mm), tmp_m_m)
    
    # 予測分布のパラメータを計算
    sigma2_star_line = np.repeat(1 / lmd, len(x_line))
    for i in range(len(x_line)):
        sigma2_star_line[i] += x_mline[:, i].T.dot(np.linalg.inv(lambda_mm)).dot(x_mline[:, i])
    mu_star_line = np.dot(m_m.T, x_mline).flatten()
    
    # 推移を記録
    if m == 0: # 初回
        # 2次元配列を作成
        sigma2_star_arr = np.array([sigma2_star_line])
        mu_star_arr = np.array([mu_star_line])
    elif m > 0: # 初回以降
        # 配列に追加
        sigma2_star_arr = np.append(sigma2_star_arr, [sigma2_star_line.copy()], axis=0)
        mu_star_arr = np.append(mu_star_arr, [mu_star_line.copy()], axis=0)

# 予測分布を作図

fig = plt.figure(figsize=(12, 8))
fig.suptitle('Predict Distribution', fontsize=20)
ax = fig.add_subplot(1, 1, 1)

# 作図処理を関数として定義
def update(m):
    
    # 前フレームのグラフを初期化
    ax.cla()
    
    # nフレーム目のグラフを描画
    ax.plot(x_line, mu_star_arr[m], color='orange', label='predict') # 予測分布の期待値
    ax.plot(x_line, mu_star_arr[m] + np.sqrt(sigma2_star_arr[m]), 
            color='#00A968', linestyle='--', label='$+\sigma$') # +sigma
    ax.plot(x_line, mu_star_arr[m] - np.sqrt(sigma2_star_arr[m]), 
            color='#00A968', linestyle='--', label='$-\sigma$') # -sigma
    ax.plot(x_line, y_line, linestyle=':', color='blue', label='model') # 観測モデル
    ax.scatter(smp_x_n, y_1n.flatten(), color='chocolate') # 観測データ
    
    # グラフの設定
    ax.set_title('M=' + str(m+1) + ', N=' + str(N), loc='left', fontsize=20)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_ylim([min(min(y_1n.flatten()), min(y_line)), max(max(y_1n.flatten()), max(y_line))])
    ax.grid(True)
    ax.legend(loc='upper right')

# gif画像を作成
ani = animation.FuncAnimation(fig, update, frames=len(mu_star_arr), interval=200)
ani.save("figname.gif")

 こちらはMをイタレーション数として、推論処理をM回繰り返します。

f:id:anemptyarchive:20201114183846g:plain
予測分布の変化:過学習


参考文献

  • 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.

おわりに

 Rで書いたものをPythonに翻訳した形です。Pythonの勉強のためにベイズ推論を扱った、という方が近いかもしれません。まだPythonらしさのようなものを全然掴めていないので、色々アレなところがあるのだろうと思います。何か気になるところがあれば、ご指摘いただければとても嬉しいです。
 せめてclassとして実装すればPythonらしかろうと試してみましたが、この例だと$x_n$を$(x_n^0, x_n^1, \cdots, x_n^{M-1})$に変換する処理を組み込むのが面倒なうえに分かりにくくなったので止めました、、

 しかし大半はコピペで済むだろうし水増しっぽくなってしまうかなと書く前には危惧していたのですが、全然そんなことなくて普通に一記事分疲れました。これまでの記事もPythonでも実装してみたいけど、まだ手が出ないかなぁ。

【次節の内容】

www.anarchive-beta.com