python按指定概率抽样_python – 没有替换概率的抽样

我对这个问题有两种解释.一个我更喜欢(“永恒”)和一个我认为技术上有效但劣等(“天真”)

永恒：

给定概率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个除以第三个

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,y’/(y’z’),z’/( y’z’)如果第一个是0 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