読者です 読者をやめる 読者になる 読者になる

Introduction to Dirichlet Process and its Applications

機械学習 言語モデル

学習とモデルの複雑さ*1

http://img.skitch.com/20100217-eak2i585drhqnhgpp2c9te6idi.jpg
http://img.skitch.com/20100217-d2m22eyufw1r4q8m7yf816acw1.jpg

  • 混合モデルにおける混合数Kや多項式回帰での次数Mをどのようにして決めるか?
  • AICなどの情報量基準
  • CVによるパラメータの決定

Nonparametric Bayesian

  • ノンパラベイジアンは違う発想をする
  • 柔軟でないモデルは間違った推論をしてしまう
    • 柔軟でないというのは混合数「5」の混合ガウス分布とか、次数「4」の多項式回帰とか
  • もっと柔軟なモデルを作ろう
    • モデルのパラメータ数をサンプル数によって可変にしよう
    • ある意味、パラメータ数を\inftyに持っていく
  • ノンパラメトリックなモデルはモデル選択をする必要性がない

パラメトリックなモデル

  • 有限個のパラメータ集合\thetaについて考えている
  • 新たなデータを予測するときには、前のデータとは独立なことを想定している
    • p(x | \theta, \mathcal{D}) = p(x | \theta)
  • 有限個のパラメータによって、データの特性全てを記述する
  • 手に入るデータの量が限られていれば、モデルの複雑さは限定されてしまう

ノンパラメトリックなモデル

  • データの分布は有限個のパラメータでは記述しきれないと仮定している
  • しかし、無限次元のパラメータ\thetaでなら記述できると考える
  • 実際には無限個あるわけではなく、データ量によってモデルの複雑さ(パラメータ数)を変化させていく
    • これにより、柔軟なモデルを得る

ノンパラは初めて?

  • そうでもない
    • ガウス過程が実はそうだった
Gaussian Process*2

A Gaussian process defines a distribution over functions, f, where f is a function mapping some input space \mathcal{X} to \mathcal{R}.
f : \mathcal{X} \rightarrow \mathcal{R}.

Notice that f can be an infinite-dimensional quantity (e.g. if \mathcal{X} = \mathcal{R})

Let's call this distribution P(f)

Let \mathbf{f} = (f(x_1), f(x_2), \cdots, f(x_n)) be an n-dimensional vector of function values evaluated at n points x_i \in \mathcal{X}. Note f is a random variable.

Definition: P(f) is a Gaussian process if for any finite subset \{x_1, \cdots, x_n\} \subset \mathcal{X}, the marginal distribution over that finite subset P(f) has a multivariate Gaussian distribution.

余談

  • GPはサンプル数に比例してパラメータ数が増えていった
    • しかも、推定のときは(工夫しないと)O(N^3)かかる
  • DPはサンプル数の対数に比例してパラメータ数が増えていく
    • \alpha \log{N}

応用例*3

  • Dirichlet Process(DP)の応用例であるDirichlet Process Mixture(DPM)を例に取って
    • クラスタリングのクラスタ数Kを適当に決めてくれる
  • クラスタ数Kをどのようにして決めるか
    • 事後確率最大化(MAP推定)
  • p(K|X) \propto p(X, K) = \int d \theta \sum_Z p(X, Z, \theta, K) = \int d \theta \sum_Z p(X|Z, \theta ,K) p(\theta|K) p(Z,K)
    • p(X|Z, \theta ,K) => モデルに依る(例:正規分布)
    • p(\theta|K) => パラメータの事前分布
    • p(Z,K) => 割り当ての事前分布
      • DPM!!

  • 「割り当ての事前分布」とは何か?
  • データx_1x_2x_3とあったときにこの組み合わせというのは以下の5通り
Z_1 = (x_1), (x_2), (x_3) K_1 = 3
Z_2 = (x_1, x_2), (x_3) K_2 = 2
Z_3 = (x_1), (x_2, x_3) K_3 = 2
Z_4 = (x_1, x_3), (x_2) K_4 = 2
Z_5 = (x_1, x_2, x_3) K_5 = 1
  • 同時確率P(Z, K)はこの5つにどのように確率を振っていくか、というのを表わしている
  • これをChinese Restaurant Process(あとで説明)で考えると新しいクラスタを作るというのが新しいテーブルに座るということに対応していることが分かる
    • [9]には"CRP is the corresponding distribution over partitions"と書いてあって、まさに分割についての分布を表わしている

