Yule過程とSimon過程

2つの古典的な確率過程―Yule過程とSimon過程について、その周辺の話題も含めて解説。MathJaxの習作として書き始めたものが結構な分量になってしまった。

 

目次

導入

1922年にWillisは生物学的属(genus)の数と属が含む種(species)の数がベキ則に従う関係を発見したref。つまり、属が含む種の数を\(n\)、\(n\)の種数を持つ属の数を\(N(n)\)としたとき、広い\(n\)のスケールにわたって \begin{equation} N(n)\propto n^{-\rho},\quad\rho=\mbox{constant} \end{equation} が成り立つ関係である。この事実を説明するために、1925年にYuleが提案した確率モデルがYule過程refである。さらに、属を単語、属が含む種の数を単語の出現頻度と読み替えれば、これは人間の言語活動に見られるZipf則の関係にも類似している。1955年にSimonが提案した単語選択の確率モデル―Simon過程refはこのYule過程をレビューしたもので、この2つの古典的な確率過程は細部が異なるものの基本的なアイデアを共有している。それは意味的には新規性の発生と、それが系の中で勢力を拡大する過程を記述したモデルであり、また数理的には2種類の指数関数的成長からポピュレーションのベキ分布を生じる分岐過程(branching process)のモデルである。本稿ではこのYule過程とSimon過程の両者に共通するアイデアと相違点について、その周辺の話題を含めて解説する。全体の流れがつかみやすいように本文中の数式は結論的なもののみとし、導出の詳細は後半のAppendixにまとめた。

Yule過程

Yule過程は遺伝的変異による属と属に含まれる種の増加を記述しており、一つ一つの変異はそれぞれレート係数一定のポアソン過程に従うと仮定する。つまり、属あるいは種の集団で見た場合には、属数あるいは種数の単位時間あたりの増加数はそれ自身の数に比例する:
  • generic mutation(属が生み出される変異)の単位時間あたりの発生数は属の総数にレート係数\(g\)で比例する。
  • specific mutation(種が生み出される変異)の単位時間あたりの発生数は属が含む種の数にレート係数\(s\)で比例する。
したがって、属の総数および属が含む種の数は時間とともにそれぞれ指数関数的に増加していく。ここでの大きな関心は、十分に長い時間が経過した後の属のサイズ(属が含む種の数)の分布である例えば、属を都市、属のサイズを都市の人口と読み替え、都市の人口動態モデルを考えることも可能だろうrefref。Yuleは、上の2つの仮定に加えて時間を連続変数とみなすことで、サイズ\(n\)の属の確率分布 \begin{equation} p(n)\propto n^{-\rho},\qquad\rho=1+\frac{g}{s} \label{eq:yule_distribution} \end{equation} が導かれることを示した。導出の詳細はAppendix(簡略版厳密版)を参照。つまり、確率分布の関数形は\(g\)と\(s\)の比によって指数が決まるベキ分布の形になる。直感的には\(s\)(属の成長速度)が\(g\)(属の生成速度)に比べて相対的に大きい場合、大きなサイズの属の存在確率が高くなると考えられる。実際に式(\ref{eq:yule_distribution})では\(s\gt g\)のとき分布のテールが長く伸びることが分かる。

修正Yule過程

Yuleは件の論文の中で、個々の属におけるgeneric mutationの発生レートが属ごとに一定ではなく、specific mutationと同様、属のサイズに比例する可能性についても言及している。つまり、多数の種を含む属においては属の変異レートそのものも高まるという考え方である。Yule自身は、属とは生態的に安定な地位を表していて、属の変異はその属に含まれるマイナーな変異種の数には依存しないと考えた。しかし、もし属のサイズにも比例すると考えると、属の増加を記述した式(\ref{eq:mean_growth_of_genera})の右辺は\(\exp(gt)\)ではなく\(\frac{g}{g+s}\exp\bigl((g+s)t\bigr)\)になり、以下、式(\ref{eq:yule_distribution})の導出と全く同様の手順でサイズ\(n\)の属の確率分布 \begin{equation} p(n)\propto n^{-\varrho},\qquad\varrho=2+\frac{g}{s} \label{eq:modified-yule_distribution} \end{equation} が導かれる。これは通常のYule過程と同じく\(g\)と\(s\)の比で指数が決まるベキ分布になっているが、指数に含まれる定数が異なる点に注意。同じ\(g\)と\(s\)の値の組を用いた場合、通常のYule過程に比べて大きなサイズの属の存在確率は小さくなる。属の総数の増加が通常のYule過程に比べて速く、個々の属の勢力が相対的に小さくなるためである。この結果は次に示すSimon過程の結果と等しい。

Simon過程

Simonは文章が書かれる(あるいは読まれる)モデルとして以下のような過程を考えた:
  1. \(N\)個の単語が既に書かれていたとして、\(N+1\)番目に書かれる単語が過去に\(n\)回使われた単語―クラス\([n]\)の単語と呼ぶ―である確率は\(nk(n,N)\)に比例するものとする。ここで\(k(n,N)\)は\(N\)個の単語が書かれた時点でのクラス\([n]\)の語彙のサイズを表す。つまり、\(nk(n,N)\)とはクラス\([n]\)の語彙の総出現数を表す。
  2. 一定確率\(\alpha\)―本稿ではnovelty rateと呼ぶ―で\(N+1\)番目の単語は過去に一度も使われたことのない新しい単語を用いる。
