二封筒問題のシミュレーション その３ 整数拡張編

■

概要

二封筒問題(参: 三浦俊彦『可能世界の哲学』)を Python を用いたシミュレーションで解く。その１では未開封と開封で答えが違うことがシミュレーションによっても確かめられた。その謎を解くため、その２では、まずは額を実数にするよう問題を拡張した。

この「その３」では整数問題に戻って拡張を試み、謎の解明を完成する。

■

整数拡張

その２の実数拡張では、額の上限が与える影響を示すことができた。整数ではさらに、額が偶数か奇数かがさらなる問題となる。具体的には、開封して奇数であれば、それは相手の二倍でないのは確実なので、交換するのがどうみても得になるわけである。

それを解析していこう。シミュレーションは実数拡張と似ているが、期待値(平均値)を直接求めるのではなく、試行をすべてメモリにたくわえてあとで処理するというある意味 Python らしいメモリをたくさん使った実装にしてみた。

プログラムは次の two_envelopes_4.py である。parse_args をはじめいろいろ省略している。ARGS.x_max すなわち X は、偶数でなければならないとする。

import random
import numpy as np
from sympy import *
import argparse
ARGS = argparse.Namespace()

ARGS.trials = 10000
ARGS.x_max = 200

# parse_args 省略

def main ():
    trials = []
    for trial in range(ARGS.trials):
        r = random.randrange(2)
        x = random.randrange(1, ARGS.x_max + 1)
        if r == 0:
            a = x * 2
            b = x
        else:
            a = x
            b = x * 2
        trials.append([a, b, r, x, b / a])
    trials = np.array(trials)
    t_a = trials[:, 0]
    t_b = trials[:, 1]
    t_r = trials[:, 2]
    t_x = trials[:, 3]
    t_a1 = trials[(t_x <= ARGS.x_max / 2) & (t_r == 0)]
    t_a2 = trials[(t_x > ARGS.x_max / 2) & (t_r == 0)]
    t_be = trials[(t_r == 1) & (t_x % 2 == 0)]
    t_bo = trials[(t_r == 1) & (t_x % 2 == 1)]
    t_c = trials[(t_a <= ARGS.x_max) & (t_a % 2 == 0)]
    t_d = trials[t_a <= ARGS.x_max]
    t_a2_ = trials[t_a > ARGS.x_max]
    t_bo_ = trials[t_a % 2 == 1]
    assert t_a1.shape[0] + t_be.shape[0] == t_c.shape[0]
    assert t_c.shape[0] + t_bo.shape[0] == t_d.shape[0]
    assert np.all(t_a2 == t_a2_)
    assert np.all(t_bo == t_bo_)

    x = Symbol('x', integer=True)
    X = Symbol('X', integer=True)
    y = Symbol('y', integer=True)

    n_a1 = summation(1, (x, 1, X/2))
    n_a2 = summation(1, (x, X/2 + 1, X))
    n_be = summation(1, (y, 1, X/2))
    n_bo = summation(1, (y, 1, X/2))

    E_a1_b = summation(x, (x, 1, X/2)) / n_a1
    E_a2_b = summation(x, (x, X/2 + 1, X)) / n_a2
    E_be_b = summation((2 * x).subs([(x, 2 * y)]), (y, 1, X / 2)) / n_be
    E_bo_b = summation((2 * x).subs([(x, 2 * y - 1)]), (y, 1, X / 2)) / n_bo
    E_c_b = (n_a1 / (n_a1 + n_be)) * E_a1_b + (n_be / (n_a1 + n_be)) * E_be_b
    E_d_b = (n_a1 / (n_a1 + n_be + n_bo)) * E_a1_b \
        + (n_be / (n_a1 + n_be + n_bo)) * E_be_b \
        + (n_bo / (n_a1 + n_be + n_bo)) * E_bo_b
    E_omega_b = ((n_a1 + n_be) / (n_a1 + n_a2 + n_be + n_bo)) * E_c_b \
        + (n_a2 / (n_a1 + n_a2 + n_be + n_bo)) * E_a2_b \
        + (n_bo / (n_a1 + n_a2 + n_be + n_bo)) * E_bo_b
    assert Eq(E_a1_b, X/4 + 1/2).simplify()
    assert Eq(E_a2_b, 3*X/4 + 1/2).simplify()
    assert Eq(E_be_b, X + 2).simplify()
    assert Eq(E_bo_b, X).simplify()
    assert Eq(E_c_b, 5*X/8 + 5/4).simplify()
    assert Eq(E_d_b, 3*X/4 + 5/6).simplify()
    assert Eq(E_omega_b, 3*X/4 + 3/4).simplify()