計算のイメージ

http://img.skitch.com/20100217-rcfqu6t44twtfyjt7k7tkq7544.jpg

  • 確率の振り方はChinese Restaurant Process(あとで説明)によって与える

Chinese Restaurant Process(CRP)

  • 日本語だと中華料理店過程。
  • DPの表現の仕方の一つ
  • このプロセスでやることはこんなこと*4
    • 最初のお客さんは必ずテーブル1に着席する
    • n番目のお客さんは、既に n_k人が着席しているテーブルkに確率\frac{n_k}{n - 1 + \alpha}で座り
    • 新しいテーブルに確率\frac{\alpha}{n - 1 + \alpha}で座る

Chinese Restaurant Processの性質

  • n人のお客さんがきた時点でのテーブルの数k_n\mathrm{E}[k_n | \alpha] = O(\alpha \log(n))となり、お客さんの数の対数オーダーでテーブル数が増えていく
  • 混んでるテーブルのほうがお客さんが座る確率が高い
    • (混んでるのに座ってくる人って、、、)
  • これをクラスタリングと見なすことができる。新たなデータは属しているデータの数が多いクラスタに配属される確率が高く、ある確率でまた新しいクラスタを生成する
    • 普通のクラスタ数を決め打ちでやるクラスタリングとは違って、データ数によってクラスタ数が可変していく
    • これがInfiniteと呼ばれるゆえんか。flexibleのほうがしっくりくるけど

交換可能性(exchangability)

  • 簡単に言うとどのお客さんから座っても構成する確率には関係がない、ということ
    • i.i.dとはまた違うらしい
  • Chinese Restaurant ProcessとDPとクラスタリングについての関係については[6]が分かりやすい
  • Z_ii番目のクラスタリングされた状態(隠れ変数)、K_ii番目のクラスタリングの結果のクラスタの数

なぜDPMなのか?

  • p(K|X) \propto p(X, K) = \int d \theta \sum_Z p(X, Z, \theta, K) = \int d \theta \sum_Z p(X|Z, \theta ,K) p(\theta|K) p(Z,K)
    • p(X|Z, \theta ,K) => モデルに依る(例:正規分布)
    • p(\theta|K) => パラメータの事前分布
    • p(Z,K) => 割り当ての事前分布
  • 割り当ての事前分布p(Z,K)は事前分布なので置き方は任意である
  • モデル化の方法としては以下の二つが使われる
    • Dirichlet Process Mixture(DPM)
    • Dirichlet分布p(Z|K) + Poisson分布p(K)
  • どっちを使うかは実装の面での使い勝手による
    • Dirichlet分布 + Poisson分布だとMCMCの詳細釣り合い条件を満たすのが面倒
    • DPMだとMCMCを使って推論を割と簡単に行うことができる => これはまたあとで

ここまでのまとめ

  • 有限個のパラメータでのモデリングは柔軟性に欠ける(とノンパラベイズの人たちは考える)
  • サンプル数によって、パラメータ数を変化させていく
  • ノンパラ、実は初めてではなかった => Gaussian Process
  • 応用例としてはクラスラリングのクラスタ数を決定してくれるようなことができる
    • Chinese Restaurant Process
    • 割り当ての事前分布
    • 推論はまた後で!!

これから

  • DPには様々な表現の仕方がある
    • Chinese Restaurant Processもその一つだった
  • そもそもの定義を見る
  • 必要事項
    • ディリクレ分布の定義とその性質

Dirichlet Distribution

http://img.skitch.com/20100217-ftkhrij6854q2ejn96134ynk2t.jpg

  • ディリクレ分布は初めて?
  • ディリクレ分布はr次元の単体上の確率分布と考えることができる
  • \mathbf{p}r次元のベクトルとする。ただし、任意のiについてp_j \geq 0で、\sum_{i=1}^r p_i = 1を満たすものとする
  • このときディリクレ分布を以下で定義する
    • p(\mathbf{p} | \alpha) = \mbox{Dir}(\alpha_1, \cdots, \alpha_r) = \frac{\Gamma(\sum_i \alpha_i)}{\prod_i \Gamma(\alpha_i)} \prod_{i=1}^r p_i^{\alpha_i - 1}