Simonは、この2つの仮定からクラス\([n]\)、つまり出現数\(n\)の単語の確率分布 \begin{equation} p(n)\propto n^{-\gamma},\qquad\gamma=1+\frac{1}{1-\alpha} \label{eq:simon_distribution} \end{equation} が導かれることを示した。導出の詳細はAppendixを参照。Yule過程と同じくベキ分布の形をしており、ベキ指数は\(\alpha\)によって決まる。例えば\(\alpha\rightarrow 1-0\)のとき\(\gamma\rightarrow+\infty\)となるが、これは追加される単語が常に過去に一度も出現したことがない単語である場合に各単語の出現回数は1回、つまり\(p(1)=1\)になることを意味する。一方、\(\alpha=0\)のとき\(\gamma=2\)となるが、これは新たな語彙は生み出されず、初期条件で存在した語彙が行き着く分布を示している。
修正Yule過程との類似を考えたとき、仮定1はspecific mutation(単語の出現数の増加)、仮定2はgeneric mutation(単語の種類の増加)に相当する。両過程のパラメータの関係について整理すると、 \begin{equation} \frac{\alpha}{1-\alpha}=\frac{g}{s},\quad\mbox{or}\quad \alpha=\frac{g}{g+s} \label{eq:yule-simon_relationship} \end{equation} が成り立つ。最初の式の両辺の分子を「種数にかかる」generic mutationのレート、分母を「種数にかかる」specific mutationのレートと考えればよいだろう。つまり、Simon過程における\(\alpha\)とは修正Yule過程におけるgeneric mutationの相対レートに他ならない。式(\ref{eq:yule-simon_relationship})を式(\ref{eq:simon_distribution})に代入してみよう。 \begin{equation} \gamma=1+\frac{1}{1-\frac{g}{g+s}}=1+\frac{g+s}{s}=2+\frac{g}{s}. \end{equation} 式(\ref{eq:modified-yule_distribution})で示した修正Yule過程の指数と同じ結果を得ることが確認できる。

Yule過程とSimon過程の違い

一つの違いは時間発展の進行手順の違いである。Yule過程では属の増加と属が含む種の増加は2つのパラメータ\(g\)と\(s\)で同時進行するのに対し、Simon過程では語彙のサイズの増加と個々の単語の出現数の増加は1つのパラメータ\(\alpha\)で確率的に選択的に進行する。そして、Yule過程では各属において同時並行かつ指数関数的に種数は増加し、全体の種数は実時間に対して指数関数の和で増加していく。Yule過程における全体の種数とはSimon過程における全ての単語の出現数の和、すなわち単語の系列の長さ=時間に対応する概念であるから、Simon過程の枠組みから見たときYule過程では時間が加速していくように見える。しかし、このいわば「時間の密度の違い」は観点の違いであり、重要な違いではない。
外形的に大きな違いは両者のサイズ分布に見られるベキ指数の違いである。この違いは修正Yule過程で成立した式(\ref{eq:yule-simon_relationship})が、通常のYule過程では \begin{equation} \frac{\alpha'}{1-\alpha'}=\frac{gN_g}{sN_s},\quad\mbox{or}\quad\alpha'=\frac{gN_g}{gN_g+sN_s} \label{eq:novelty_rate_in_yule_process} \end{equation} となることに起因する。ここで\(N_g\)は属の総数、\(N_s\)は種の総数を表す。つまり、式(\ref{eq:yule-simon_relationship})の最初の式では右辺の分子の\(g\)と分母の\(s\)がともに種数\(N_s\)にかかっていたのに対し、式(\ref{eq:novelty_rate_in_yule_process})の\(g\)にかかるのは\(N_s\)ではなく\(N_g\)である。この違いの考え方については修正Yule過程で述べた。いま\(g\)と\(s\)を固定したとき、\(N_g\leq N_s\)より\(\alpha\ge\alpha'\)である。つまり、Simonのnovelty rateの観点から見れば、Yule過程に対応する(pseudo) novelty rateはSimon過程よりも小さな値になる。式(\ref{eq:novelty_rate_in_yule_process})の詳しい解はAppendixを参照。
もう一つの重要な違いは、Simon過程ではクラス\([n]\)の語彙が選択されたとしても、その中から具体的にどの単語が選択されるかについては任意性が残されている点であるつまり、Simon過程のみでは個々の単語の成長を計算アルゴリズムとして実装することができない。。Yule過程では種数が等しい属のspecific mutationの発生レートが互いに等しかったことを思い出そう。これをSimon過程の枠組みで捉えるならば、それは出現数が等しい単語の成長レートが互いに等しいことを意味する。そのような単語の選択は、クラス\([n]\)の語彙が次に出現する単語の母体として選ばれたとき、その中から例えば一様乱数で単語を一つ選び出すことで実現できる厳密にはクラス\([n]\)に属するすべての単語が等しく選び出される必要があるが、状態(出現数)は離散的なのでそれは不可能である。。比較のため、Simon過程とnovelty rateの扱いが等価な修正Yule過程について、Simon過程の枠組みを用いて単語の成長過程で起こっていることを眺めてみよう。この枠組みを本稿では「Yule-Simon過程」と呼ぶことにする。いま、単語の系列の長さ\(N\)のときの単語\(i\)の出現数を\(n_i\)とおくと、修正Yule過程では単語\(i\)が次に選び出される確率\(P_i(N)\)は \begin{equation} P_i(N)=(1-\alpha)\cdot\frac{n_ik(n_i,N)}{N}\cdot\frac{n_i}{n_ik(n_i,N)}\quad\biggl[=(1-\alpha)\frac{n_i}{N}\biggr] \label{eq:probability_of_word_selection} \end{equation} と書ける。右辺の第一項はnovelty rateの補確率、つまり次の単語が新しい単語ではなく既存の中から選択される確率(Simonの仮定2)、第二項はクラス\([n_i]\)の語彙が選択される確率(Simonの仮定1)、第三項はクラス\([n_i]\)の中から単語\(i\)が選択される確率である。第三項はSimon過程が要求しない「余分な」項であり、クラスからの等確率な単語の選択を意味する。第二項と第三項を合わせれば、式(\ref{eq:probability_of_word_selection})角括弧のように\(n_i/N\)と書け、これは各単語がそのサイズに比例して選択されることを意味する。このように、Yule過程にはクラスから単語を一つ選び出すルールが暗黙的に導入されている。言い換えれば、Simon過程には具体的な単語の選択に等確率以外の新たなルールを導入する余地が残されており例えば、FILO(First-In, Last-Out。クラスに最後に追加された単語が必ず選ばれる)のようなルールは「より新しいものほど選ばれやすい」という効果を生み出すだろう。、より弱い前提条件でベキ則に従うポピュレーションの確率分布を生成することができるのである。
まとめると、Yule過程とSimon過程の相違点は以下の2点に集約される:
  1. 新規性の発生則に関して、Yule過程はSimon過程よりもnovelty rateが相対的に小さく、その結果、ポピュレーションのベキ指数が小さい(テールが長くなる)。
  2. 新規性の成長則に関して、Yule過程はSimon過程よりも強い仮定を置いている。ポピュレーションのベキ分布を生じるにはSimonの仮定で十分。

