表が出る確率が0.5の硬貨を10枚投げたとき、表が $r$ 枚出る確率 $q_r$ を $r = 0,1,2,\ldots,10$ について求めるには、次のようにします(2項検定参照)。
from scipy import stats stats.binom.pmf(range(11), 10, 0.5)
array([0.00097656, 0.00976563, 0.04394531, 0.1171875 , 0.20507812,
0.24609375, 0.20507812, 0.1171875 , 0.04394531, 0.00976563,
0.00097656])
対称性 $q_r = q_{10-r}$ を満たしているようです。
さて、「表が $r$ 枚出た」という事象の(両側)$p$ 値は、$q_i \leq q_r$ を満たすすべての $q_i$ の和と定義されます。これを $r = 0,1,2,\ldots,10$ について求めてみましょう。
import numpy as np q = stats.binom.pmf(range(11), 10, 0.5) np.array([q[q <= q[r]].sum() for r in range(11)])
array([9.76562500e-04, 2.14843750e-02, 1.09375000e-01, 3.43750000e-01,
5.48828125e-01, 1.00000000e+00, 7.53906250e-01, 3.43750000e-01,
1.09375000e-01, 2.14843750e-02, 1.95312500e-03])
あれ? 対称であるはずの $p$ 値が対称でなくなってしまいました。
np.array にすると丸められた値が表示されるのですが、実際は数値計算の誤差があったのです。
list(q)
[ 0.0009765624999999993, 0.009765625000000009, 0.043945312500000076, 0.11718750000000004, 0.20507812499999997, 0.2460937500000002, 0.205078125, 0.11718750000000004, 0.043945312500000076, 0.009765625000000009, 0.0009765625 ]
このような状況で、<= を使って比較することは危険です。ほんの少し違っていることも許容するように、右辺に非常に小さい値を加えましょう。よく使われる値は np.finfo(np.float64).eps(または同じことですが sys.float_info.epsilon)で、$2.22 \times 10^{-16}$ ほどの値です。
np.array([q[q <= q[r] + np.finfo(np.float64).eps].sum() for r in range(11)])
array([0.00195312, 0.02148438, 0.109375 , 0.34375 , 0.75390625,
1. , 0.75390625, 0.34375 , 0.109375 , 0.02148438,
0.00195312])
うまくいきました。これが正しい値です。
場合によっては、右辺をそれより大きい最小の値で置き換えることがあります。
np.array([q[q <= np.nextafter(q[r], np.inf)].sum() for r in range(11)])
array([9.7656250e-04, 2.1484375e-02, 1.0937500e-01, 3.4375000e-01,
7.5390625e-01, 1.0000000e+00, 7.5390625e-01, 3.4375000e-01,
1.0937500e-01, 2.1484375e-02, 1.9531250e-03])
でも、これでは一番最初の値が正しくありません。
ちなみに、先ほどの np.finfo(np.float64).eps は1より大きい最小の値と1との差に相当し、機械エプシロン(machine epsilon、計算機エプシロン、マシンエプシロン)と呼ばれます。
np.nextafter(1, np.inf) - 1
np.float64(2.220446049250313e-16)
さらに、2つの数が近いかどうかを判断する np.isclose() を使う方法もあります。
np.array([q[(q < q[r]) | np.isclose(q, q[r])].sum() for r in range(11)])
これでも正しい値が得られました。np.isclose(a, b, rtol=1e-05, atol=1e-08) は
np.absolute(a - b) <= atol + rtol * np.absolute(b)
と同じ意味です。デフォルトの相対許容誤差 rtol、絶対許容誤差 atol は大きめに設定されているので、もっと小さくていいかもしれません。