没有替换概率的抽样
Sampling Without Replacement Probabilities
我正在使用 np.random.choice 进行无放回采样。
我希望以下代码选择 0 50% 的时间、1 30% 的时间和 2 20% 的时间。
import numpy as np
draws = []
for _ in range(10000):
draw = np.random.choice(3, size=2, replace=False, p=[0.5, 0.3, 0.2])
draws.append(draw)
result = np.r_[draws]
如何正确选择 np.random.choice 的参数以获得我想要的结果？
我想要的数字代表事件被排在第 1 位或第 2 位的概率。
print(np.any(result==0, axis=1).mean()) # 0.83, want 0.8
print(np.any(result==1, axis=1).mean()) # 0.68, want 0.7
print(np.any(result==2, axis=1).mean()) # 0.47, want 0.5
我对这个问题给出了两种解释。一种我更喜欢 ("Timeless") 一种我认为技术上有效但较差的 ("Naive")
永恒：
给定概率 x, y, z 这种方法计算 x', y', z' 这样如果我们独立绘制两次并丢弃所有相等的对 0, 1, 2 的频率是 x, y, z.
这在两次试验中给出了正确的总频率，并且在第一次和第二次试验等效的意义上具有简单和永恒的额外好处。
要做到这一点，我们必须
(x'y' + x'z') / [2 (x'y' + x'z' + y'z')] = x
(x'y' + y'z') / [2 (x'y' + x'z' + y'z')] = y (1)
(y'z' + x'z') / [2 (x'y' + x'z' + y'z')] = z
如果我们将其中两个相加并减去第三个，我们得到
x'y' / (x'y' + x'z' + y'z') = x + y - z = 1 - 2 z
x'z' / (x'y' + x'z' + y'z') = x - y + z = 1 - 2 y (2)
y'z' / (x'y' + x'z' + y'z') = -x + y + z = 1 - 2 x
将其中的 2 个乘以除以第 3 个
x'^2 / (x'y' + x'z' + y'z') = (1 - 2 z) (1 - 2 y) / (1 - 2 x)
y'^2 / (x'y' + x'z' + y'z') = (1 - 2 z) (1 - 2 x) / (1 - 2 y) (3)
z'^2 / (x'y' + x'z' + y'z') = (1 - 2 x) (1 - 2 y) / (1 - 2 z)
因此高达一个常数因子
x' ~ sqrt[(1 - 2 z) (1 - 2 y) / (1 - 2 x)]
y' ~ sqrt[(1 - 2 z) (1 - 2 x) / (1 - 2 y)] (4)
z' ~ sqrt[(1 - 2 x) (1 - 2 y) / (1 - 2 z)]
因为我们知道 x', y', z' 总和为 1 这足以求解。
但是：我们实际上不需要完全解决 x', y', z'。由于我们只对不等对感兴趣，因此我们只需要条件概率 x'y' / (x'y' + x'z' + y'z')、x'z' / (x'y' + x'z' + y'z') 和 y'z' / (x'y' + x'z' + y'z')。这些我们可以使用等式（2）计算。
然后我们将它们中的每一个减半以获得有序对的概率，并从具有这些概率的六个合法对中抽取。
天真：
这是基于（在我看来是武断的）假设，即在第一次抽签概率 x', y', z' 之后，如果第一次抽签的概率为 0，则第二次抽签的条件概率必须为 0, y' / (y'+z'), z' / (y'+z') x' / (x'+z'), 0, z' / (x'+z') 如果第一个是 1 和概率 x' / (x'+y'), y' / (x'+y'), 0) 如果第一个是 2.
这有一个缺点，据我所知，没有简单的、封闭形式的解决方案，而且第二次和第一次抽签有很大不同。
优点是可以直接用np.random.choice；然而，这太慢了，以至于在下面的实现中我给出了一个避免这个函数的解决方法。
经过一些代数运算后发现：
1/x' - x' = c (1 - 2x)
1/y' - y' = c (1 - 2y)
1/z' - z' = c (1 - 2z)
其中 c = 1/x' + 1/y' + 1/z' - 1。这个我只能用数字来解决。
实施和结果：
这是实现。
import numpy as np
from scipy import optimize
def f_pairs(n, p):
p = np.asanyarray(p)
p /= p.sum()
assert np.all(p <= 0.5)
pp = 1 - 2*p
# the following two lines show how to compute x', y', z'
# pp = np.sqrt(pp.prod()) / pp
# pp /= pp.sum()
# now pp contains x', y', z'
i, j = np.triu_indices(3, 1)
i, j = i[::-1], j[::-1]
pairs = np.c_[np.r_[i, j], np.r_[j, i]]
pp6 = np.r_[pp/2, pp/2]
return pairs[np.random.choice(6, size=(n,), replace=True, p=pp6)]
def f_opt(n, p):
p = np.asanyarray(p)
p /= p.sum()
pp = 1 - 2*p
def target(l):
lp2 = l*pp/2
return (np.sqrt(1 + lp2**2) - lp2).sum() - 1
l = optimize.root(target, 8).x
lp2 = l*pp/2
pp = np.sqrt(1 + lp2**2) - lp2
fst = np.random.choice(3, size=(n,), replace=True, p=pp)
snd = (
(np.random.random((n,)) < (1 / (1 + (pp[(fst+1)%3] / pp[(fst-1)%3]))))
+ fst + 1) % 3
return np.c_[fst, snd]
def f_naive(n, p):
p = np.asanyarray(p)
p /= p.sum()
pp = 1 - 2*p
def target(l):
lp2 = l*pp/2
return (np.sqrt(1 + lp2**2) - lp2).sum() - 1
l = optimize.root(target, 8).x
lp2 = l*pp/2
pp = np.sqrt(1 + lp2**2) - lp2
return np.array([np.random.choice(3, (2,), replace=False, p=pp)
for _ in range(n)])
def check_sol(p, sol):
N = len(sol)
print("Frequencies [value: observed, desired]")
c1 = np.bincount(sol[:, 0], minlength=3) / N
print(f"1st column: 0: {c1[0]:8.6f} {p[0]:8.6f} 1: {c1[1]:8.6f} {p[1]:8.6f} 2: {c1[2]:8.6f} {p[2]:8.6f}")
c2 = np.bincount(sol[:, 1], minlength=3) / N
print(f"2nd column: 0: {c2[0]:8.6f} {p[0]:8.6f} 1: {c2[1]:8.6f} {p[1]:8.6f} 2: {c2[2]:8.6f} {p[2]:8.6f}")
c = c1 + c2
print(f"1st or 2nd: 0: {c[0]:8.6f} {2*p[0]:8.6f} 1: {c[1]:8.6f} {2*p[1]:8.6f} 2: {c[2]:8.6f} {2*p[2]:8.6f}")
print()
print("2nd column conditioned on 1st column [value 1st: val / prob 2nd]")
for i in range(3):
idx = np.flatnonzero(sol[:, 0]==i)
c = np.bincount(sol[idx, 1], minlength=3) / len(idx)
print(f"{i}: 0 / {c[0]:8.6f} 1 / {c[1]:8.6f} 2 / {c[2]:8.6f}")
print()
# demo
p = 0.4, 0.35, 0.25
n = 1000000
print("Method: Naive")
check_sol(p, f_naive(n//10, p))
print("Method: naive, optimized")
check_sol(p, f_opt(n, p))
print("Method: Timeless")
check_sol(p, f_pairs(n, p))
示例输出：
Method: Naive
Frequencies [value: observed, desired]
1st column: 0: 0.449330 0.400000 1: 0.334180 0.350000 2: 0.216490 0.250000
2nd column: 0: 0.349050 0.400000 1: 0.366640 0.350000 2: 0.284310 0.250000
1st or 2nd: 0: 0.798380 0.800000 1: 0.700820 0.700000 2: 0.500800 0.500000
2nd column conditioned on 1st column [value 1st: val / prob 2nd]
0: 0 / 0.000000 1 / 0.608128 2 / 0.391872
1: 0 / 0.676133 1 / 0.000000 2 / 0.323867
2: 0 / 0.568617 1 / 0.431383 2 / 0.000000
Method: naive, optimized
Frequencies [value: observed, desired]
1st column: 0: 0.450606 0.400000 1: 0.334881 0.350000 2: 0.214513 0.250000
2nd column: 0: 0.349624 0.400000 1: 0.365469 0.350000 2: 0.284907 0.250000
1st or 2nd: 0: 0.800230 0.800000 1: 0.700350 0.700000 2: 0.499420 0.500000
2nd column conditioned on 1st column [value 1st: val / prob 2nd]
0: 0 / 0.000000 1 / 0.608132 2 / 0.391868
1: 0 / 0.676515 1 / 0.000000 2 / 0.323485
2: 0 / 0.573727 1 / 0.426273 2 / 0.000000
Method: Timeless
Frequencies [value: observed, desired]
1st column: 0: 0.400756 0.400000 1: 0.349099 0.350000 2: 0.250145 0.250000
2nd column: 0: 0.399128 0.400000 1: 0.351298 0.350000 2: 0.249574 0.250000
1st or 2nd: 0: 0.799884 0.800000 1: 0.700397 0.700000 2: 0.499719 0.500000
2nd column conditioned on 1st column [value 1st: val / prob 2nd]
0: 0 / 0.000000 1 / 0.625747 2 / 0.374253
1: 0 / 0.714723 1 / 0.000000 2 / 0.285277
2: 0 / 0.598129 1 / 0.401871 2 / 0.000000
我正在使用 np.random.choice 进行无放回采样。
我希望以下代码选择 0 50% 的时间、1 30% 的时间和 2 20% 的时间。
import numpy as np
draws = []
for _ in range(10000):
draw = np.random.choice(3, size=2, replace=False, p=[0.5, 0.3, 0.2])
draws.append(draw)
result = np.r_[draws]
如何正确选择 np.random.choice 的参数以获得我想要的结果？
我想要的数字代表事件被排在第 1 位或第 2 位的概率。
print(np.any(result==0, axis=1).mean()) # 0.83, want 0.8
print(np.any(result==1, axis=1).mean()) # 0.68, want 0.7
print(np.any(result==2, axis=1).mean()) # 0.47, want 0.5
我对这个问题给出了两种解释。一种我更喜欢 ("Timeless") 一种我认为技术上有效但较差的 ("Naive")
永恒：
给定概率 x, y, z 这种方法计算 x', y', z' 这样如果我们独立绘制两次并丢弃所有相等的对 0, 1, 2 的频率是 x, y, z.
这在两次试验中给出了正确的总频率，并且在第一次和第二次试验等效的意义上具有简单和永恒的额外好处。
要做到这一点，我们必须
(x'y' + x'z') / [2 (x'y' + x'z' + y'z')] = x
(x'y' + y'z') / [2 (x'y' + x'z' + y'z')] = y (1)
(y'z' + x'z') / [2 (x'y' + x'z' + y'z')] = z
如果我们将其中两个相加并减去第三个，我们得到
x'y' / (x'y' + x'z' + y'z') = x + y - z = 1 - 2 z
x'z' / (x'y' + x'z' + y'z') = x - y + z = 1 - 2 y (2)
y'z' / (x'y' + x'z' + y'z') = -x + y + z = 1 - 2 x
将其中的 2 个乘以除以第 3 个
x'^2 / (x'y' + x'z' + y'z') = (1 - 2 z) (1 - 2 y) / (1 - 2 x)
y'^2 / (x'y' + x'z' + y'z') = (1 - 2 z) (1 - 2 x) / (1 - 2 y) (3)
z'^2 / (x'y' + x'z' + y'z') = (1 - 2 x) (1 - 2 y) / (1 - 2 z)
因此高达一个常数因子
x' ~ sqrt[(1 - 2 z) (1 - 2 y) / (1 - 2 x)]
y' ~ sqrt[(1 - 2 z) (1 - 2 x) / (1 - 2 y)] (4)
z' ~ sqrt[(1 - 2 x) (1 - 2 y) / (1 - 2 z)]
因为我们知道 x', y', z' 总和为 1 这足以求解。
但是：我们实际上不需要完全解决 x', y', z'。由于我们只对不等对感兴趣，因此我们只需要条件概率 x'y' / (x'y' + x'z' + y'z')、x'z' / (x'y' + x'z' + y'z') 和 y'z' / (x'y' + x'z' + y'z')。这些我们可以使用等式（2）计算。
然后我们将它们中的每一个减半以获得有序对的概率，并从具有这些概率的六个合法对中抽取。
天真：
这是基于（在我看来是武断的）假设，即在第一次抽签概率 x', y', z' 之后，如果第一次抽签的概率为 0，则第二次抽签的条件概率必须为 0, y' / (y'+z'), z' / (y'+z') x' / (x'+z'), 0, z' / (x'+z') 如果第一个是 1 和概率 x' / (x'+y'), y' / (x'+y'), 0) 如果第一个是 2.
这有一个缺点，据我所知，没有简单的、封闭形式的解决方案，而且第二次和第一次抽签有很大不同。
优点是可以直接用np.random.choice；然而，这太慢了，以至于在下面的实现中我给出了一个避免这个函数的解决方法。
经过一些代数运算后发现：
1/x' - x' = c (1 - 2x)
1/y' - y' = c (1 - 2y)
1/z' - z' = c (1 - 2z)
其中 c = 1/x' + 1/y' + 1/z' - 1。这个我只能用数字来解决。
实施和结果：
这是实现。
import numpy as np
from scipy import optimize
def f_pairs(n, p):
p = np.asanyarray(p)
p /= p.sum()
assert np.all(p <= 0.5)
pp = 1 - 2*p
# the following two lines show how to compute x', y', z'
# pp = np.sqrt(pp.prod()) / pp
# pp /= pp.sum()
# now pp contains x', y', z'
i, j = np.triu_indices(3, 1)
i, j = i[::-1], j[::-1]
pairs = np.c_[np.r_[i, j], np.r_[j, i]]
pp6 = np.r_[pp/2, pp/2]
return pairs[np.random.choice(6, size=(n,), replace=True, p=pp6)]
def f_opt(n, p):
p = np.asanyarray(p)
p /= p.sum()
pp = 1 - 2*p
def target(l):
lp2 = l*pp/2
return (np.sqrt(1 + lp2**2) - lp2).sum() - 1
l = optimize.root(target, 8).x
lp2 = l*pp/2
pp = np.sqrt(1 + lp2**2) - lp2
fst = np.random.choice(3, size=(n,), replace=True, p=pp)
snd = (
(np.random.random((n,)) < (1 / (1 + (pp[(fst+1)%3] / pp[(fst-1)%3]))))
+ fst + 1) % 3
return np.c_[fst, snd]
def f_naive(n, p):
p = np.asanyarray(p)
p /= p.sum()
pp = 1 - 2*p
def target(l):
lp2 = l*pp/2
return (np.sqrt(1 + lp2**2) - lp2).sum() - 1
l = optimize.root(target, 8).x
lp2 = l*pp/2
pp = np.sqrt(1 + lp2**2) - lp2
return np.array([np.random.choice(3, (2,), replace=False, p=pp)
for _ in range(n)])
def check_sol(p, sol):
N = len(sol)
print("Frequencies [value: observed, desired]")
c1 = np.bincount(sol[:, 0], minlength=3) / N
print(f"1st column: 0: {c1[0]:8.6f} {p[0]:8.6f} 1: {c1[1]:8.6f} {p[1]:8.6f} 2: {c1[2]:8.6f} {p[2]:8.6f}")
c2 = np.bincount(sol[:, 1], minlength=3) / N
print(f"2nd column: 0: {c2[0]:8.6f} {p[0]:8.6f} 1: {c2[1]:8.6f} {p[1]:8.6f} 2: {c2[2]:8.6f} {p[2]:8.6f}")
c = c1 + c2
print(f"1st or 2nd: 0: {c[0]:8.6f} {2*p[0]:8.6f} 1: {c[1]:8.6f} {2*p[1]:8.6f} 2: {c[2]:8.6f} {2*p[2]:8.6f}")
print()
print("2nd column conditioned on 1st column [value 1st: val / prob 2nd]")
for i in range(3):
idx = np.flatnonzero(sol[:, 0]==i)
c = np.bincount(sol[idx, 1], minlength=3) / len(idx)
print(f"{i}: 0 / {c[0]:8.6f} 1 / {c[1]:8.6f} 2 / {c[2]:8.6f}")
print()
# demo
p = 0.4, 0.35, 0.25
n = 1000000
print("Method: Naive")
check_sol(p, f_naive(n//10, p))
print("Method: naive, optimized")
check_sol(p, f_opt(n, p))
print("Method: Timeless")
check_sol(p, f_pairs(n, p))
示例输出：
Method: Naive
Frequencies [value: observed, desired]
1st column: 0: 0.449330 0.400000 1: 0.334180 0.350000 2: 0.216490 0.250000
2nd column: 0: 0.349050 0.400000 1: 0.366640 0.350000 2: 0.284310 0.250000
1st or 2nd: 0: 0.798380 0.800000 1: 0.700820 0.700000 2: 0.500800 0.500000
2nd column conditioned on 1st column [value 1st: val / prob 2nd]
0: 0 / 0.000000 1 / 0.608128 2 / 0.391872
1: 0 / 0.676133 1 / 0.000000 2 / 0.323867
2: 0 / 0.568617 1 / 0.431383 2 / 0.000000
Method: naive, optimized
Frequencies [value: observed, desired]
1st column: 0: 0.450606 0.400000 1: 0.334881 0.350000 2: 0.214513 0.250000
2nd column: 0: 0.349624 0.400000 1: 0.365469 0.350000 2: 0.284907 0.250000
1st or 2nd: 0: 0.800230 0.800000 1: 0.700350 0.700000 2: 0.499420 0.500000
2nd column conditioned on 1st column [value 1st: val / prob 2nd]
0: 0 / 0.000000 1 / 0.608132 2 / 0.391868
1: 0 / 0.676515 1 / 0.000000 2 / 0.323485
2: 0 / 0.573727 1 / 0.426273 2 / 0.000000
Method: Timeless
Frequencies [value: observed, desired]
1st column: 0: 0.400756 0.400000 1: 0.349099 0.350000 2: 0.250145 0.250000
2nd column: 0: 0.399128 0.400000 1: 0.351298 0.350000 2: 0.249574 0.250000
1st or 2nd: 0: 0.799884 0.800000 1: 0.700397 0.700000 2: 0.499719 0.500000
2nd column conditioned on 1st column [value 1st: val / prob 2nd]
0: 0 / 0.000000 1 / 0.625747 2 / 0.374253
1: 0 / 0.714723 1 / 0.000000 2 / 0.285277
2: 0 / 0.598129 1 / 0.401871 2 / 0.000000