Barabási-AlbertグラフとYule-Simon過程のアナロジー

1999年にBarabásiらが提案した成長するグラフモデル(Barabási-Albertグラフ)refとのアナロジーについて考える。Barabásiらは、ノードが次数\(k\)(接続されたエッジの数)に比例してさらなるエッジを獲得していくことで、次数の確率分布が \begin{equation} p(k)\propto k^{-3} \label{eq:solution_of_ba_graph} \end{equation} となることを示した。導出の詳細はAppendixを参照。このモデルはインターネットやウェブの成長メカニズムを議論するために導入され、次数分布に見られるベキ則(スケールフリー性)が大きな関心を集めたことから「スケールフリーグラフ」とも呼ばれた。ノードが次数に比例してエッジを獲得していく過程は"Preferential Attachment"と呼ばれ、しばしば"The Rich Get Richer"といったレトリックで語られる。
いまグラフのノードを属、ノードに接続されたエッジを属に含まれる種と読み替えれば、ノードにエッジが接続されていく過程は属が種数を増やしていく過程とみなすことができる。そのような視点に立てば、Preferential AttachmentとはYule過程における種数増加のルール(属が含む種数に比例してspecific mutationが発生)と等価であることが分かる。さらに、次数の獲得イベントが新規ノードと既存ノードで等しい割合で発生する条件は、Simon過程における単語の生成と成長が等しい割合で発生する条件\(\alpha=0.5\)と等価である。つまり、Barabási-AlbertグラフとはYule-Simon過程の亜種に他ならない。次数分布のベキ指数3は式(\ref{eq:simon_distribution})と\(\alpha=0.5\)から率直に得られる帰結である。Barabásiらの論文の1年後,BornholdtらはSimon過程における単語選択ルールの任意性を再検討することで、Barabási-Albertグラフよりも一般性の高いグラフの生成モデルを提案したref

Yule-Simon過程における単語の成長

ここまでは主に生物学的属や語彙が示すサイズの確率分布について議論してきた。ここではYule-Simon過程(つまりPreferential Attachmentに従うSimon過程)を基礎に、個々の単語がどのように成長するか?つまり単語の個性について議論するこれはCitation Dynamics(学術論文が引用件数をどのように獲得していくか)、あるいはそれに関連してAttention Decay(新規性はどのように忘れられていくか)といった分野にもつながるトピックである。
単語の系列の長さ\(N\)の増加に伴い単語\(i\)はどのように出現数を増やしていくのか?これはBarabási-Albertグラフにおけるノードの成長と同じ問題設定であり、式(\ref{eq:growth_of_degree_in_ba_graph})の導出で用いた議論がそのまま当てはまる。時刻\(N\)のときの単語\(i\)の出現数を\(n_i(N)\)とすると、式(\ref{eq:growth_of_degree_in_ba_graph})にならい \begin{equation} \Delta n_i(N)=n_i(N+\Delta N)-n_i(N)=(1-\alpha)\frac{n_i(N)}{N}\Delta N \end{equation} という単語の出現数の時間発展式が立てられる。これを解くと \begin{equation} n_i(N)=\left(\frac{N}{N_i}\right)^{1-\alpha} \label{eq:growth_of_vocabulary} \end{equation} を得る。\(N_i\)は単語\(i\)が初めて使われた時刻、つまり\(n_i(N_i)=1\)である。つまり、時間に対してベキで成長していく。ただし、この成長関数は\(n_i\)と\(N\)に対して連続量を仮定しており、実際の成長にはゆらぎが存在する。このゆらぎの大きさを見積もってみよう。
単語の成長速度は\(N_i\)が大きいほど遅くなるため、観測時間を固定した場合、\(N_i\)の大きな単語はゆらぎの絶対値も小さくなる。この初出時刻によるゆらぎの大きさの違いを正規化するため、ゆらぎの観測時間のスケールを決める定数\(\lambda(\gt 1)\)と\(N_i\)を用いて、観測時間を\((\lambda-1)N_i\)と定義する。つまり、ゆらぎの観測時刻を\(\lambda N_i\)とし、初出時刻が遅い語彙ほどゆらぎの観測時間を長くする。そして時刻\(\lambda N_i\)において式(\ref{eq:growth_of_vocabulary})が与える期待値からの逸脱を評価する。結論を述べれば、出現数が期待値\(n_i(\lambda N_i)=\lambda^{1-\alpha}\)の\(x\)倍になる確率密度は、\(\alpha\rightarrow 0\), \(\lambda\rightarrow \infty\)の極限で \begin{equation} p(x)\propto \exp(-x) \end{equation} となる。導出の詳細は文献refまたはその要約を参照。つまり、ゆらぎの大きさは指数関数的に減衰し、期待値(\ref{eq:growth_of_vocabulary})からオーダーが外れるようなゆらぎはほとんど存在しないことが分かる。