Dirichlet Distributionの分かりやすい例

http://img.skitch.com/20100217-t4damdw9r13ht7shnpkegfyn9r.jpg
http://img.skitch.com/20100217-d6bmbnts9efsxrjsnn7u9hufaj.jpg

  • [15]のスライドより

Property of Dirichlet Distribution

  • 多項分布と共役な分布である(有名な事実)
    • c | \mathbf{p} \sim \mbox{Multinomial(\cdot | \mathbf{p})}
    • p(c = j | \mathbf{p}) = p_jを表わしている
  • 事後分布もまた以下のディリクレ分布となる
    • p(\mathbf{p} | c = j, \alpha) = \frac{p(c = j | \mathbf{p}) p(\mathbf{p} | \alpha)}{p(c = j | \alpha)} = \mbox{Dir}(\alpha^\prime)
    • ここで\alpha_j^\prime = \alpha_j + 1l \neq jでは\alpha_l^\prime = \alpha_l

Definition of Dirichlet Process

G_0\Theta上の関する分布とし、\alphaを正の実数とする。\Thetaの任意の(可測な)分割A_1, \cdots, A_rに対して*5、ベクトル(G(A_1), \cdots, G(A_r))は確率変数のベクトルである。ここで、GGP(\alpha, G_0)に従う確率変数である。G \sim GP(\alpha, G_0)であるとは、\Thetaの任意の(可測な)分割A_1, \cdots, A_rに対して
(G(A_1), \cdots, G(A_r)) \sim \mbox{Dir}(\alpha G_0(A_1), \cdots, \alpha G_0(A_r))
とであるということ。

Relationship between GP and DP

  • ディリクレ過程は無限次元版ディリクレ分布と考えることができる
    • どこで周辺化してきてもディリクレ分布が出てくる
  • ガウス過程と非常によく似ている。
  • ガウス過程 => distribution over function
  • ディリクレ過程 => distribution over distribution
  • GP : どんな(\mathbf{x}_1, \cdots, \mathbf{x}_{\infty})と取ってきても、対応する(y_1, \cdots, y_{\infty})がガウス分布に従う
  • DP : 空間のどんな離散化(X_1, \cdots, X_{\infty})についても、対応する離散分布がディリクレ分布\mbox{Dir}(\alpha(X_1), \cdots, \alpha(X_{\infty}))に従う

What Dirichlet Process Means?

http://img.skitch.com/20100217-f3ticf5r46stcxusj4tki5tud2.jpg

  • 上の図が基底関数G_0
  • 真ん中の図が基底関数G_0\Thetaのとある分割A_1, A_2, A_3, A_4, A_5, A_6
    • G_0(A_4)は区間A_4で基底関数G_0を積分した面積に相当
  • 下の図は、区間における線分の長さの和(と考えてもたぶんよい)
  • cf. (G(A_1), \cdots, G(A_r)) \sim \mbox{Dir}(\alpha G_0(A_1), \cdots, \alpha G_0(A_r))

定義を満たすもの

このディリクレ過程の定義を満たすような可測な分割の仕方の数というのは非可算無限にある。ここで起こる疑問は「このような定義を満たすようなものはあるのか、あるとすれば条件はなんなのか?」ということである。これを解決するようなapproachがいくつもある。「満たすようなものがあるのか?」について、つまりwell-defineかどうかについてのほうはwikipedia:en:De_Finetti's_theoremを使って示してある。「直接的にGPになるような確率過程を構成してやればいいじゃない」という考えのもので作られているほうは結構たくさんあって、有名なものとしては

  • Gamma Process
  • Chinese Restaurant Process(CRP)
  • Urn Model
  • Stick Breaking Representation

などがある。Chinese Restaurant ProcessはMCMCとStick Breaking Representationは変分ベイズと相性がよいことから推論のときはこの2つがよく使われるようである。

Stick Breaking Representation

http://img.skitch.com/20100216-gnmgutb7hfan64772jg4y4a5a8.jpg

  • G(\theta)は離散分布G(\theta) = \sum_{k=1}^{\infty} \pi_k \delta_{\theta_k}(\theta)と書けることが知られている(詳細は[7])
  • このように書けることを知った上でDPを構成していこうというのがStick Breaking Representationのやり方

