働きアリの法則と微分方程式 - 現実と数学の区別が付かない

こんな面白い記事がはてブに上がっていました。

ざっくり言うと，350回以上も引用された，「ロサダの法則」というものを提唱した心理学の論文があるのですが，そこで使われている微分方程式のモデルに妥当性が全く無いことが明らかになったというお話です。その事を指摘した論文では，微分方程式のモデルが妥当であるために満たすべき条件が5つ (VA1)～(VA5) 挙げられており，それらは一つでも満たされなかったらモデルに妥当性が無いことになるのに，「ロサダの法則」の論文はそのいずれも満たしていない（VA1のほんの一部を除いて）とされています。

http://www.physics.nyu.edu/sokal/complex_dynamics_final_clean.pdf

In the preceding section we listed five criteria, VA1–VA5, all of which must be met if an application of differential equations is to be valid. Here we shall demonstrate that none of these criteria have been met by Losada’s work, except arguably a small part of VA1.

社会科学の分野における現象を微分方程式のモデルで解析する場合は，特にそのモデルの妥当性に注意を払わなければなりません。使い方を誤るとリーマンショックのような痛い目にあってしまいます。まあ金融工学よく知らんけど（免責）。

例えばあるサービスの利用者が年々指数的に増加していた場合，利用者数を $X(t)$ ， $a>0$ を定数として，

\begin{align} \cfrac{dX}{dt}=aX \end{align}

のような連続モデルを考えるでしょう。これは利用者の人数に比例して利用者が増加するというモデルです。しかし利用者が自己分裂で増殖していない限り，なぜそう言えるのか，その根拠を明らかにしなければなりません。単純にたくさんの人が利用していると，他の人も利用し出すと考えていいのか？他の原因（たとえばサービスの質の向上や価格の低下など）はないのか？常に慎重にあるべきです。

これも「現実と数学の区別が付かない」案件でもあるので，拾い食いせずにはいられません。それっぽい微分方程式のモデルをでっちあげるという遊びを私もしてみたくなりました。

働きアリの法則

テーマに選んだのは働きアリの法則と呼ばれるものです。

働きアリの法則 - Wikipedia

働きアリの中で，よく働いているアリと，普通に働いているアリと，ずっとサボっているアリの割合は，2:6:2になるという法則です。よく働いているアリだけを集めても，サボっているアリだけを集めても，最終的には2:6:2になるというのです。

この働きアリたちの労働意欲を表した微分方程式のモデルを考えてみましょう。参考にするのはSIRモデルです。

さらにこのブログの「追記：数値計算」にあるRのコードをまるパクリして数値実験もしてみることにします。

時刻 $t$ において，

$X(t)$ を働き者のアリの数
$Y(t)$ を普通のアリの数
$Z(t)$ をサボっているアリの数

とします。ここで以下のような仮定を置きます。

アリは生活の中で出会うアリの労働意欲に影響されて，労働意欲が上下する
サボっているアリに出会うと，ある一定の確率で労働意欲が上がる
働き者のアリに出会うと，ある一定の確率で労働意欲が下がる
普通のアリが普通のアリに出会うと，ある一定の確率で気まぐれに，働き者になったり，サボるアリになったりする
働き者のアリとサボっているアリは普通のアリに出会っても影響を受けない

周りがサボっていたら危機感を感じて働き出し，周りが働いていたら安心してサボってしまうというモデルです。本来ならばこの仮定が本当に正しいのか慎重になるべきですが，今回はでっち上げなので別にどうでもいいのです。なんかそれっぽい感じもするし。

この仮定に従って微分方程式を立ててみましょう。アリたちの数の時間的な変化の割合 $\cfrac{dX}{dt}, \cfrac{dY}{dt}, \cfrac{dZ}{dt}$ はどうなるでしょうか。
\begin{align}
\cfrac{dX}{dt}&=&???\\
\cfrac{dY}{dt}&=&???\\
\cfrac{dZ}{dt}&=&???
\end{align}

まず，働き者のアリ同士が出会った時のことを考えてみます。働き者のアリ同士が出会う確率は $X^2$ に比例し，働き者のアリは働き者のアリと出会ったことで，一定の確率で労働意欲が下がり，普通のアリになってしまいます。そこで，「一定の確率」から決まるある定数 $a>0$ を用意して， $\cfrac{dX}{dt}=$ の右辺に $-aX^2$ ， $\cfrac{dY}{dt}=$ の右辺に $aX^2$ を加えます。

\begin{align}
\cfrac{dX}{dt}&=&-aX^2\\
\cfrac{dY}{dt}&=&aX^2\\
\cfrac{dZ}{dt}&=&
\end{align}

次に働き者のアリと普通のアリが出会った時のことを考えてみます。働き者のアリと普通のアリが出会う確率は $XY$ に比例し，普通のアリは働き者のアリと出会ったことで，一定の確率で労働意欲が下がり，サボるアリになってしまいます。よってある定数 $b>0$ を用意して， $\cfrac{dY}{dt}=$ の右辺に $-bXY$ ， $\cfrac{dZ}{dt}=$ の右辺に $bXY$ を加えます。

