マシュマロを食べることができるだろうか？

本日(2019/09/08) の日本テレビ，世界の果てまで行ってQ で，100m の高さから 10 個くらいのマシュマロを一度に投げて，下で待ち受ける8人ほどが口に入れることができるだろうかという実験をやっていた。結果は失敗だったけど，どんな条件だと成功する確率がどれくらいになるのかシミュレーションしてみた。

library(MASS)

# 9 人を 40 センチ間隔の格子点に配置（単位はセンチメートル）
x = c(-40, -40, -40, 0, 0, 0, 40, 40, 40)
y = c(-40, 0, 40, -40, 0, 40, -40, 0, 40)
person = length(x)

sigma = 4000 # 二次元正規分布の分散
marshmallow = 10 # 一度に投げるマシュマロの数
trial = 1000 # 実験回数
result = integer(trial) # マシュマロが何人の口に入ったかの記録
for (loop in 1:trial) {
    # 二次元正規乱数の発生
    d = mvrnorm(n=marshmallow, mu=c(0, 0), Sigma=matrix(c(sigma, 0, 0, sigma), 2), empirical=TRUE)
    # plot(d, cex = 0.5)
    # points(x, y, pch=19, col=2, cex = 0.5)
    success = 0 # 何人の口に入ったか
    for (i in 1:marshmallow) { # それぞれのマシュマロの着地点について
        for (j in 1:person) { # 各人の口から半径 5 センチ以内なら口に入るだろうと判定
            if ((d[i, 1] - x[j])^2 + (d[i, 2] - y[j])^2 < 5^2) {
                success = success + 1
            }
        }
    }
    result[loop] = success # 結果の記録
}
table(result) # 結果集計
print(mean(result)) # 実験 1 回あたりの，マシュマロが口に入る回数

マシュマロを10個投げる。下で 9 人が待ち受ける。マシュマロの着地点が，待ち受ける人口の半径5センチ以内なら，口に入れることができるとの条件で計算してみた。

ちなみに，マシュマロの着地点と 9 人の位置関係は，下図のような感じ。ここでは，マシュマロを 100 個投げたときの状況を示す。白抜きの丸がマシュマロ。赤丸が待ち受けている人。これだと，けっこう入りそうな気がする。

結果は，ほぼ 20% の成功確率。

マシュマロの数を倍の 20 個にすると当然（？？）倍くらい（40%）の成功確率となる。

判定条件を各人の口から 3 センチ以内とすると，成功確率は 6% ぐらいになる。

まあ，大まかでもいいけどシミュレーションしておかないと，「成功しました！！」とはなかなか行かないものだということだ。

Python って奴は，なんて自由な奴なんだ！

プログラミングの初歩の段階で，「1 以上，10 以下の場合は 'OK'，そうでない場合は 'NO' を書くプログラムを作りなさい」などというのがあり，初心者の中には，

if (1 < x < 10) print('OK')

などと書きたくなる人がいる。大抵の言語はこれはだめで，

if (1 < x && x < 10) print('OK')

としないといけない。

 

しかし，である。Python はこれを許す。なんですと!! ちっとも知らなかった。

def if_test(x):
    if 1 <= x <= 10:
        print('OK')
    else:
        print('NG')

if_test(0) # NG
if_test(5) # OK
if_test(10) # OK

ビットカウント（ビットが 1 である個数）

「ビットカウントする高速アルゴリズムをPythonで実装しながら詳しく解説してみる」というのを見つけた。

Python の整数は 64 ビット整数，R は 32 ビット整数であることを考慮し，以下のように翻訳した。

bitcount = function(n) {
  # nの立っているビット数をカウントする関数(nは64bit整数)
  # 2bitごとの組に分け、立っているビット数を2bitで表現する
  n = n - bitwAnd(bitwShiftR(n, 1), 0x55555555)
  # 4bit整数に 上位2bit + 下位2bit を計算した値を入れる
  n = bitwAnd(n, 0x33333333) + bitwAnd(bitwShiftR(n, 2), 0x33333333)
  n = bitwAnd((n + bitwShiftR(n, 4)), 0x0f0f0f0f) # 8bitごと
  n = n + bitwShiftR(n, 8) # 16bitごと
  n = n + bitwShiftR(n, 16) # 32bitごと = 全部の合計
  return(bitwAnd(n, 0x0000007f))
}

正しく移植できたか，確認してみる。

> bitcount(7)
[1] 3
> bitcount(127)
[1] 7
> bitcount(2^31-1)
[1] 31
> bitcount(-1)
[1] 32

R の intToBits() は，整数の二進表記を得るものである。これを用いてビット 1 を数えるのは自然である。 intToBits() のクラスは "raw" なので，as.integer() で整数化してから合計する。

bitcount2 = function(n) {
  return(sum(as.integer(intToBits(n))))
}

同じ結果になることが確認できる。

> bitcount2(7)
[1] 3
> bitcount2(127)
[1] 7
> bitcount2(2^31-1)
[1] 31
> bitcount2(-1)
[1] 32

では，両者の速度比較をしてみる。intToBits() がベクトルに対応していないので，心ならずも for ループで検証する。

> mx = 1000000000
> n = sample(mx, 1000000) # 100 万個の整数
> a = integer(mx)
> b = integer(mx)
> system.time({for (i in n) a = bitcount(i)})
   ユーザ   システム       経過  
     6.915      0.198      7.118
> system.time({for (i in n) b = bitcount2(i)})
   ユーザ   システム       経過  
     1.620      0.148      1.769
> all(a == b) # どちらの結果も同じか？
[1] TRUE

結論は，「intToBits() を使う方が 4 倍ほど速い」ということだが，100 万個の整数を処理して 2 〜 7 秒なので，1 個や 100 個なら，どちらの方法でも，ゴミのような差しかない。

======================================

ついでに，Python でも代替関数について調べてみた。

Python には二進表記をするための bin() がある。これは文字列として二進数を返すので，文字列の中に '1' がいくつあるか数える。

大元のビットカウントのプログラム

def bitcount(x):
    '''xの立っているビット数をカウントする関数
    (xは64bit整数)'''
    # 2bitごとの組に分け、立っているビット数を2bitで表現する
    x = x - ((x >> 1) & 0x5555555555555555)
    # 4bit整数に 上位2bit + 下位2bit を計算した値を入れる
    x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)
    x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f # 8bitごと
    x = x + (x >> 8) # 16bitごと
    x = x + (x >> 16) # 32bitごと
    x = x + (x >> 32) # 64bitごと = 全部の合計
    return x & 0x0000007f

対する，bin() を使う関数。短い。

def bitcount2(x):
    return sum([i == '1' for i in bin(x)])

100 万個の整数を対象とする。

import scipy as sp
a = sp.random.randint(-10, 100000000, 1000000)


def func1(): # 4.383 sec.
    x = [bitcount(i) for i in a]

def func2(): # 3.416 sec.
    y = [bitcount2(i) for i in a]

やはり，bin() を使う方が速い。
sp.all(x == y) は True である。