このページで扱う手法の帰属について

ここで解説する方法論は高橋卓也氏（本論文の第一著者）のものです。アクセント類の統合を表す突然変異モデルの設計、それを空間拡散つきのベイズ系統推定の枠組みに統合する定式化、および下記の C++ 実装は、すべて高橋氏の仕事です。小野原彩香は第二著者として共同研究に参加しました。

このページは、その方法論を授業の「7つの観点」のうち観点3（ベイズ）に接続して読むために小野原が整理したものであり、手法そのものの考案者・実装者ではありません。引用・参照の際は下記の原論文（Takahashi, Onohara & Ihara 2024）と Zenodo のコードリポジトリをご参照ください。

1何のための手法か

言語の系統推定で最大の難点は「借用」である。アクセント類の統合は、その難点を回避できる数少ない現象である。

語彙に基づく言語系統推定では、似ている語が共通祖先に由来するのか、それとも接触によって借りられたのかを区別できない。この曖昧さは、推定される樹形を系統ではなく接触の地図に近づけてしまう。

日本語のピッチアクセントには、この難点を回避できる性質がある。アクセント類の統合（accentual class merger）である。古い段階で区別されていた複数のアクセント類が、ある方言ではひとまとまりに合流して区別を失う。いったん失われた区別は、外から単語を借りても復活しない。したがって「どの類とどの類が統合しているか」という情報は、借用によって攪乱されにくい系統の信号になる。

たとえば二モーラ名詞は5つのアクセント類（第1類〜第5類）を持つが、東京方言では第2類と第3類が、第4類と第5類がそれぞれ統合している。これを 1/23/45 と書く。京都方言は 1/23/4/5、伊吹島方言は 1/2/3/4/5（どの類も統合していない）である。従来の方言学は、こうした統合状態の分布から系統関係を論じてきたが、統計モデルとして書かれたことはほとんどなかった。この手法は、統合という現象を確率過程として定式化し、ベイズ系統推定の枠組みに載せる。

関連：系統ネットワーク（Neighbor-Net）　—　距離行列から樹形を仮定せずに描く方法。こちらは樹形そのものを事後分布として推定する。

2観測データ：何を「観測」と見なすか

観測されるのは表層のアクセントパターンと、そこから読み取れる統合状態だけである。

対象は本土方言である。琉球諸語と北海道は除かれている。本土の12地点に位置づけられた15の分類群（現代方言に加え、歴史文書に記録されたアクセント体系を含む）が、系統樹の葉になる。地理空間は日本列島を覆う691個のセル（およそ20km四方、経度 0.25°・緯度 0.167°）に分割される。

アクセントの記述は6つの語類 \(C_1,\dots,C_6\) に分かれる。モーラ数・品詞・活用が異なるものをまとめず、別々に扱うためである。各語類 \(C_i\) が持つアクセント類の数を \(l(i)\) と書くと、\(l(1)=3,\ l(2)=5,\ l(3)=2,\ l(4)=3,\ l(5)=2,\ l(6)=2\) である。

各分類群について観測されるのは2つある。ひとつは表層アクセントパターン \(D_i\)（H/L/M/F/R の並びで記述される、実際に聞こえる型）。もうひとつは表層統合状態 \(M_i\)（1/23/45 のようにスラッシュで書かれる、どの類が合流しているか）。この2つを合わせた \(Y_i = D_i \cap M_i\) が、モデルに与えられるデータである。

3モデル①：なぜ潜在変数を置くのか

観測をそのままマルコフ過程にすると状態数が爆発する。状態数を抑えるための設計が、このモデルの核心にある。

各セルには、語類ごとに2つの潜在変数が割り当てられる。

潜在変数
\[ N_i : \text{拡張統合状態} \qquad b_i = (b_{i1},\dots,b_{il(i)}) : \text{潜在アクセントパターン} \]

\(N_i\) は、表層の統合状態 1/234/5 に「どの類が他を上書きしたか」の情報を足したものである。'1','234','5' と書き、各文字列の先頭の文字が先導クラス（他を上書きした類）を表す。'1','234','5' と '1','324','5' は、表層では同じ 1/234/5 だが、拡張統合状態としては別物である。過去にどの類がどの類を呑み込んだかは観測できないので、これは潜在変数になる。
\(b_i\) は各アクセント類に割り当てられた「本来の」型で、その語類について経験的に記録されている型のいずれかを取る。