Heaps則とZipf則との関係

Simon過程が文章生成の確率モデルになることは既に述べた。人間の言語活動において広く観測される統計的傾向にはZipf則以外にもHeaps則と呼ばれる法則がある。Heaps則とは書かれた単語の総数とその語彙のサイズがベキ則に従う関係であり、 \begin{equation} K(N)\propto N^{\beta} \label{eq:heaps_law} \end{equation} と表される。\(N\)は書かれた単語の数、\(K(N)\)は\(N\)個の単語が書かれるまでに出現した単語の種類の数である。一般に\(\beta\)は1より小さい定数になることが知られており、これはSub-Linear則とも呼ばれている。ここではHeaps則とZipf則の関係、そしてSimon過程との関わりについて明らかにする。
まず、Zipf則を \begin{equation} n_r=A r^{-\theta},\qquad A=N/\sum_{r=1}^Kr^{-\theta} \label{eq:zipf_law} \end{equation} と表記する。\(K\)は語彙サイズ、\(r\)は単語を出現頻度の降順で並べた時のランク(ランク1が最頻出単語、ランク\(K\)が最希少単語)、\(n_r\)はランク\(r\)の単語の出現数、\(A\)は規格化定数、\(\theta\)はベキ指数である。Baeza-Yatesらは現象論的な考察によりHeaps則とZipf則の2つのベキ指数は逆数の関係にあることを示したref。考え方は以下のとおり:
  1. Zipf則が成り立つと仮定し、最希少な単語は1回程度しか出現しないと考える。すると、ランク\(K\)の単語の出現数は\(AK^{-\theta}\approx 1\)となり、\(K\)は\(O(A^{1/\theta})\)でスケールすることが分かる。
  2. \(\theta\gt 1\)を仮定すると、式(\ref{eq:zipf_law})の規格化定数\(A\)の分母は有限の値をとり、したがって\(A\)は\(O(N)\)でスケールする。
  3. 以上より、\(K\)は\(O(N^{1/\theta})\)でスケールする。つまり\(\theta=1/\beta\)である。
さて、Simon過程では\(N\)個の単語が使われたときの語彙サイズの期待値は\(K(N)=\alpha N\)であり、このとき\(\beta=1\)のHeaps則(\ref{eq:heaps_law})が成り立つのは自明である。一方、式(\ref{eq:simon_distribution})より、Simon過程から得られるZipf則のベキ指数は\(\theta=1/(\gamma-1)=1-\alpha\)である。つまり、Simon過程ではHeaps則とZipf則の2つのベキ指数は逆数の関係にはなっていない。これはどう理解できるだろうか?この不整合はBaeza-Yatesらの考え方で示した2の仮定(\(\theta\gt 1\))がSimon過程では成立しない(\(\theta=1-\alpha\))ためで、これは\(\alpha\)を定数ではなく時間的に減衰していく形に一般化したSimon過程を考えることで矛盾なく理解できる。時間減衰する\(\alpha(N)\)を用いたSimon過程から得られるZipf則の指数についてはSimon自身が求めており、 \begin{equation} \theta(N)=\frac{K(N)}{\alpha(N)N}\bigl(1-\alpha(N)\bigr) \end{equation} となる。導出の詳細はAppendixを参照。これに式(\ref{eq:heaps_law})と\(\alpha(N)=dK(N)/dN\)を代入すると、 \begin{equation} \theta(N)=\frac{1-\alpha(N)}{\beta},\qquad\alpha(N)\propto N^{\beta-1} \end{equation} が得られる。つまり、Sub-Linear則の下で十分大きな\(N\)については\(\alpha(N)\ll 1\)となり、\(\theta(N)\sim 1/\beta\)という関係が導かれる。機構論的なモデルとしてのSimon過程にHeaps則を仮定することで、Heaps則の指数と逆数の関係にあるZipf則が導かれるという考え方は、仮定されたZipf則からHeaps則の指数を見積もるBaeza-Yatesらの現象論的な考え方とはちょうど逆のロジックになっている。

Appendix

Yule過程の解の導出(簡略版)ref

