Rでゾンビの増え方のシミュレーションを雑にしてみた

ゾンビ化する人の数は？

柳田理科雄氏のこの記事を読んで、ゾンビの増え方について考えてみた。

news.yahoo.co.jp

1体のゾンビが1人に噛みつくと、ゾンビは2体になる。

その2体がそれぞれ1人を噛むと4体になる。

その4体がそれぞれ1人を襲って8体に、さらに16体に、32体に……とゾンビは倍増していく。

ゾンビの恐ろしさは、この「倍々」という増え方で、これによって短時間でゾンビだらけになる。

とゾンビの増え方が定義（？）されており、

ゆきたちが暮らす巡ヶ丘市の人口を10万人と仮定して、市民全員がゾンビになる過程を考えてみよう。

最初の1体が1人に噛みつくのを「1段階」とすれば、前述の32体になるのが5段階め。それ以降のゾンビ数を記すと、こうなる。

6段階……64人

7段階……128人

8段階……256人

9段階……512人

10段階……1024人

（中略）

16段階……6万5536人

続く17段階では、13万1072人と、人口10万人を超えてしまうから、この時点で市民の総ゾンビ化は完了することになる。

と記述されている。

これは単に、n段階目のゾンビ数は2のn乗の結果と考えているとみられる。

つまり、10段階でのソンビ数は2の10乗で1024人ということ。

そこから、人口10万人の都市が完全にゾンビ化するのには17段階目で完了するとしている。

これは前提として、最初に柳田氏が示した条件に加えて「必ず次の段階では前の段階でゾンビになっている人すべてが一人ずつ、ゾンビになっていない人を一人ずつみつけてゾンビ化する」という条件が必要となる。

こういう条件ならそれもそれでいいと思うが、次に閉鎖空間（学校）での例をあげ、

ここまでは、巡ヶ丘市の状況である。では、ゆきたちの高校はどうなのか？

学校は外との出入りがない「閉鎖空間」だ。 （中略） ゾンビが増えると、非ゾンビ人口が減ってきて、噛んでもゾンビが増えないケースが出てくる（ゾンビは「共噛み」もする）からだ。

巡ヶ丘学院高校の生徒と教職員の数が千人だったと仮定しよう。

その場合、巡ヶ丘市が完全ゾンビ化した「17段階」の時点では、学校内でゾンビになってしまった者は993人。

（中略） 18段階……996人（残された人間は4人）

19段階……998人（残された人間は2人）

20段階……999人（残された人間は1人）

なんとなんと、20段階の時点でも、たった1人だけゾンビ化を免れている！

としている。

ゾンビがソンビと共噛みする可能性を示唆している。

この条件下だど、次第に非ゾンビを発見しにくくなり、ゾンビ化が鈍化することが示されている。

ここで、単純な2のn乗ではないことは明らかだが、疑問としてどうやって計算したのかということが気になる。

おそらく良いモデルがあるのかもしれないが、難しいことはわからないので、簡単なシミュレーションをRでしてみることにする。

追記

SIRモデルは、ちょっと違うかなぁと思っていたら調べたらSZDモデルというSIRを拡張したモデルがあった。

journals.aps.org

S は感受性の高い集団、つまり感染していない人間を表している。 Z は感染状態であるゾンビ、そして Rは除去された状態、この場合は人間によって除去されたゾンビを表しています（正規には脳を破壊して動作不能にする）。人間はゾンビに噛まれると感染し、ゾンビは人間の直接行動によって破壊される。

面白い。

追記終了

ゾンビ数のRでシミュレーション

条件をきちんと整理する。

• 最初は1人のゾンビから始まる。
• ゾンビに噛まれるとゾンビになる。
• ゾンビは共噛みをする可能性がある（ゾンビ数は増えず）。
• 人数は有限。
• ゾンビが次に誰を噛むのかは完全にランダム（無作為）。

この条件で考えると、単純化すると次のように考えられる。

• 人口-1分の長さで値がすべて0でできているベクトルを考える。
• このベクトルに1（最初のゾンビ）を結合する。
• 1以上の値をもつものがゾンビと定義する。
• 1以上の値の個数だけ、無作為にベクトル中の値に1を加える。例：c(1, 1, 0, 0)の場合、1以上の値が2個あるのでだれかが2回1を加えられる。これは0の人とは限らず、1の値を持つ人が1を加えられる場合もありえる。
• ただし、自分で自分を噛むことはないので、自分の値の影響で自分の値が+1することはない。例：c(1, 1, 0, 0)の場合、1つ目の値の影響で2つ目の値に1が加えられる可能性はあるが、1つ目の値の影響で1つ目の値に1が加えられることはないようにする。
• ゾンビは必ず誰かを次の段階で噛む。必ずどこかの値をゾンビの数だけ+1する。
• 複数のゾンビが一人（非ゾンビ、ゾンビ関わらず）を噛むのは妨げない。