# E_c_a, E_d_ba などを求めるのを省略。

    f = lambda q: q.subs([(X, ARGS.x_max)]).simplify().evalf()

    print("Realm: c")
    av_c = np.mean(t_c, axis=0)
    print("valid:", t_c.shape[0], "/", ARGS.trials)
    print("E(A) =", av_c[0], "(==", f(E_c_a), ")")
    print("E(B) =", av_c[1], "(==", f(E_c_b), ")")
    print("E(X) =", av_c[3])
    print("E(B/A) =", av_c[4], "(==", f(E_c_ba), ")")
    print("E(B)/E(A) =", av_c[1] / av_c[0], "(==", f(E_c_b / E_c_a), ")")

# 領域 a2, bo, d に関する print を省略

    assert Eq(E_d_b/E_d_a, (9 * X + 10) / (6 * X + 8)).simplify()
    #plot(E_d_b/E_d_a, (X, 0, 100))

# 領域 Omega に関する print を省略

if __name__ == '__main__':
    # parse_args 省略
    main()

ソースの説明の前に実行結果は次のようになる。

$ python two_envelopes_4.py
Realm: c
valid: 5040 / 10000
E(A) = 101.04166666666667 (== 101.000000000000 )
E(B) = 125.7 (== 126.250000000000 )
E(X) = 75.58055555555555
E(B/A) = 1.2407738095238094 (== 1.25000000000000 )
E(B)/E(A) = 1.244041237113402 (== 1.25000000000000 )

Realm: a2
valid: 2494 / 10000
E(A) = 301.0208500400962 (== 301.000000000000 )
E(B) = 150.5104250200481 (== 150.500000000000 )
E(X) = 150.5104250200481
E(B/A) = 0.5 (== 0.500000000000000 )
E(B)/E(A) = 0.5 (== 0.500000000000000 )

Realm: bo
valid: 2466 / 10000
E(A) = 99.07948094079481 (== 100.000000000000 )
E(B) = 198.15896188158962 (== 200.000000000000 )
E(X) = 99.07948094079481
E(B/A) = 2.0 (== 2.00000000000000 )
E(B)/E(A) = 2.0 (== 2.00000000000000 )

Realm: d
valid: 7506 / 10000
E(A) = 100.39701572075673 (== 100.666666666667 )
E(B) = 149.5054622968292 (== 150.833333333333 )
E(X) = 83.30082600586198
E(B/A) = 1.4902078337330136 (== 1.50000000000000 )
E(B)/E(A) = 1.4891424931659545 (== 1.49834437086093 )

Realm: Omega
valid: 10000
E(A) = 150.4326 (== 150.750000000000 )
E(B) = 149.7561 (== 150.750000000000 )
E(X) = 100.0629
E(B/A) = 1.24325 (== 1.25000000000000 )
E(B)/E(A) = 0.9955029694361461 (== 1.00000000000000 )

ソースの説明をする。

最初に for 文でシミュレーションをしている。これまでの違いは、試行を足し合わせるのではなく trials というリストに残している点である。

そして後から、numpy の機能を使って、trials を領域の条件を満たすかどうかで分類している。基本となる領域 a1, a2, be, bo は次のようになる。なお、X は ARGS.x_max のことである。r == 0 で a が二倍、r == 1 で b が二倍になる。

● 領域 a1: r == 0 で、x <= X/2。 ● 領域 a2: r == 0 で、x > X/2。 ● 領域 be: r == 1 で、x が偶数(even)。 ● 領域 bo: r == 1 で、x が奇数(odd)。

そしてそれらの組み合わせと等しくなるのが領域 c, d, Omega である。

● 領域 c: a <= X かつ a が偶数。 ● 領域 d: a <= X。 ● 領域 Omega: 全領域。

なお、領域 c は 領域 a1 と be を組み合わせたものになり、領域 d は 領域 c と bo を組み合わせたものになる。また、領域 a2 は、a > X の領域と合致し、領域 bo は a が奇数の領域と合致する。それは、assert で結果から確かめてもいる。

次にそれぞれの領域における b と a と b/a の期待値を求めている。

例えば、領域 c における b に交換したときの期待値は E_c_b として式を建てている。それは、期待値の分割を使って E_a1_b と E_be_b から求められる。

  E_c_b = (n_a1 / (n_a1 + n_be)) * E_a1_b + (n_be / (n_a1 + n_be)) * E_be_b

そして例えば領域 a1 における b の期待値は E_a1_b は、n_a1 が領域 a1 のすべての場合の数だとして

  E_a1_b = summation(x, (x, 1, X/2)) / n_a1