簡略版では時間に加えて種数を連続変数とみなす。時刻\(t\)における属の総数の期待値\(\langle N_g(t)\rangle\)は \begin{equation} \langle N_g(t)\rangle=\exp(gt) \label{eq:mean_growth_of_genera} \end{equation} となる。適当な時刻\(T\)において年齢が\(\tau\)より小さい属の確率は \begin{equation} P_g(\tau'\lt\tau)=1-\frac{\langle N_g(T-\tau)\rangle}{\langle N_g(T)\rangle}=1-\frac{\exp\bigl(g(T-\tau)\bigr)}{\exp(gT)}=1-\exp(-g\tau) \end{equation} となる(\(T\)に依存しない)。したがって、年齢\(\tau\)の属の確率密度は \begin{equation} p(\tau)=\frac{dP_g(t\lt\tau)}{d\tau}=g\exp(-g\tau) \label{eq:age_distribution_of_genera} \end{equation} となる。一方、年齢\(\tau\)の属に含まれる種数を\(n(\tau)\)とおくと、 \begin{equation} \frac{dn(\tau)}{d\tau}=sn(\tau),\qquad n(0)=1 \end{equation} より \begin{equation} n(\tau)=\exp(s\tau) \label{eq:number_of_species_in_age-tau-genera} \end{equation} を得る。\(\tau\)と\(n\)が互いに一価関数の関係にあると仮定すると、式(\ref{eq:number_of_species_in_age-tau-genera})を変形して \begin{equation} \tau=\frac{1}{s}\ln(n) \label{eq:tau-ns_relationship} \end{equation} が得られ、また \begin{equation} p(\tau)d\tau=p(n)dn \label{eq:probability_of_age-tau_genera} \end{equation} が成り立つ。式(\ref{eq:tau-ns_relationship})を式(\ref{eq:age_distribution_of_genera})に代入すると \begin{equation} p(\tau)=g\exp\left(-g\cdot\frac{1}{s}\ln(n)\right)=gn^{-g/s}. \label{eq:age_distribution_of_genera2} \end{equation} また式(\ref{eq:tau-ns_relationship})より \begin{equation} \frac{d\tau}{dn}=\frac{1}{sn}. \label{eq:tau-ns_derivative} \end{equation} 式(\ref{eq:probability_of_age-tau_genera})(\ref{eq:age_distribution_of_genera2})(\ref{eq:tau-ns_derivative})より \begin{equation} p(n)=p(\tau)\frac{d\tau}{dn}=\frac{gn^{-g/s}}{sn}=\frac{g}{s}n^{-\left(1+\frac{g}{s}\right)}. \end{equation} したがって \begin{equation} p(n)=(\rho-1)n^{-\rho},\qquad\rho=1+\frac{g}{s} \end{equation} となり、式(\ref{eq:yule_distribution})が導かれる。

Yule過程の解の導出(厳密版)refref

厳密版では簡略版で追加した種数の連続性などの仮定を用いない導出を行う。年齢\(\tau\)の属が種数1である確率、つまり一度もspecific mutationが発生しない確率を\(q_1(\tau)\)とおくと、ポアソン分布より \begin{equation} q_1(\tau)=\frac{(s\tau)^0\exp(-s\tau)}{0!}=\exp(-s\tau) \end{equation} となる。同様に、年齢\(\tau\)の属が種数2である確率を\(q_2(\tau)\)とおくと、 \begin{align} q_2(\tau)&=\int_0^{\tau}d\tau_1\left[\exp(-s\tau_1)s\exp\bigl(-s(\tau-\tau_1)\bigr)\exp\bigl(-s(\tau-\tau_1)\bigr)\right]\nonumber\\ &=\exp(-s\tau)\bigl(1-\exp(-s\tau)\bigr) \end{align} となる。これは年齢\(0\)から\(\tau\)のどこかの点\(\tau_1\)で1度だけspecific mutationが発生する確率を意味する。つまり
  1. \(\tau_1\)まで変異が発生しない確率:\(\exp(-s\tau_1)\)
  2. \(\tau_1\)で2つに増えた種がそれぞれ\(\tau\)まで変異が発生しない確率:\(\exp\bigl(-s(\tau-\tau_1)\bigr)\)
  3. 微小区間で変異が1回発生する確率:\(sd\tau_1\)
の積について、可能な\(\tau_1\)の範囲で積分した値である。同様の考え方で種数\(n\)についても \begin{equation} q_n(\tau)=\exp(-s\tau)\bigl(1-\exp(-s\tau)\bigr)^{n-1} \label{eq:q_n} \end{equation} と求まる。年齢\(\tau\)の属の存在確率が式(\ref{eq:age_distribution_of_genera})で与えられることを使えば、種数\(n\)を持つ属の存在確率は式(\ref{eq:q_n})を式(\ref{eq:age_distribution_of_genera})で重み付け積分することで \begin{align} p(n)&=\int_0^{\infty}p(\tau)q_n(\tau)d\tau\nonumber\\ &=\int_0^{\infty}g\exp\bigl(-(g+s)\tau\bigr)\bigl(1-\exp(-s\tau)\bigr)^{n-1}d\tau \label{eq:p_n} \end{align} と求まる。この式を展開していくと階乗が現れ、Stirlingの近似を用いることでベキ関数の式(\ref{eq:yule_distribution})が導かれる。Simkinらの簡明な方法とBacaërの丁寧な方法の2つを紹介する。

Simkinらの方法ref

変数変換 \(x=\exp(-s\tau)\) によって積分範囲が0~1に変わり、\(dx/d\tau=-sx\)より \begin{align} p(n)&=\int_0^1gx^{\frac{g+s}{s}}(1-x)^{n-1}\frac{dx}{sx}\nonumber\\ &=\frac{g}{s}\int_0^1x^{\frac{g}{s}}(1-x)^{n-1}dx \end{align} となる。これはそのままベータ関数の形をしているので \begin{align} p(n)&=\frac{g}{s}\ \mathrm{B}\left(1+\frac{g}{s}, n\right)\nonumber\\ &=\frac{g}{s}\ \frac{\Gamma(n)\ \Gamma\left(1+\frac{g}{s}\right)}{\Gamma\left(n+1+\frac{g}{s}\right)} \end{align} を得る。以下、式(\ref{eq:stirling})以降を参照。

Bacaërの方法ref

部分積分によって \begin{align} p(n)&=\left[-\frac{g}{g+s}\exp\bigl(-(g+s)\tau\bigr)\bigl(1-\exp(-s\tau)\bigr)^{n-1}\right]_0^{\infty}\quad\mbox{(This term disappears.)}\nonumber\\ &\quad -\int_0^{\infty}d\tau\left[-\frac{g}{g+s}\exp\bigl(-(g+s)\tau\bigr)(n-1)\bigl(1-\exp(-s\tau)\bigr)^{n-2}s\exp(-s\tau)\right]\nonumber\\ &=\frac{(n-1)s}{g+s}\int_0^{\infty}d\tau\Bigl[g\exp\bigl(-(g+2s)\tau\bigr)\bigl(1-\exp(-s\tau)\bigr)^{n-2}\Bigr] \label{eq:p_n_exp} \end{align} を得る。また、式(\ref{eq:p_n})を \begin{align} p(n)&=\int_0^{\infty}d\tau\Bigl[g\exp\bigl(-(g+s)\tau\bigr)\bigl(1-\exp(-s\tau)\bigr)^{n-2}\bigl(1-\exp(-s\tau)\bigr)\Bigr]\nonumber\\ &=\int_0^{\infty}d\tau\Bigl[g\exp\bigl(-(g+s)\tau\bigr)\bigl(1-\exp(-s\tau)\bigr)^{n-2}\Bigr]-\int_0^{\infty}d\tau\Bigl[g\exp\bigl(-(g+2s)\tau\bigr)\bigl(1-\exp(-s\tau)\bigr)^{n-2}\Bigr] \end{align} と変形すると、右辺第1項はそのまま\(p(n-1)\)、右辺第2項は式(\ref{eq:p_n_exp})の積分部分と等しいことが分かる。これを利用すると、 \begin{equation} p(n)=p(n-1)-p(n)\frac{g+s}{(n-1)s} \end{equation} となり、 \begin{equation} p(n)=\frac{(n-1)s}{g+ns}p(n-1) \end{equation} を得る。\(p(1)=g/(g+s)\)を用いれば、 \begin{equation} p(n)=\frac{g}{g+s}\cdot\frac{s}{g+2s}\cdot\frac{2s}{g+3s}\cdot\frac{3s}{g+4s}\cdots\frac{(n-1)s}{g+ns} \label{eq:yule-solution} \end{equation} を得る。この式は、\(n\)が大きくなるにつれて乗算される係数が1に近づいていくため、確率の減衰が遅くなる、つまり確率分布のテールが伸びることを示唆している。
\(u=g/s\)とおき、式(\ref{eq:yule-solution})をさらに変形すると、 \begin{align} p(n)&=\frac{u}{1+u}\cdot\frac{1}{2+u}\cdot\frac{2}{3+u}\cdot\frac{3}{4+u}\cdots\frac{n-1}{n+u}\nonumber\\[6pt] &=\frac{u\ \Gamma(n)\ \Gamma(1+u)}{\Gamma(n+1+u)} \label{eq:yule:bacaer-solution-Gamma} \end{align} となり、Simkinらの方法と同様の結果を得る。最後に大きな\(n\)におけるStirlingの近似 \begin{equation} \ln\Gamma(n)\simeq n\ln(n)-n+\frac{1}{2}\ln(2\pi n) \label{eq:stirling} \end{equation} を用いれば \begin{equation} \ln\bigl(p(n)\bigr)\simeq-(1+u)\ln(n)+\mbox{constant}. \end{equation} \(u\)を戻すと \begin{equation} p(n)\propto n^{-\rho},\qquad\rho=1+\frac{g}{s} \end{equation} を得る。

Simon過程の解の導出refref

ここでは\(k(n,N)\)を\(N\)個の単語が書かれた時点におけるクラス\([n]\)の語彙サイズの期待値とみなす。すると語彙サイズに関する時間発展式は \begin{align} k(n,N+1)-k(n,N)&=(1-\alpha)\left\{\frac{(n-1)k(n-1,N)}{N}-\frac{nk(n,N)}{N}\right\}\qquad (n\geq 2), \label{eq:simon:rate_equation_of_class-n}\\ k(1,N+1)-k(1,N)&=\alpha-(1-\alpha)\frac{k(1,N)}{N} \label{eq:simon:rate_equation_of_class-1} \end{align} というマスター方程式の形で書ける。両式の左辺はそれぞれ\(N+1\)番目の単語が書かれたときのクラス\([n]\)およびクラス\([1]\)の語彙の増加量の期待値を表す。また、\(nk(n,N)\)はクラス\([n]\)の語彙の総出現数であるから、既存語彙の中からクラス\([n]\)の語彙が選択される相対確率は\(nk(n,N)/N\)である。したがって、式(\ref{eq:simon:rate_equation_of_class-n})の右辺第一項はクラス\([n-1]\)の語彙が選択されクラス\([n]\)に移動することによるクラス\([n]\)の増加量、第二項はクラス\([n]\)の語彙が選択されクラス\([n+1]\)に移動することによるクラス\([n]\)の減少量を表す。式(\ref{eq:simon:rate_equation_of_class-1})の右辺第一項はクラス\([1]\)の増加が新しい語彙の生成(確率\(\alpha\))によってのみ生じることを意味する。
いま、すべての語彙の中でクラス\([n]\)の語彙の割合を\(p(n,N)\)とおくと、 \begin{equation} p(n,N)=\frac{k(n,N)}{\alpha N} \label{eq:simon:probability_of_class-n} \end{equation} と書ける。右辺の分母は語彙のサイズである。系が定常状態、つまり新しい単語が書かれてもクラスに属する語彙の割合が変わらない状態に達していると仮定すると、 \begin{equation} p(n)\equiv p(n,N+1)=p(n,N) \label{eq:simon:stationary_state} \end{equation} が成り立つ。式(\ref{eq:simon:probability_of_class-n})(\ref{eq:simon:stationary_state})を式(\ref{eq:simon:rate_equation_of_class-n})(\ref{eq:simon:rate_equation_of_class-1})に代入することで \begin{align} p(n)&=\frac{(1-\alpha)(n-1)}{1+(1-\alpha)n}p(n-1)\qquad (n\geq 2),\\ p(1)&=\frac{1}{2-\alpha} \end{align} を得る。簡単のため\(u=1/(1-\alpha)\)とおくと、 \begin{align} p(n)&=\frac{1}{1+(1-\alpha)}\cdot\frac{1-\alpha}{1+2(1-\alpha)}\cdot\frac{2(1-\alpha)}{1+3(1-\alpha)}\cdot\frac{3(1-\alpha)}{1+4(1-\alpha)}\cdots\frac{(1-\alpha)(n-1)}{1+n(1-\alpha)}\\[6pt] &=\frac{u}{1+u}\cdot\frac{1}{2+u}\cdot\frac{2}{3+u}\cdot\frac{3}{4+u}\cdots\frac{n-1}{n+u} \label{eq:simon:solution}\\[6pt] &=\frac{u\ \Gamma(n)\Gamma(1+u)}{\Gamma(n+1+u)} \label{eq:simon:solution-gamma}\\ \end{align} となり、Yule過程の解の導出で出現した式(\ref{eq:yule:bacaer-solution-Gamma})と同じ形の解を得る。\(u\)を戻して大きな\(n\)の漸近近似を用いることで \begin{equation} p(n)\propto n^{-\gamma},\qquad\gamma=1+\frac{1}{1-\alpha} \label{eq:simon:solution-approx} \end{equation} となり、式(\ref{eq:simon_distribution})が導かれる。

Yule過程におけるnovelty rateの評価ref

Simon過程における\(\alpha\)がYule過程の2つのmutation rateを用いてどのように表されるかについてもう少し細かく見てみると、 \begin{equation} \alpha'\simeq\left\{ \begin{array}{ll} 1-\frac{s}{g}&\qquad (g\gt s),\\ \frac{g(s-g)}{s^2}\exp\bigl(-(s-g)t\bigr)&\qquad (g\lt s) \label{eq:solution_of_simons_alpha_in_yule_process} \end{array} \right. \end{equation} という関係が導かれる。\(g\gt s\)の場合に\(\alpha'\)は定数、\(g\lt s\)の場合には時間とともに0に漸近していく。導出は以下:
Yule過程において、時刻\(t\)での種数の期待値\(\langle N_s(t)\rangle\)は \begin{equation} \langle N_s(t)\rangle=\int_0^tg\langle N_g(t_1)\rangle\exp\bigl(s(t-t_1)\bigr)dt_1+\exp(st) \end{equation} と表せる。右辺第1項の積分について、\(g\langle N_g(t_1)\rangle\)は時刻\(t_1\)に新たに生成される属数、\(\exp\bigl(s(t-t_1)\bigr)\)はその新しい属が時刻\(t\)までに生み出す種数を表す。したがって、\(t_1\)を可能な範囲で積分した値が、時刻\(t\)までに生成された属が生み出す種数の総数になる。右辺第2項は最初の属が時刻\(t\)までに生み出した種数である。これを部分積分で解くと、 \begin{align} \langle N_s(t)\rangle&=\frac{g}{g-s}\bigl(\exp(gt)-\exp(st)\bigr)+\exp(st)\nonumber\\[6pt] &=\frac{g}{g-s}\exp(gt)+\frac{s}{s-g}\exp(st) \end{align} を得る。大きな\(t\)の極限を考えたとき、\(g>s\)ならば\(\exp(st)\)の項が落ち、\(s>g\)ならば\(\exp(gt)\)の項が落ちる。したがって \begin{equation} \langle N_s(t)\rangle\simeq\left\{ \begin{array}{ll} \frac{g}{g-s}\exp(gt) &\quad (g\gt s),\\ \frac{s}{s-g}\exp(st) &\quad (g\lt s) \end{array}\right. \label{eq:number_of_species_in_yule_process} \end{equation} となり、これを式(\ref{eq:novelty_rate_in_yule_process})に代入することで式(\ref{eq:solution_of_simons_alpha_in_yule_process})を得る。

Barabási-Albertグラフの解の導出ref

Barabási-Albertグラフとは、単位時間ごとにネットワークにノードを1つ追加し、新規ノードから既存のノードにその次数に比例した確率でエッジを張るネットワークの生成モデルである。本来は時間も次数も離散量だが、ここでは連続変数とみなして解を導出する。
新規ノードと既存のノードの間に張るエッジの数を\(m\)、時刻\(t\)でのノード\(i\)の次数を\(k_i(t)\)とすると、\(\Delta t\)経過後の次数の増分\(\Delta k_i(t)\)は \begin{equation} \Delta k_i(t)=\frac{mk_i(t)}{\sum_i k_i(t)}\Delta t \end{equation} となる。右辺が次数に比例したエッジの獲得を表し、Simonの枠組みから見た修正Yule過程の式(\ref{eq:probability_of_word_selection})の\(n_i/N\)に対応している。\(\sum_ik_i(t)=2mt\)という関係を用いて、\(\Delta k_i\rightarrow 0\), \(\Delta t\rightarrow 0\)の極限を考えれば、 \begin{equation} \frac{dk_i}{dt}=\frac{k_i}{2t} \end{equation} が得られる。積分し、 \begin{equation} \int\frac{1}{k_i}dk_i=\int\frac{1}{2t}dt. \end{equation} これを解いて \begin{equation} k_i(t)=At^{1/2}. \label{eq:growth_of_degree_in_ba_graph_} \end{equation} \(A\)は積分定数。ノード\(i\)がネットワークに追加された時刻\(t_i\)を式(\ref{eq:growth_of_degree_in_ba_graph_})に代入すれば、\(A{t_i}^{1/2}=m\)より、ノード\(i\)についての次数の成長関数 \begin{equation} k_i(t)=m\left(\frac{t}{t_i}\right)^{1/2} \label{eq:growth_of_degree_in_ba_graph} \end{equation} が求まる。ある次数\(k\)に対して\(k_i(t)\lt k\)となる確率を考えると、式(\ref{eq:growth_of_degree_in_ba_graph})より \begin{equation} P\bigl(k_i(t)\lt k\bigr)=P\biggl(t_i\gt\left(\frac{m}{k}\right)^2t\biggr) \end{equation} が得られる。時間軸\(t\)の直線上にノードの参加時刻は等間隔で並んでいるので、\((0,t)\)の範囲で\(t_i\)が時刻\(t_0\)よりも大きい確率は\(1-t_0/t\)となる。したがって、 \begin{equation} P\bigl(k_i(t)\lt k\bigr)=1-\left(\frac{m}{k}\right)^2. \end{equation} 次数\(k\)の確率密度\(p(k)\)は \begin{equation} p(k)=\frac{dP\bigl(k_i(t)\lt k\bigr)}{dk}=2m^2k^{-3} \end{equation} となり、式(\ref{eq:solution_of_ba_graph})を得る。

novelty rateが時間減衰するSimon過程の解の導出ref

新しい単語の生成確率\(\alpha\)が定数ではなく、時間的に減衰する場合のSimon過程について考える。まず、式(\ref{eq:simon:rate_equation_of_class-n})(\ref{eq:simon:rate_equation_of_class-1})の\(\alpha\)を\(\alpha(N)\)に置き換える。また、式(\ref{eq:simon:probability_of_class-n})の右辺分母(\(N\)個のワードが書かれたときの語彙サイズ)はもはや\(\alpha N\)とは書けないため、\(K(N)\)と表記する。つまり \begin{equation} p(n,N)=\frac{k(n,N)}{K(N)},\qquad dK(N)/dN=\alpha(N). \label{eq:simon:probability_of_class-n_time-decaying-alpha} \end{equation} 定常状態(\ref{eq:simon:stationary_state})を仮定し、式(\ref{eq:simon:probability_of_class-n_time-decaying-alpha})を式(\ref{eq:simon:rate_equation_of_class-n})(\ref{eq:simon:rate_equation_of_class-1})に代入することで \begin{align} p(n)&=\frac{\bigl(1-\alpha(N)\bigr)(n-1)}{\frac{\alpha(N)N}{K(N)}+\bigl(1-\alpha(N)\bigr)n}p(n-1)\qquad (n\geq 2),\\ p(1)&=\frac{1}{1+\frac{K(N)}{\alpha(N)N}\bigl(1-\alpha(N)\bigr)} \end{align} を得る。簡単のため\(u=\frac{\alpha(N)N}{K(N)\bigl(1-\alpha(N)\bigr)}\)とおくと、式(\ref{eq:simon:solution-gamma})と全く同じ形になり、 \begin{equation} p(n)\propto n^{-\gamma},\qquad\gamma=1+\frac{\alpha(N)N}{K(N)\bigl(1-\alpha(N)\bigr)} \label{eq:simon:solution-approx_time-decaying-alpha} \end{equation} を得る。もちろん、\(\alpha\)が定数であれば\(\alpha(N)N=K(N)\)となるので、式(\ref{eq:simon:solution-approx_time-decaying-alpha})は式(\ref{eq:simon:solution-approx})と同一となる。

参考文献

 
Willis, J. Christopher: Age and area, a study in geographical distribution and origin of species, Cambridge University Press, 1922.
Yule, G. Udny: A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis, F.R.S., Phil. Trans. Roy. Soc. B, 213, pp.21-87, 1925.
Bacaër, Nicolas: Yule and evolution (1924), A short history of mathematical population dynamics, pp.81-88, Springer & Business Media, 2011.
Simkin, Mikhail V. & Roychowdhury, Vwani P.: Re-inventing Willis, Physics Reports, 502, pp.1-35, 2011.
Simon, Herbert A.: On a class of skew distribution functions, Biometrika, pp.425-440, 1955.
Barabási, Albert-László & others: Mean-field theory for scale-free random networks, Physica A, 272, 1, pp.173-187, 1999.
Bornholdt, Stefan & Ebel, Holger: World Wide Web scaling exponent from Simon’s 1955 model, Phys. Rev. E, 64, 3, 035104, 2001.
Baeza-Yates, Ricardo & Navarro, Gonzalo: Block addressing indices for approximate text retrieval, Journal of the American Society for Information Science, 51, 1, pp.69-82, 2000.
Gabaix, Xavier: Zipf's Law for Cities: An Explanation, The Quarterly Journal of Economics, 114, 3, pp.739-767, 1999.
Hashimoto, Yasuhiro: Growth fluctuation in preferential attachment dynamics, Phys. Rev. E, 93, 4, 042130, 2016.
Tria, Francesca and others: The dynamics of correlated novelties, Scientific Reports, 4, 5890, 2014.
Kauffman, Stuart A.: Investigations, Oxford University Press, 2000.
例えばFrank, S. A.: George Price’s contributions to Evolutionary Genetics: J. Theor Biol, 175, 3, pp.373-388, 1995.

脚注

 

更新履歴

20170605 v2.1.
20170106 v2.
20150401 v1.