とした場合、1000人の学校（閉鎖空間）を想定した場合、100回シミュレーションするとしたら次のようになる。

# 結果を保存するオブジェクト
# 何段階目で全員ゾンビになったか
simulate_times <- NULL
# 各段階での非ゾンビ数
simulate_save_n <- NULL
# 1000回ループする
for(i in 1:100){
  # 非ゾンビ999人
  city <- c(rep(0, 999))
  # ゾンビ1人
  zombie <- 1
  # ゾンビ1人を非ゾンビの中に投入し1000人の閉鎖空間を考える。
  zombie_city <- c(zombie, city)
  # ゾンビを含めた集団の人数
  city_n <- length(zombie_city)
  # 各シミュレーション実行回ごとの各段階での生存者数を残しておくオブジェクト
  save_n <- NULL
  # 各シミュレーション実行回ごとの全員がゾンビになった段階数
  times <- 0
  # 集団が全員ゾンビになる、つまり1以上の値を取るまでループさせる
  while(any(zombie_city == 0)){
    # 集団におけるゾンビの数を計算する
    zombie_n <- sum(zombie_city >= 1)
    # 一番最初はゾンビが一人で、ベクトルの一番最初につけているので、自分以外を噛む
    # つまり、2番目以降の誰かを噛むので2番目以降のだれを噛むか決定する
    if(zombie_n == 1) next_zombie <- sample(2:city_n, size = 1, replace = TRUE)
    # ゾンビが2人以上いる場合は、誰を噛むのかを決定する。
    # 自分を噛まないようにする。
    else{
      # カウンター
      j <- 1
      # ゾンビ化している人の位置を求めて、その人の数だけループする
      for(k in which(zombie_city >= 1)){
        # 次にゾンビ化する人を決める
        # 自分自身を除く全ての人から一人選ぶ（ゾンビ、非ゾンビは区別しない）
        # カウンターをつかって、値を格納していく
        next_zombie[j] <- sample(c(1:city_n)[-k], size = 1, replace = TRUE)
        # カウンターに1を加えてもう一度ループする
        j <- j + 1
      }
    } 
    # 先程のソンビ化する人の位置に+1を加える
    zombie_city[next_zombie] <- zombie_city[next_zombie] + 1
    # 段階数に1を加える
    times <- times + 1
    # 生存者数を段階数ごとの位置に入れておく
    save_n[times] <- sum(zombie_city == 0)
    # 全員がゾンビ化する（すべての値が1以上になる）までループ
  }
  # 全員がゾンビ化したらループを抜ける
  # その時の結果を回数と生存者数で保存する
  simulate_times[i] <- times
  # 生存者数は、シミュレーションごとに段階数（ベクトルの長さ）が異なるのでリストで残しておく
  simulate_save_n <- c(simulate_save_n, list(save_n))
}

1000人くらいだと100回やっても割りとすぐ。

平均何段階目で全員ゾンビ化したかを見るなら次のようにするとわかる。

mean(simulate_times)

## [1] 17.95

最小値と最大値を確認してみる。

min(simulate_times)

## [1] 16

max(simulate_times)

## [1] 23

平均、最小、最大は実行するたびに微妙に変わるが100回もやってるのでだいたい同じくらい（平均18、最小16、最大24くらい）に落ち着く。

箱ひげ図を書いて分布を確認してみる。

boxplot(simulate_times)

ボリュームゾーンがだいたい17～19段階目にある印象。柳田氏の想定よりも早い。

生存者数をみてみる。

リスト形式だと取り扱いにくいので、データフレーム化する。

長さの異なるリストをデータフレーム化するのはちょっとめんどくさいので、githubに公開している{infun}パッケージのlist2data.frame_rbind()を使う。

github.com

これを使うと長さの異なるリストをデータフレームにでき、かつ、長さが足りない部分はNAとしてくれる。

simulate_save_n_df <- infun::list2data.frame_rbind(simulate_save_n)

ここから生存者数の各段階での平均値を見てみる。

n <- apply(simulate_save_n_df, 2, mean, na.rm = TRUE)
n

##          X1          X2          X3          X4          X5          X6 
## 998.0000000 996.0000000 992.0300000 984.1800000 968.6500000 938.8000000 
##          X7          X8          X9         X10         X11         X12 
## 883.0500000 785.9600000 634.9600000 442.6700000 254.2500000 121.3700000 
##         X13         X14         X15         X16         X17         X18 
##  50.3100000  19.5700000   7.4600000   2.6300000   1.0752688   0.5862069 
##         X19         X20         X21         X22         X23 
##   0.3928571   0.4444444   0.5000000   0.5000000   0.0000000