Stick Breaking Representationの構成方法*6

  • まず、長さ1の棒をv_1 : (1-v_1)で折って、\pi_1 = v_1とする
  • 次に残った長さ(1-v_1)の棒をv_2 : (1-v_2)で折って、\pi_2 = v_2(1 - v_1)とする
  • これを繰り返すと\pi_k = v_k \prod_{j=1}^{k-1}(1 - v_j)となる
    • v_kはベータ分布v_k \sim \mbox{Beta}(1, \alpha)から生成される
    • 普通、\alphaは超事前ガンマ分布で与えられる(形状、尺度も割と適当でうまくいく)

ベータ分布

http://img.skitch.com/20100216-c1g6rspn8khf14gi8ja9rmcdsk.jpg

  • ベータ分布は定義域が0から1の間で、かなり自在に変形できる分布
    • ただし、単峰性の分布に限定される
  • 無情報変則事前分布としても使わることがある
  • Stick Breaking Representationでは1変数固定

有限で打ち切る

  • 十分大きなkに大しては\theta_kは急激に小さくなっていくため、閾値Kを設けて打ち切って近似しよう、ということが[8]ではやられていた
  • ここではベータ分布の最初の引数のほうを固定していたが、[9]のスライド*7を見ていると以下のような変形があるらしい
    • Dirichlet Process
    • Beta Two-parameter Process
    • Pitman-Yor Process
  • 特にPitman-Yor Processは言語の特性によくfitすることから言語モデルのほうでよく使われるようだ*8
    • この場合は、dに超ベータ分布を与え、補助変数を通して最適化する

Finite Mixtures*9

http://img.skitch.com/20100218-pg8r5uqsy331edscuiejftdxjp.jpg

  • \pi \sim \mbox{Dir}(\alpha / K, \cdots, \alpha / K)
  • c_n \sim \mbox{Multinomial}(\pi)
  • \eta_k \sim G_0
  • y_n | c_n, \eta_1, \cdots, \eta_k \sim F(\cdot|\eta_{c_n})
  • ここでのK\inftyに持って行きたい

Dirichlet Process Mixtures

http://img.skitch.com/20100216-gnmgutb7hfan64772jg4y4a5a8.jpg
DPはG(\theta) = \sum_{k=1}^{\infty} \pi_k \delta_{\theta_k}(\theta)と書けることからも分かるように、離散確率分布であった。それゆえ、連続分布に対して事前分布として使うことができない。これをどうにかしたい。Dirichlet Process Mixturesはこれを解決するためのものである。

Dirichlet Process Mixturesでは、DPから混合モデルのパラメータを取ってくる。

  • G \sim DP(\cdot | G_0, \alpha)
  • \theta_i \sim G(\cdot)
  • x_i \sim p(\cdot | \theta_i)
  • 例えばp(\cdot | \theta)をパラメータ\theta = (\mu, \Sigma)のガウス分布とすれば、これはガウス分布をDirichlet Processによって混合したものになる
    • もちろんp(\cdot | \theta)は任意の密度関数でよい

グラフィカルモデル 1

http://img.skitch.com/20100211-1yrbrpm4erac87d6g629da4ftn.jpg
Dirichlet Process Mixturesでは、DPから混合モデルのパラメータを取ってくる。

  • G \sim DP(\cdot | G_0, \alpha)
  • \theta_i \sim G(\cdot)
  • x_i \sim p(\cdot | \theta_i)
  • サンプルx_iごとにモデルパラメータ\theta_iが違うかもしれない、ということ
  • 現実的には、いくつかのパラメータは同じであることが多い
    • 例 : \theta_1 = \theta_3 = \theta_5 \neq \theta_2 = \theta_4
    • クラスタ効果

グラフィカルモデル 2*10

http://img.skitch.com/20100218-6e1jhj4afppuu4sd35hb6wa7g.jpg

  • 有限混合モデル、DPMとそれと等価なモデルにおけるグラフィカルモデルでの表現の違い
  • 右の3つは出てくるものは同じだけど、中でやっているプロセスが違う
  • SBPは有限混合モデルとほとんど似たようなグラフィカルモデルをしているが、混合要素の数のところが有限ではなく無限まで考えている