表層アクセントパターンは、この2つから決まる。統合されたブロックの内部では、先導クラスの潜在パターンが、そのブロックに属する全ての類の表層パターンとして観測される。二モーラ名詞の例を挙げる。拡張統合状態が \(N_2 = \) '1','23','54'、潜在パターンが \(b_2 = (\mathrm{HHM},\ \mathrm{HLL},\ \mathrm{HHL},\ \mathrm{LHL},\ \mathrm{LLH})\) のとき、表層は

潜在 → 表層
\[ D_2 = (\mathrm{HHM},\ \mathrm{HLL},\ \mathrm{HLL},\ \mathrm{LLH},\ \mathrm{LLH}) \]

ブロック '23' では第2類が先導なので第3類も HLL になり、ブロック '54' では第5類が先導なので第4類も LLH になる。

この設計が捉えている2つの性質

第一に、統合した類は同時に変化する。表層が (LHH, HLL, HLL) から (LHH, LHL, LHL) に変わるのは、2回の突然変異ではなく、先導クラスの潜在パターンが1回変わっただけで説明される。第二に、統合した類では常に同じ型が観測される。この2つは統合という現象の定義そのものであり、二値形質を並べる標準的な系統モデルでは表現できない。

状態数：なぜ 196 で済むのか

二モーラ名詞は5つのアクセント類を持つ。5要素の集合を分割する方法は52通り（ベル数）だが、拡張統合状態はさらに各ブロックの先導クラスを区別するので、分割ごとに「ブロックの大きさの積」だけの場合がある。分割型ごとに数えると、5：\(1\times5=5\)、4+1：\(5\times4=20\)、3+2：\(10\times6=60\)、3+1+1：\(10\times3=30\)、2+2+1：\(15\times4=60\)、2+1+1+1：\(10\times2=20\)、1+1+1+1+1：\(1\times1=1\) で、合計 196通りである。これが \(N_i\) の取りうる値の全てになる。

もし表層アクセントパターン \(D_i\) を直接マルコフ過程にすると、状態数は数千から数万以上になり、枝刈りアルゴリズムの行列べき乗が現実的な時間で終わらない。潜在変数を置いたのは、この計算量の問題を解くためである。

4モデル②：樹形の事前分布を空間から作る

系統樹の事前分布を、人口規模と地理的距離に基づく「どのセルから受け継いだか」の確率から構成する。

離散時間モデルで、各タイムステップの最初に、すべてのセルがどこか1つのセルからピッチアクセント体系全体を受け継ぐ。変数の一部を別のセルから取ってくることは許さない。これにより、進化史が系統ネットワークではなく系統樹として表現されることが保証される。

式 (2.1) — 継承確率
\[ a_{ij} \;=\; K_i^{-1}\,\pi_j \exp\!\left(-\frac{d_{ij}^{\,2}}{2\sigma^{2}}\right) \]
式 (2.2) — 正規化定数
\[ K_i \;=\; \sum_{k=1}^{n} \pi_k \exp\!\left(-\frac{d_{ik}^{\,2}}{2\sigma^{2}}\right) \]

\(a_{ij}\) はセル \(P_i\) がセル \(P_j\) からアクセント体系を受け継ぐ確率。\(\pi_j\) は \(P_j\) の人口規模、\(d_{ij}\) は2セル間の大円距離、\(\sigma\) はガウスカーネルの標準偏差である。\(K_i\) は \(\sum_j a_{ij}=1\) を保証するための定数にすぎない。

この式が言っているのは2つである。近いセルから学びやすい（ガウスカーネル）ことと、人口の多いセルの影響が大きい（\(\pi_j\) の因子）ことである。結果として、地理的に近い分類群がまとまる樹形ほど事前確率が高くなる。人口データは延喜式などの令制国別推定を 1995年の1kmメッシュ人口で空間配分して作られている。

過去 \(\tau\) ステップ分の継承元を \(\tau \times n\) の行列 \(G\) に保存しておくと、\(G\) から分類群の系統樹 \(T\) を切り出せる。\(G\) の事前確率は式 (2.1) の積で計算できるので、セルの空間的相互作用のパターンがそのまま樹形の事前分布になる。