変化をグラフで可視化してみる。

plot(n, type = "b")

最初、生存者数はたくさんいるが、ある段階になると急激に生存者数が少なく、最後の数段階は数人しか生存していないことがわかる。

次に市全体の場合を考える。

柳田氏は市全体の場合は、単純に指数関数的に増えることを想定しているが、市の場合も、実際には学校と同様に人口は有限であり、（流入流出を想定しなければ）学校の場合と同じ計算ができる。

人口10万人の都市を想定し、100回シミュレーションする。

ところで、ゾンビが各段階で誰かを噛む必要性はないのではないかとも思う。つまり、自分以外に誰かがいればその人を噛むが誰もいなければ噛まない（または自分を噛む）みたいなのでも良いと思う。

どの程度の確率で自分を噛むのかというのを決めるとそれでも計算できそうだが、とりあえず自分を含めた全員を筋斗雲確率で噛むとする。

正直これは、各ゾンビが自分を噛まないという計算をすると時間がかかりすぎるというのもある。

simulate_times <- NULL
simulate_save_n <- NULL

for(i in 1:100){
  # 非ゾンビ99999人
  city <- c(rep(0, 99999))
  # ゾンビ1人
  zombie <- 1
  # ゾンビ1人を非ゾンビの中に投入し100000人の閉鎖空間を考える。
  zombie_city <- c(zombie, city)
  city_n <- length(zombie_city)
  save_n <- NULL
  times <- 0
  while(any(zombie_city == 0)){
    zombie_n <- sum(zombie_city >= 1)
    if(zombie_n == 1) next_zombie <- sample(2:city_n, size = 1, replace = TRUE)
    else next_zombie <- sample(1:city_n, size = zombie_n, replace = TRUE)
    zombie_city[next_zombie] <- zombie_city[next_zombie] + 1
    times <- times + 1
    save_n[times] <- sum(zombie_city == 0)
  }
  simulate_times[i] <- times
  simulate_save_n <- c(simulate_save_n, list(save_n))
}

ここでも同様に何段階で全員ゾンビ化したのか平均、最大値、最小値を求める。

mean(simulate_times)

## [1] 29.4

min(simulate_times)

## [1] 27

max(simulate_times)

## [1] 34

箱ひげ図を書いてみる。

boxplot(simulate_times)

ランダム要素があるので実行するたびに結果は変わるが、概ね柳田氏が想定したよりも全員がゾンビ化するまでは時間がかかりそう。

各段階の生存者数の平均値も求める。

simulate_save_n_df <- infun::list2data.frame_rbind(simulate_save_n)
n <- apply(simulate_save_n_df, 2, mean, na.rm = TRUE)
n

##           X1           X2           X3           X4           X5           X6 
## 9.999800e+04 9.999600e+04 9.999200e+04 9.998400e+04 9.996800e+04 9.993600e+04 
##           X7           X8           X9          X10          X11          X12 
## 9.987206e+04 9.974426e+04 9.948947e+04 9.898278e+04 9.798128e+04 9.602345e+04 
##          X13          X14          X15          X16          X17          X18 
## 9.228215e+04 8.543413e+04 7.385806e+04 5.687156e+04 3.694566e+04 1.968380e+04 
##          X19          X20          X21          X22          X23          X24 
## 8.824310e+03 3.541230e+03 1.352170e+03 5.028700e+02 1.844300e+02 6.798000e+01 
##          X25          X26          X27          X28          X29          X30 
## 2.525000e+01 9.570000e+00 3.530000e+00 1.353535e+00 7.534247e-01 3.636364e-01 
##          X31          X32          X33          X34 
## 4.000000e-01 3.333333e-01 5.000000e-01 0.000000e+00

変化をグラフで可視化してみる。

plot(n, type = "b")

ただ、最初に書いたとおり、今回は条件としてゾンビが誰を噛むのかは完全にランダムなので距離も何も考えていない、通常はゾンビは近い人から噛むと思うし、自分の周りがゾンビだらけになったゾンビは新たにゾンビを生み出すゾンビにはなり得ないなどの条件が存在し得ると思う。そもそもゾンビ側の条件（性能？）も不明だし。

調べれば良いモデルがあるのかもしれない。

追記

調べてみたらSIRモデルやその拡張を使って計算しているものが複数あった。

結局どのモデルを使うか、結局どの程度広がるのか、実際には観測し得ないので条件をどうするのかが難しそうだと思った。空間の問題もあると思った。

追記終了