Infinite Mixtures 1*11

  • データセット\mathcal{D} = \{\mathbf{x}^{(1)}, \cdots, \mathbf{x}^{(n)}\}に関して、K個の混合要素からなる確率モデルを考える
  • つまりp(\mathbf{x}^{(i)} | \theta) = \sum_{j=1}^K \pi_j p_j(\mathbf{x}^{(i)} | \theta_j) = \sum_{j=1}^K P(s^{(i)} = j | \pi) p_j(\mathbf{x}^{(i)} | \theta_j, s^{(i)} = j)について考える
    • ここで\mathbf{s} = (s^{(1)}, \cdots, s^{(n)})\piがgivenのもとでの*12多項分布P(s^{(1)}, \cdots, s^{(n)} | \pi) = \prod_{j=1}^K \pi_j^{n_j}, \, n_j = \sum_{i=1}^n \delta_{s^{(i)}}(j)である
  • そして、\piの共役事前分布として、Dirichletを持ってくる
    • p(\pi | \alpha) = \frac{\Gamma(\alpha)}{\Gamma(\alpha / K)^K} \prod_{j=1}^K \pi_j^{\alpha / K - 1}
  • 最初はP(s^{(1)}, \cdots, s^{(n)} | \pi)について考えていたが、\piを積分で消して、P(s^{(1)}, \cdots, s^{(n)} | \alpha)について考えると、以下の事後予測分布を得る
    • P(s^{(1)}, \cdots, s^{(1)} | \alpha) = \int d \pi P(\mathbf{s} | \pi) P(\pi | \alpha) = \frac{\Gamma(\alpha)}{\Gamma{(n + \alpha)}} \prod_{j=1}^K \frac{\Gamma(n_j + \alpha / K)}{\Gamma(\alpha / K)}

Infinite Mixtures 2

  • \mathbf{s}_{-i}\mathbf{s}_{-i} = \{s^{(1)}, \cdots, s^{(i-1)}, s^{(i+1)}, \cdots, s^{(n)}\}と定義する
  • このとき条件付き確率P(s^{(i)} = j | \mathbf{s}_{-i}, \alpha)P(s^{(i)} = j | \mathbf{s}_{-i}, \alpha) = \frac{n_{-i, j} + \alpha / K}{n - 1 + \alpha}と書ける
    • ここでn_{-i, j} = \sum_{l \neq i}\delta_{s^{(l)}}(j)
    • つまり、i以外で、クラスタjにいる数
  • クラスタ数Kを無限大に持っていったとき、この条件付き確率P(s^{(i)} = j | \mathbf{s}_{-i}, \alpha)は場合分けが必要
  • jがすでに存在しているクラスタだった場合
    • P(s^{(i)} = j | \mathbf{s}_{-i}, \alpha) = \frac{n_{-i, j}}{n - 1 + \alpha}となる。これは上の式でKを素直に無限大に持って行った場合
  • もう一つのケースは、jが新しくできたクラスタである場合
    • この場合はP(s^{(i)} = j | \mathbf{s}_{-i}, \alpha) = \frac{\alpha}{n - 1 + \alpha}となる、らしいんだけどよく分からない><*13

Dirichlet Process Mixturesの近似

近似するための方法は色々あって[9]では

  • Gibbs sampling [3]
  • Variational approximation [10]
  • Expectation propagation [11]
  • Hierarchical clustering [12]

などが紹介されていた。

Gibbs sampling*14

  • 事後分布p(\mathbf{Z} | \mathbf{X})をサンプリングして近似する方法
  • \mathbf{Z} = (\mathbf{z}_1, \cdots, \mathbf{z}_n)とし、\mathbf{z}_i \sim p(\mathbf{z}_i | \mathbf{Z}_{- i}, \mathbf{X})を繰り返すことでp(\mathbf{Z}|\mathbf{X})を近似する。
    • \mathbf{z}_i \sim p(\mathbf{z}_i | \mathbf{Z}_{- i}, \mathbf{X}) \propto p(\mathbf{X}, \mathbf{Z})はDPMに関しては(割と簡単に)できる
    • 導出は簡単ではないとか

DPMでのGibbs sampling

  • モデルp(\mathbf{X} | \mathbf{Z}, \theta)を指数分布族に決める
    • 例えば正規分布
  • p(\theta)を共役事前分布に、p(\mathbf{Z})をDPMにする
  • p(\mathbf{X}, \mathbf{Z}) = \int d \theta p(\mathbf{X}, \mathbf{Z}, \theta)を解析的に計算する
    • DPMにすることでp(\mathbf{X}, \mathbf{Z})を解析的に計算することができる*15
  • \mathbf{z}_i \sim p(\mathbf{z}_i | \mathbf{Z}_{- i}, \mathbf{X}) \propto p(\mathbf{X}, \mathbf{Z})を繰り返す