\begin{align}
\cfrac{dX}{dt}&=&-aX^2\\
\cfrac{dY}{dt}&=&aX^2-bXY\\
\cfrac{dZ}{dt}&=&bXY
\end{align}

以下，同様な感じで計算していくと，次の微分方程式が得られます。

$a,b,c,d,e>0$ は定数。
\begin{align}
\cfrac{dX}{dt}&=&-aX^2+cYZ+eY^2\\
\cfrac{dY}{dt}&=&aX^2-bXY-cYZ+dZ^2-2eY^2\\
\cfrac{dZ}{dt}&=&bXY-dZ^2+eY^2
\end{align}

定数を適当な値に固定して，いろいろな初期値に対して数値計算をしてみましょう。 $a=0.1, b=0.028, c=0.028,d=0.1,e=0.002$ ， $X+Y+Z=10$ として，先ほどのSIRモデルのブログを参考にRで数値計算をしてみます。

library(deSolve)

parameters <- c(a=0.1, b=0.028, c=0.028,d=0.1,e=0.002)
state <- c(X=0, Y=9, Z=1)

Ants <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    dX <- -a*X*X+c*Y*Z+e*Y*Y
    dY <- a*X*X-b*X*Y-c*Y*Z+d*Z*Z-2*e*Y*Y
    dZ <- b*X*Y-d*Z*Z+e*Y*Y
    list(c(dX, dY, dZ))
  })
}

times <- seq(0, 30, by=0.1)

ode_out <- ode(y=state, times=times, func=Ants, parms=parameters)

###ggplot2###
library(tidyr)
library(dplyr)
library(ggplot2)
library(scales)
ode_out2 <-as.data.frame(ode_out) %>% 
  gather(variable, value, -time)

theme_set(theme_bw(20))

ggplot(ode_out2)+ 
  geom_line(aes(x=time,y=value,colour=variable),size=1.3)+
  scale_colour_manual(values = c("royalblue","orange2","forestgreen" ),name="")+
  scale_y_continuous(breaks=seq(0,10,1))

初期値が $(X,Y,Z)=(0,9,1)$ のときはこうなります。

あっ！ $X:Y:Z$ が最終的に $2:6:2$ になってますね！

他の初期値に対しても計算してみましょう

###初期値の変更###
state <- c(X=0, Y=0, Z=10)
###再計算###
ode_out <- ode(y=state, times=times, func=Ants, parms=parameters)
###描画###
ode_out2 <-as.data.frame(ode_out) %>% 
  gather(variable, value, -time)

theme_set(theme_bw(20))

ggplot(ode_out2)+ 
  geom_line(aes(x=time,y=value,colour=variable),size=1.3)+
  scale_colour_manual(values = c("royalblue","orange2","forestgreen" ),name="")+
  scale_y_continuous(breaks=seq(0,10,1))

初期値が $(X,Y,Z)=(0,0,10)$ のとき。

初期値が $(X,Y,Z)=(7,3,0)$ のとき。

いずれの場合も $X:Y:Z$ が最終的に $2:6:2$ になってます！働きアリの法則の微分方程式モデルをでっち上げることに成功しました。

実は $c=b, d=a$ のときは，初期値が全て $0$ 以上でどれかが正のとき， $X(t),Y(t),Z(t)$ は $t\to \infty$ で収束し， $\displaystyle{ \lim_{t\to \infty} \cfrac{Z(t)}{X(t)}=1,~\lim_{t\to \infty} \cfrac{Y(t)}{X(t)}=\cfrac{\sqrt{b^2+4ae}-b}{2e}}$ となります。つまり $\cfrac{\sqrt{b^2+4ae}-b}{2e}=3\Leftrightarrow a-3b-9e=0,e\neq 0$ となるように定数 $a,b,e$ を定めると $2:6:2$ のモデルになります。

考察

なんかでっち上げとか言ってきましたが，こんなにうまくいくとひょっとしたらモデルとして合ってるんじゃないの？って気もしてきます。Wikipedia 働きアリの法則の解説の項を読んでみましょう。

北海道大学の長谷川英祐が進化生物学の見地から詳しく研究し、一般向けの解説書を出している。それによると、働くアリと働かないアリの差は「腰の重さ」、専門的に言うと「反応閾値」によるという。アリの前に仕事が現れた時、まず最も閾値の低い（腰の軽い）アリが働き始め、次の仕事が現れた時には次に閾値の低いアリが働く、と言う形で、仕事の分担がなされている。仕事が増えたり、最初から働いていたアリが疲れて休むなどして仕事が回ってくると、それまで仕事をしていなかった反応閾値の高い（腰の重い）アリが代わりに働きだす。

あれ全然違うこと言ってる。てへぺろ☆(・ω<)
都合よく作ったモデルで都合のいい結果が出たからと言って，気持ちよくなって正しいと思い込まないように注意しましょう。