\(\sigma\) は推定されていない

\(\sigma\) を自由パラメータにすると MCMC が収束しなかったため、本分析では \(\sigma = 70\,\mathrm{km}\) に固定されている。これは推定結果ではなく分析上の設定である。異なる \(\sigma\) の値での結果は感度分析として電子補足資料に置かれており、樹形のトポロジーはおおむね変わらなかったと報告されている。

関連：ICAR 空間ランダム効果　—　どちらも「空間的な近さ」を事前分布に入れる発想だが、ICAR が隣接グラフ上の平滑化であるのに対し、こちらは系統の継承元そのものを空間で重みづける。

5モデル③：突然変異を2種類に分ける

統合が起こる過程（\(q\)）と、型が置き換わる過程（\(\alpha\)）を、独立な2つの過程として書く。

受け継いだ後、各セルは突然変異を起こしうる。語類が異なれば突然変異は独立に起こると仮定される（\(N_i,b_i\) の進化は \(i \neq j\) なら \(N_j,b_j\) と独立）。

統合の突然変異 — パラメータ \(q\)
\[ \Pr[\text{1ステップで } N_i \text{ に統合が起こる}] \;=\; \binom{u_i}{2} q \;=\; \tfrac{1}{2}\,u_i\,(u_i-1)\,q \]

\(u_i\) は \(N_i\) に含まれるブロックの数。どのブロックの対も等しい確率 \(q\) で統合すると仮定するので、起こりうる対の数 \(\binom{u_i}{2}\) を掛ける。2つのブロックが統合するとき、それぞれを表す文字列はランダムな順序で連結され、どちらか一方の先導クラスがもう一方を上書きする。たとえば \(N_i=\) '1','23','45' は、'123','45'、'231','45'、'145','23'、'451','23'、'1','2345'、'1','4523' のいずれかに、それぞれ確率 \(q/2\) で変異する。

