続・サイコロゲームと貧富の差

# はじめに

[サイコロゲームと貧富の差](https://qiita.com/kaityo256/items/69668bf5a92d538c24ae)という記事で、こんなゲームを紹介しました。

1. 最初、6人がリンゴを一つずつ持っている
2. 6人はそれぞれ1〜6までの番号が振られている
3. サイコロを二回振り、最初に出た目の人が、次に出た目の人にリンゴを一つあげる。ただし、最初の人がリンゴを持っていない場合はそのまま
4. 以上を繰り返す

これが等重率を仮定することでボルツマン分布になること示したのですが、

* 等重率の仮定が何を指すかが不明瞭
* 厳密解との比較がない

という問題がありました。この記事では上記についてもう少し踏み込んだ内容を書いてみようと思います。

# 遷移行列

この問題はマルコフ過程とみなすことができますが、リンゴをk個持っている人の人数を$p_k$として、$p_k$に関する遷移過程を書こうとすると非線形項がでてきてうまくいきません。しかし、[tsekineさん](https://twitter.com/tsekine)の指摘の通り、「何人が何個のリンゴを持っているか」という状態を基底にとると遷移過程が線形になり、遷移行列を書き下すことができます。例えば、リンゴをみんなが1つずつ持っている状態を(1,1,1,1,1,1)、一人が全部もっている状態を(0,0,0,0,0,6)と表記します[^sort]。すると、この系がとり得る状態は以下の11通りあります。

[^sort]: 後で`sort`を使う関係で昇順にしておきます。

```rb
[0,0,0,0,0,6]
[0,0,0,0,1,5]
[0,0,0,0,2,4]
[0,0,0,1,1,4]
[0,0,0,0,3,3]
[0,0,0,1,2,3]
[0,0,1,1,1,3]
[0,0,0,2,2,2]
[0,0,1,1,2,2]
[0,1,1,1,1,2]
[1,1,1,1,1,1]
```

これらそれぞれに0から10までの状態番号を振ると、11×11の遷移行列を書き下すことができます。

例えば現在の状態が[0,0,0,0,0,6]であった時、状態が変化するのはリンゴを6つ持っている人がリンゴを上げる人に選ばれ,かつリンゴをもっていない人がもらう人に選ばれる時だけです。従って、サイコロの目36通りのうち、[0,0,0,0,0,0,6]の状態にとどまるのは31通り、[0,0,0,0,1,5]に移るのが5通りとなります。全部考えるのは鬱陶しいのでスクリプトにやらせましょう。手抜きですがこんな感じにかけるでしょうか。

```rb:test.rb
def makematrix(h,a,m)
  ma = h[a]
  6.times do |i|
    6.times do |j|
      b = a.clone
      if b[i]!=0
        b[i] = b[i] - 1
        b[j] = b[j] + 1
      end
      mb = h[b.sort]
      m[mb*11+ma] = m[mb*11+ma] + 1
    end
  end
end

a = Array.new
a.push [0,0,0,0,0,6]
a.push [0,0,0,0,1,5]
a.push [0,0,0,0,2,4]
a.push [0,0,0,1,1,4]
a.push [0,0,0,0,3,3]
a.push [0,0,0,1,2,3]
a.push [0,0,1,1,1,3]
a.push [0,0,0,2,2,2]
a.push [0,0,1,1,2,2]
a.push [0,1,1,1,1,2]
a.push [1,1,1,1,1,1]

h = Hash.new
a.size.times do |i|
  h[a[i]] = i
end

m = Array.new(a.size**2){0}
a.size.times do |i|
  makematrix(h,a[i],m)
end

m.each_with_index do |v,i|
  puts if i % 11 == 0
  print "#{v} "
end
puts 
```

実行結果はこんな感じになります。

```shell-session
$ ruby test.rb

31 1 0 0 0 0 0 0 0 0 0 
5 30 1 2 0 0 0 0 0 0 0 
0 1 26 2 2 1 0 0 0 0 0 
0 4 4 27 0 1 3 0 0 0 0 
0 0 1 0 26 1 0 0 0 0 0 
0 0 4 2 8 26 6 6 4 0 0 
0 0 0 3 0 3 22 0 2 4 0 
0 0 0 0 0 1 0 21 2 0 0 
0 0 0 0 0 3 3 9 24 12 0 
0 0 0 0 0 0 2 0 4 19 30 
0 0 0 0 0 0 0 0 0 1 6 

```

# 定常状態

遷移行列としては各要素を36で割るべきですが、固有値や固有ベクトルは変わらないので気にしないことにします。これはマルコフ行列であるため、最大固有値は必ず1となり、それに対応する固有ベクトルが定常状態となります。

さて、これを例えばMathematicaに突っ込むと、最大固有値1に対応する固有ベクトルは[6, 30, 30, 60, 15, 120, 60, 20, 90, 30, 1]となります。これは、例えば[0,0,0,0,0,6]の状態に対応する重みが6、[0,0,0,0,1,5]に対応する重みが30...となります。これは、それぞれのリンゴの数の配分の場合の数となります。

これは6つの箱があり、6つのリンゴをわける時、リンゴを区別せずに入れるやり方全てが等しい確率で出現する、ということを意味します。これがこの系における「等重率」の意味です。

この固有ベクトルを使うことで「リンゴをx個持っている人の人数の期待値」はこんな感じに求められるでしょう。

```rb:test2.rb
a = Array.new
a.push [0,0,0,0,0,6]
a.push [0,0,0,0,1,5]
a.push [0,0,0,0,2,4]
a.push [0,0,0,1,1,4]
a.push [0,0,0,0,3,3]
a.push [0,0,0,1,2,3]
a.push [0,0,1,1,1,3]
a.push [0,0,0,2,2,2]
a.push [0,0,1,1,2,2]
a.push [0,1,1,1,1,2]
a.push [1,1,1,1,1,1]

s = Array.new(7){0}

eigs = [6, 30, 30, 60, 15, 120, 60, 20, 90, 30, 1]
sum = eigs.inject{|v,i| v+i}
eigs.map!{|v| v.to_f/sum}

a.size.times do |j|
  7.times do |i|
    s[i] = s[i] + a[j].count(i)*eigs[j]
  end
end

s.each_with_index do |v,i|
  puts "#{i} #{v}"
end
```

実行結果は

```shell-session
$ ruby test2.rb
0 2.727272727272727
1 1.6363636363636365
2 0.9090909090909091
3 0.45454545454545453
4 0.1948051948051948
5 0.06493506493506493
6 0.012987012987012988
```

前回の記事のシミュレーション結果は

```shell-session
0 2.729575
1 1.634563
2 0.907818
3 0.455346
4 0.193464
5 0.065497
6 0.013737
```

でしたので、まぁまぁ合ってますかね。

# まとめ

サイコロゲームについて遷移行列を書き下すことで、定常状態を厳密にもとめてみました。一般に「リンゴを０個もっている人が何人、1個もっている人が何人...」という数字列で状態を指定すると遷移行列を書き下すことができ、かつ最大固有値に対応する固有ベクトルから、「リンゴを6つの箱にわけるやり方」が全て等しい確率で出現することがわかりました。これが「等重率の原理」に対応するものです。

後は任意の人数について(0,0,1,2,3,..)みたいな「リンゴの分け方」の列挙と、その分け方の場合の数を書き下すスクリプトがかければ`test2.rb`みたいな計算で「リンゴをx個持っている人数の期待値」を求めることができます。

さて、遷移行列を書き下すと、たしかに「リンゴを区別せずに、6つの箱にいれるやり方が全て等しい確率で出現する」ことがわかるのですが、このゲームの帰結として「なぜそうなるか？」の直感的な説明は私にはできません。また、これが「なぜ$f \ln f$というエントロピーを最大化した状態として実現されるのかも、簡単な説明は思いつきませんでした。気になった方は身近にいる熱・統計力学に詳しい方に聞いてみてはどうでしょうか。

なおこのモデルは、参加人数が増えると緩和時間が飛躍的に伸びます。$N$が一桁台ならともかく、$N=100$とかの場合、ナイーブなシミュレーションで定常状態を求めるのはかなり厳しいんじゃないかなぁ、と思います。

# 追記:等重率の説明

あまり直感的ではないですが、一応「等重率」が実現されることの説明を思いついたので書いておきます。

今、例えばA,B,C,D,E,Fの6人が、それぞれ「4,1,1,0,0,0」個のリンゴを持っているとします。さて、いま「リンゴあげる人」にAさんが、「リンゴをもらう人」にDさんが選ばれる確率は1/36です。この時、リンゴの数は「3,1,1,1,0,0」になります。

逆に、「3,1,1,1,0,0」の状態から「4,1,1,0,0,0」の状態に戻るには、「リンゴあげる人」にDさんが、「リンゴをもらう人」にAさんが選ばれる場合ですから、やはり1/36です。

ちなみにここでは「人は区別するが、リンゴは区別しない」状態を考えるため、Aさんが6個全部もっている状態「6,0,0,0,0,0」とBさんが6個全部もっている状態「0,6,0,0,0,0」は異なる状態として区別して考えます。

このように「人は区別するが、リンゴは区別しない」状態間の遷移を考えると、任意の二つの状態間の遷移確率が等しくなります。つまり、ある状態Xからある状態Yへ行く確率$P(Y \leftarrow X)$と、その逆の確率$P(X \leftarrow Y)$は等しくなります。さて、定常状態において状態$X, Y$にいる確率をそれぞれ$\pi(X)$、$\pi(Y)$とすると、[詳細釣り合い条件](https://qiita.com/kaityo256/items/f05f9914eb0ad16afe05)から

```math
\pi(X) P(Y \leftarrow X) = \pi(Y) P(X \leftarrow Y)
```

が成り立ちます。いま、$P(Y \leftarrow X) = P(X \leftarrow Y)$ ですから、$\pi(X) = \pi(Y)$であることがわかります。任意の二状態についてこれが示されるため、全ての状態の定常状態における存在確率が等しくなります。これがこの系における「等重率の原理」です。

というわけで一応「等重率」の説明をしてみましたが、いまいち直感的では無い気がしますね。良い説明がありましたらお知らせください。