アルゴリズム

K := 現在のクラスタ数
for k in 1:K+1 # z_i \in {1, \cdots, K+1} でK+1は新しいクラスタを意味する
p(k) := p(\mathbf{X}, \{\mathbf{Z}_{-i} \cup \mathbf{z}_i\})
end
p(k) = 1となるようにp(k)を正規化する
p(k)からziをサンプリングする

変分ベイズ*16

その他応用例*17

http://img.skitch.com/20100218-tk1q412irfretipkm6xwqn3rm7.jpg
http://img.skitch.com/20100218-fai8nyrbri3rtu1p54kg3ys34h.jpg

  • HDP
  • HMM

参考文献

[1] Ferguson (1973) "A Bayesian Analysis of Some Nonparametric Problems", The Annals of Statistics,Vol.1, No.2, 1973. pdf
[2] Antoniak (1974) "Mixtures of Dirichlet Processes with Applications to Bayesian Nonparametric Problems", The Annals of Statistics,Vol.2, No.6, 1974. pdf
[3] Escobar and West (1995) "Bayesian Density Estimation and Inference using Mixtures", Journal of the American Statistical Association,Vol. 90, No. 430, 1995. pdf
[4] West (1992) "Hyperparameter estimation in Dirichlet process mixture models", ISDS discussion paper 1992-03, 1992. ps
[5] Zhang (2008) "A Very Gentle Note on the Construction of Dirichlet Process", Technical Note. pdf
[6] 栗原 (2008) "Dirichlet Processを用いたクラスタリング", MIRU2008 ワークショップ: 画像処理とその周辺における確率的情報処理の展開. pdf
[7] J. Sethuraman. A constructive definition of Dirichlet priors. Statistica Sinica, 4:693–650, 1994. pdf
[8] 持橋大地, 菊井玄一郎. "無限混合ディリクレ文書モデル", 情報処理学会研究報告 2006-NL-172, pp.47-53, 2006. pdf
[9] Ghahramani (2005) "Non-parametric Bayesian Methods", Advanced Machine Learning Tutorial (Based on UAI 2005 Conference Tutorial), 2005. pdf
[10] Blei, D.M. and Jordan, M.I. (2005) Variational methods for Dirichlet process mixtures. Bayesian Analysis. pdf
[11] Minka, T.P. and Ghahramani, Z. (2003) Expectation propagation for infinite mixtures. NIPS’03 Workshop on Nonparametric Bayesian Methods and Infinite Models. pdf
[12] Heller, K.A. and Ghahramani, Z. (2005) Bayesian Hierarchical Clustering. Twenty Second International Conference on Machine Learning (ICML-2005) pdf
[13] Muller, P. and Quintana, F.A. (2003) Nonparametric Bayesian Data Analysis. pdf
[14] C. E. Rasmussen. The Infinite Gaussian Mixture Model. in Advances in Neural Processing Systems 12, pp. 554-560, MIT Press (2000). pdf
[15] 山本幹雄, 持橋大地. Topicに基づく統計的言語モデルの最前線 ーPLSIからHDPまで― 言語処理学会第12回年次大会チュートリアル資料 pp.11-28, 2006. pdf
[16] pdf
[17] pdf
[18] pdf
[19] pdf
[20] ps
[21] pdf

*1:[9]より

*2:[9]より

*3:[6]より

*4:Chinese Restaurant Process - Beyond 3 Sigmaより

*5:分割は有限

*6:分かりやすかったので、Stick-Breaking Process (SBP) - Beyond 3 Sigmaから図を借りてくる

*7:このスライドは結構分かりやすい

*8:自然言語処理特論最終回 - Seeking for my unique color.

*9:[18]より

*10:[16]より

*11:[9]のスライドを見つつ

*12:\piがパラメータ(隠れ変数的な)。完全にベイズの世界なのでこういう書き方

*13:[16]に書いてあるっぽい

*14:[6]を参考にしつつ

*15:計算過程は追い切れてない。。。

*16:[10]を参考

*17:[17]より