型の置換 — パラメータ \(\alpha\)
\[ \Pr[\,b_{ij} \to b'\,] = \alpha \quad (b' \neq b_{ij},\ b' \text{ はその語類で記録されている型}) \]

潜在アクセントパターンは要素ごとに変化し、その語類について経験的に記録されている他の型へ等確率で移る。これは Mk モデルの離散時間版である（\(k\) 個の可能な状態のあいだを等しい率で移る）。単モーラ名詞には HL, LH, HM, HH, LL, LM, MH の7型が記録されているので、ある類の潜在パターンが変化する確率は \(6\alpha\) になる。

統合は不可逆で、型の置換は可逆である

この非対称性がモデルの要である。\(q\) の過程はブロックを併合する一方向の過程で、分かれたブロックが元に戻ることはない。\(\alpha\) の過程は型のあいだを行き来する可逆な過程である。系統の信号を担っているのは主に前者で、後者は同じ統合状態を持つ方言のあいだの違いを説明する。

6事後分布と、尤度を計算する仕組み

同時事後分布を書き下し、Felsenstein の枝刈りアルゴリズムで樹の尤度を計算する。

式 (2.3) — 同時事後分布
\[ P(G,q,\alpha,B^{A},N \mid Y) \;\propto\; P(G)\,P(q)\,P(\alpha) \prod_{i=1}^{6} P(b_i^{A})\, P(N_i \mid T,q)\, P(M_i \mid N_i)\, P(D_i \mid T,\alpha,b_i^{A},N_i) \]

\(B^A=(b_1^A,\dots,b_6^A)\) は樹の根における潜在アクセントパターン。根では全ての類が未統合なので、根の潜在パターンは表層パターンと一致する。\(P(q)\), \(P(\alpha)\) はいずれも一様分布。\(P(G)\) は式 (2.1) の \(a_{ij}\) の積。\(P(b_i^A)\) は一様で、たとえば単モーラ名詞なら3つの類に7型から重複なく割り当てるので \(7\times6\times5=210\) 通り、各々 \(1/210\) になる。

\(P(M_i \mid N_i)\) の項は、フィルタとして働く。拡張統合状態 \(N_i\) が表層統合状態 \(M_i\) と整合すれば 1、しなければ 0 である。つまり MCMC は、観測された \(M_i\) と矛盾しない \(N_i\) の空間だけを探索すればよい。

枝刈りアルゴリズム：欠測を「どの状態でもよい」として扱う

最後の項 \(P(D_i \mid T,\alpha,b_i^{A},N_i)\) は、Felsenstein の枝刈りアルゴリズムを不規則な形で実行して計算される。ポイントは、\(N_i\) の下で先導クラスでない類は表層に現れない、つまりその葉では観測が欠測しているという点にある。

葉での初期化
\[ L_{ij}(p,v) = \begin{cases} 1 & \text{第 } j \text{ 類が } v \text{ で先導クラスで、かつ } p = D_{ij} \\ 0 & \text{第 } j \text{ 類が } v \text{ で先導クラスで、かつ } p \neq D_{ij} \\ 1 & \text{第 } j \text{ 類が } v \text{ で先導クラスでない} \end{cases} \]
内部ノードでの再帰（子ノードを \(s_1, s_2\) とする）
\[ L_{ij}(p,v) = \left[\sum_{p'} L_{ij}(p',s_1)\,P(s_1 \text{ が } p' \mid v \text{ が } p)\right] \times \left[\sum_{p'} L_{ij}(p',s_2)\,P(s_2 \text{ が } p' \mid v \text{ が } p)\right] \]
樹全体の尤度
\[ P(D_i \mid T,\alpha,b_i^{A},N_i) = \prod_{j=1}^{l(i)} L_{ij}\!\left(b_{ij}^{A},\ \text{root of } T\right) \]

先導クラスでない葉に 1 を入れるのが、欠測の扱いである。「その葉が何であっても矛盾しない」という意味になり、通常の枝刈りと同じ再帰がそのまま使える。

実装：2状態の遷移確率は閉じた形で書かれている

高橋氏の C++ 実装では、枝刈りの2状態版が行列べき乗ではなく閉じた式で計算されている。遷移率 \(q_{12}, q_{21}\) の離散時間2状態マルコフ連鎖について、\(\lambda = q_{12}/(q_{12}+q_{21})\) と置くと

実装で使われている閉形式（pruning.cpp の prob_two）
\[ P(X_t = 2) \;=\; (1-q_{12}-q_{21})^{t}\,\bigl(P(X_0=2) - \lambda\bigr) + \lambda \]

遷移行列の固有値が \(1\) と \(1-q_{12}-q_{21}\) であることから直ちに出る式である。数値的な行列べき乗を避けることで、MCMC の内側ループが速く、かつ丸め誤差に強くなる。実装では long double が一貫して使われており、葉から根へ向かう動的計画法の配列（\(2M-1\) 個のノード）が明示的に確保されている。

7時間校正と MCMC の設定

樹の長さは本来タイムステップ単位で出てくるので、実時間（年）に対応づける必要がある。1タイムステップ = 25年（ヒトのおよそ1世代）と仮定されている。1ステップ = 10年とした場合の MCMC も実行されており、得られる系統樹はほとんど変わらなかったと報告されている。

分類群（葉）の観測年は、現代方言の記述の調査年・出版年を考慮して西暦1950年とされた。歴史文書に記録されたアクセント体系は、それぞれの成立年代に置かれる。根（共通祖先）の年代の事前分布は西暦450年から1100年の一様分布である。上限を西暦450年としたのは、琉球諸語と本土日本語の分岐時期（近年の研究では3〜7世紀頃とされる）を考慮したためで、琉球方言には本土方言に見られないアクセント類が存在することから、本土のアクセント体系の共通祖先はその分岐より新しいと考えられる。

MCMC の設定値
独立な連鎖の数3
1連鎖あたりの反復数4 × 106
バーンインとして破棄1 × 106
標本抽出の間隔（thinning）103
合計標本サイズ9,000
ガウスカーネルの \(\sigma\)（固定）70 km
空間セル数691

初期樹（\(G\) の初期値）を変えた追加の MCMC が2回実行され、結果が確認されている。MCC 樹の作成には BEAST に同梱の TreeAnnotator、可視化には FigTree が使われている。

関連：ベイズ推定（Stan / brms）　—　事前分布・収束診断・事後分布の読み方という、ベイズ推定に共通する実務はこちらに整理してある。

8何が言えたか

3つのクレード、共通祖先の位置、そして突然変異率の事後分布。

推定された MCC 樹は東京型・京阪型・N型（九州）の3つのクレードから成り、これは従来のピッチアクセント体系の分類と一致する。注目すべきは、列島の西側（広島・大分）と東側（東京・弘前）に分かれて分布する東京型アクセントが、単系統群を形成すると示唆された点である。空間拡散を事前分布に入れてなお、地理的に離れた東京型がまとまった。

最近共通祖先までの年代95% 信用区間
現代方言すべての共通祖先（根）西暦 450 – 825
東京型西暦 1000 – 1550
京阪型（類聚名義抄を含む）西暦 825 – 1100
N型西暦 1050 – 1700

根の年代は古墳時代中期から平安時代初期にあたる。ただし下記の「限界」に記したとおり、この事後分布は事前分布に強く依存している。

共通祖先の「故郷」

行列 \(G\) には過去の各ステップで系統がどのセルにいたかが記録されているので、MRCA の地理的位置の事後分布が得られる。全分類群の共通祖先は現代の京都・大阪・奈良・神戸のあたりにある可能性が高く、京阪型の MRCA も同じ地域に置かれる。N型の MRCA は九州北部または中部に位置する可能性が高い。東京型の MRCA は比較的広く散らばり、近畿周辺に最も高いピーク、東京周辺に2番目のピークがある（東北・九州は可能性が低い）。

ただしこの推定には構造的な偏りがある。式 (2.1) は人口規模の大きいセルからの継承を高く評価するので、MRCA は人口の多いセルに置かれやすい。京都のセルの事後確率が最も高いのは、近畿で最大の人口規模を持つセルであることと、語源に最も近い類聚名義抄のアクセントパターンがこのセルに割り当てられたことの両方が効いている。

突然変異率

統合率 \(q\) の事後中央値は年あたり \(3.37 \times 10^{-4}\)（95%信用区間 \(2.25 \times 10^{-4}\) – \(4.86 \times 10^{-4}\)）で、これはある2つの未統合ブロックの対が統合するのにおよそ 2,967年かかることを意味する（\(1/3.37\times10^{-4}\)）。型の置換率 \(\alpha\) の事後中央値は年あたり \(9.89 \times 10^{-5}\)（95%信用区間 \(6.97 \times 10^{-5}\) – \(1.39 \times 10^{-4}\)）で、単モーラ名詞では7型から他の6型へ移るので率は \(6\alpha\)、すなわち潜在パターンが1回変わるのにおよそ 1,685年かかる。

従来の仮説との食い違い

この結果は、方言学が提唱してきたいくつかの仮説と整合しない。第一に、外輪型（大分・弘前）は金田一によって中輪型（東京・広島）から系統的に離れたものと仮定されてきたが、推定樹は東京型が単系統群を形成することを示唆し、外輪型と中輪型は直接の共通祖先を共有しない側系統ないし多系統のまとまりとして現れた。第二に、中輪型が補忘記に記録された体系から派生したという有力な仮説は支持されなかった。第三に、現代京都方言のアクセント体系は補忘記の体系に由来するのではなく、室町時代から存続していた未記載の系統（統合状態 1/2/3/4/5）に由来すると示唆された。京都方言と補忘記に至る2つの系統は、第2類と第3類の統合をそれぞれ独立に経たという解釈になる。

また、語彙に基づく先行研究（Lee & Hasegawa のベイズ系統解析など）が示す本土方言の東西分裂とも一致しない。これは「本土方言のアクセント体系と語彙が異なる歴史を持つ可能性」を示唆するが、どちらかが誤っているという結論を直ちに導くものではない。

9限界

年代は事前分布に強く依存する

論文自身が繰り返し注意しているとおり、共通祖先の年代の事後分布は事前分布に大きく依存する。上に挙げた「西暦450–825年」という区間は、根の事前分布を西暦450–1100年の一様分布とした前提の下での結果であり、本土アクセントの MRCA が西暦450年より遡らないという仮定に依存している。樹形（トポロジー）が感度分析を通じて安定していたのとは対照的に、年代は結論として弱い。

三モーラ名詞が除外されている

可能な拡張統合状態の数はアクセント類の数とともに急増するため、6〜7類を持つ三モーラ名詞は計算上の理由で分析から除かれている。先行研究のいくつかは三モーラ名詞に基づいて系統を再構成しており、この分析はそこにある系統信号を取りこぼしている可能性がある。ただし、三モーラ名詞のアクセントパターンは多くの方言で二モーラ名詞のそれと独立ではないため、失われた情報がどれだけかは自明ではない。

データが複数の出典の寄せ集めである

日本語方言のアクセントパターンを一貫した基準で網羅したデータベースが公開されていないため、アクセントパターンは複数の異なる研究から収集されている。収録された方言間で記述の基準が揃わないことは避けられず、結果が出典の選択によって偏っている可能性がある。

語類の独立性の仮定は現実には成り立たない

モデルは語類ごとに突然変異が独立に起こると仮定するが、実際のデータでは二モーラ名詞と動詞の第1類のアクセントパターンはほとんどの現代方言で同じである。こうした相関した変化はモデルの下では二重の突然変異を要求するため、樹形に偏りを生じうる。名詞のみ（\(C_1, C_2\)）のサブセットを使った感度分析が行われており、トポロジーはおおむね変わらなかったと報告されている。

潜在変数の解釈は発見的である

「潜在アクセントパターンがマルコフ過程に従う」という設定は、計算を可能にするための工夫であって、言語変化の実体についての主張ではない。論文自身がこの設定を heuristic（発見的）であり解釈がやや難しいと述べている。この手法が推定しているのは、このモデルの下での事後分布であって、言語変化そのものではない。

10出典とコード

方法論と実装は高橋卓也氏によるものです

本ページで解説した突然変異モデルの設計、ベイズ系統推定の枠組みへの統合、および C++ 実装のすべては高橋卓也氏（第一著者）の仕事です。小野原彩香は第二著者として共同研究に参加しました。引用・参照の際は、下記の原論文とコードリポジトリを直接ご参照ください。

実装の構成

MCMC は C++ で書かれており、モデルの各部分がファイルに分かれている。責務の切り分けがモデルの構造にそのまま対応している。

main.cpp — MCMC の本体（サンプリングループ）
parameter_set.cpp — パラメータ集合の提案と受容
genealogy.cpp — 行列 \(G\) の管理。シミュレーションによる初期化、樹の切り出し（cut_tree_exact）、\(G\) の変更提案（mutate_diagram）、式 (2.1) に基づく事前確率（log_prior_from_matrix）
group_integration.cpp — 拡張統合状態 \(N_i\) の表現と統合の突然変異
pruning.cpp — Felsenstein の枝刈り。2状態の遷移確率を閉じた式で計算し、long double で動的計画法を回す
tree.cpp — 樹の構造（子ノードと枝長）
transmission_rate.cpp — 継承確率 \(a_{ij}\) の計算
lexical_class.cpp、location.cpp、id_converter.cpp、csv_input.cpp — 語類・地点・ID 変換・データ読み込み

手法の文献

  • Felsenstein, J. (1981) Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of Molecular Evolution, 17, 368–376.（枝刈りアルゴリズム）
  • Lewis, P. O. (2001) A likelihood approach to estimating phylogeny from discrete morphological character data. Systematic Biology, 50(6), 913–925.（Mk モデル）
  • Bouckaert, R. et al. (2012) Mapping the origins and expansion of the Indo-European language family. Science, 337, 957–960.（言語系統への空間モデルの適用）
  • Lee, S. & Hasegawa, T. (2011) Bayesian phylogenetic analysis supports an agricultural origin of Japonic languages. Proceedings of the Royal Society B, 278, 3662–3669.

この手法がどの観点に属するか

この手法は、樹形 \(T\)・突然変異率 \(q, \alpha\)・根の状態 \(B^A\) のすべてを確率変数として扱い、事前分布と尤度から事後分布を構成する。樹形そのものが事後分布を持つ点が特徴で、「どの系統仮説がどれだけもっともらしいか」を確率として比較できる。これは7つの観点を式で読むの観点3（ベイズ）そのものである。同時に、空間拡散を事前分布に組み込んでいる点で観点3の ICAR と、距離から樹を描く発想の対極にある点で観点1の Neighbor-Net と対比して読める。

このページは小野原が整理したものですが、解説している手法の考案者・実装者は高橋卓也氏です。数値はすべて公開された論文から引用しており、推定値や概算は記載していません。
手法の位置づけの全体像は 7つの観点を式で読む、どの授業回がどの観点に当たるかは 手法と数学の地図 にあります。