「データの見えざる手」のU分布を、Rでシミュレート

さて、本筋とは関係のないところでケチばっかりつけていた、↓前回と前々回の記事でしたが、

「データの見えざる手」の正規分布の図が一様分布に見えたのでRで試した

「データの見えざる手」の図が分かりにくかったので、Rで一次元プロット

今回は、著者の矢野和男さんの言うところの「U分布」なるものを、計算機シミュレーションで作り出してみましょう。

次に、このようにランダムに玉を分配した後で、マス目間で玉をやりとりさせてみよう。
ランダムにマス目を二つ選んで、一方から他方に玉を１個移す。そして、これを繰り返してみよう。もともと、ランダムに置いた玉なのだから、そこからランダムにマス目を選んで、玉を動かしても、結果は変わらない、と思うだろう。この問題を多くの人に出題してみたが、全員が「結果は変わらない」と答えた。

たしかに、直感的には、ランダムに配置後にランダムに交換しても、マクロな状況は変わらないような気がしますね。でも、そうじゃないのが興味深いところ。

書籍では初期値はランダムとありましたが、手を抜いて、1マス80個の「平等」状態からスタートさせてみました。（結果は同じになりますよね、たぶん）

指定回数の交換を行ったあと、それぞれのマスが持っている個数でソートして、少ない方が左になるようにプロットしています。

1万回、2万回、・・・と実行しながら、プロット結果を画像として出力していきます。jのところのループ回数を変えて、1億回まで実行してみました。


n <- 72000  # 点の個数
m <- 900    # マスの個数
masu <- rep(n/m, m) # 平等に配分

for( i in 1:9 ) {
  for( j in 1:10000 ) {
    s <- sample(1:m, 2) # ランダムに2つのマスを選ぶ
    if( masu[s[1]] != 0 ) { # 無い袖は振れないケースへの対処
      masu[s[1]] <- masu[s[1]] - 1 # 1つ目のマスから取って、
      masu[s[2]] <- masu[s[2]] + 1 # 2つ目のマスへ入れる
    }
  }
  s1 <- paste(i, "万回（リニア）", sep="")
  s2 <- paste(i, "万回（対数）", sep="")

  # リニアでプロット
  png(paste(s1,"png",sep="."))
  plot(masu[order(masu)], main=s1)
  dev.off()

  # y軸を片対数でプロット
  png(paste(s2,"png",sep="."))
  plot(log(masu[order(masu)],2), main=s2)
  dev.off()
}


まず、１万回↓

1万回 リニアプロット

全然、指数関数っぽくないですよね。まだ「交換」が十分に行われていないので、80個付近のマスが多いため、真ん中あたりが踊り場のようになっています。そんな中でも、0個近くの負け組と200個以上の勝ち組が現れてきています。

次、10万回↓

10万回 リニアプロット

あまり変わらないですね。（機械的にやったのでグラフのタイトル表記が変ですが、気にしないでね）

そして、100万回↓

100万回 リニアプロット

踊り場が無くなってきました。

さらに、1000万回↓

1000万回 リニアプロット

おお、指数関数っぽくなってきました。

最後に、1億回

1億回 リニアプロット

1000万回とあまり変わっていません、収束してきたのかな。


さて、曲線を見る限りは指数関数っぽいですが、対数をとってグラフが直線になっているのかどうか見てみましょう。

1万回 対数プロット

10万回 対数プロット

100万回 対数プロット

1000万回 対数プロット

1億回 対数プロット

う～ん、直線にはなりませんねえ。ということは、このシミュレーション結果は厳密には指数関数にはならないものなのでしょうか、それとも私の組んだロジックに不備があるのでしょうか。

でも、中間のあたりに関しては直線になっていると言えなくもないですね。