…になる。これは assert 文でわかるように結局、E_a1_b == X/4 + 1/2 となる。assert の読み方はその２で説明したのでそちらを参照していただきたい。そして、E_c_b == 5*X/8 + 5/4 となる。

そして最後に a の開封によって観測できる各領域ごとに E(a), E(b), E(x), E(b/a), E(b)/E(a) を表示している。

蛇足だが、領域 d における E(B)/E(A) はほぼ 1.5 になるが、厳密には 1.5 に満たなくなるようだ。E_d_b/E_d_a == (9 * X + 10) / (6 * X + 8) となる。これを plot すると次のようになる。100 ですでに 1.5 近いから 10000 あたりになれば、ほぼ 1.5 になるのだが、ご参考まで。

以上が、two_envelopes_4.py の説明である。

その結果、領域 c では、b への交換が 1.25 倍有利になっている。しかし、領域 a2 では逆に 0.5 倍となり不利になっている。また、領域 bo では 2.0 倍になり 1.25 倍よりさらに有利になっている。結果全領域ではいろいろな効果が打ち消しあって 1.0 倍…すなわち変わらないということになるが、これはかなり偶然が重なってのことであるという印象を受ける。

■	「二封筒問題」の結論

さて、開封して 1万円だとわかった場合という二封筒問題の開封バージョンの条件は、一目して受ける印象と違ってかなり条件が制限されたものだというのがわかる。

それは領域 c に属するということ、つまり、最初に決める額の上限額 X (本当の上限は 2*X) より開封して見た額 a が小さく、かつ、偶数の額であるということ…それが 1.25 倍という結果を導いているのである。

もし、上限 X が例えば 9000円であれば、開封した a が1万円になることは可能だが、交換は 1.25 倍とはならならず 0.5 倍となる。そういう不利な条件がまず排除されている。

さらに開封した額 a が 9999円のような奇数なら交換はからなず 2 倍となる。そういう有利な条件もまた排除されているということになる。

未開封なら 1倍というところから、そういう有利なところと不利なところを刈りこんだら 1.25 倍の領域が残ったというのが、この問題の謎の解になろう。

排除された部分に 10002円がないとして、1万円だったときに交換する者は、 10002円だからと言って交換しない理由はないだろうし、1万円を見て交換しない者は、10002円なら交換する理由はあるまい。そういう意味では1万円で交換する者の期待値は、領域 c についての期待値になると見ても良いと思われる。

ただ、厳密には、領域 c にあることと、1万円だとわかることには違いがある。開封して 1万円だとわかったときの期待値の計算は、三浦の最初の計算で正しいのだろう。謎は、開封した場合ではなく未開封の場合で、それは今回領域 Omega として計算したようになる。…ということなのではないか。

おまけに two_envelopes_5.py というプログラムを付けておいた。これは a <= ARGS.q なる領域 d と a == ARGS.q なる領域 q について期待値を計算するものである。あわよくば a <= ARGS.q に関する期待値の差分を取れば、a == ARGS.q の期待値のヒントが得られるかもと思ったが、E_d_b_ の計算結果を見る限りまったくのあてはずれだった。領域 q の計算は三浦の最初の計算とほぼ同じものにするしかなかった。

■

今後に向けて

その２と同じく、中身を見れず額がわからなくとも、上限額よりはるかに小さい額しか入っていないという情報があれば、交換が有利になるとは言えるだろう。額が偶数だろうが奇数だろうがそれが言える…と。ただ、分布が一様でないなら難しくはなるが…。

その２の最後でも書いたが、ところで「開封」とは、どういう意味を数学・論理学的に持つのだろうか？

開封によって知ることができ、かつ、その知ったことによって変化させることができる事象というものが、実は、はっきり列挙されているわけではない。知ること／知らないことをどう表すか…それが問われていると言えるかもしれない。

領域の分け方はこれ以外ありえないのか、それは証明されていない。この問題に関しては、これ以上の場合わけの方法はないと私には思われるが、それで十分という保証があるわけではない。

領域 c に属すことと、1万円とわかることには違いがある。領域 c に属す各要素について、同じ期待値を保証するということができれば、例えば特定のサイズの魚が池にいた確率(参: 三浦俊彦『ゼロからの論証』)のような実数値の議論にも応用できるのかもしれない。それはどうすればできるのか、私はその方法を知っていない。その一般的性質・どういうときできるかも興味がある。それを知るのも今後の課題だろう。

■	参考＆著者＆ライセンス＆配布物

この記事には、その１ 基礎編、その２ 実数拡張編、その３ 整数拡張編があるが、その１を参照のこと。

更新：	2020-11-05--2020-11-12
初公開：	2020年11月05日 23:05:33
最新版：	2020年11月12日 23:02:02