与值和位置无关的位转换
可以在多种情况下计算某些位状态转换到另一个位状态的概率。
最简单的一种是当某个位转换到不同状态的给定概率p 时，这与位值、位在位状态中的位置以及其他位转换无关.
当然，位不翻转的概率由q == 1 - p给出。
（n有两个结果的独立事件的统计数据为studied extensively。）
对于更多位，可以通过乘法组合多个位转换的概率。
从a 到b 的转换概率（其中a 和b 是相同长度的两个位配置n）取决于位转换t_ab 和非转换的数量s_ab == n - t_ab:
p(a, b) == (p ** t_ab) * (q ** s_ab)
例如，转换：0b00011 到 0b00101 由以下公式给出：
p(0b00011, 0b00101) == (q ** 3) * (p ** 2)
请注意，这与例如0b011 到 0b101 的转换概率，因为要考虑的位数起作用。
给定一个计算数字中1个数的函数：
def count_set_bits(num):
result = 0
while num:
result += num & 1
num >>= 1
return result
计算t 的一种简单方法是通过xor 运算符：
t = count_set_bits(a ^ b)
因此，可以通过简单的循环“手动”计算转移概率矩阵w_bits。
除非加速显式循环，否则计算速度非常慢。
此用例最简单的加速之一是使用Numba。
所有_nb-ending 函数都使用它加速。
可以将 fastmath 标志 nb.njit(fastmath=True) 设置为可能将执行时间减少几个百分点。
import numpy as np
import numba as nb
@nb.njit
def count_set_bits(num):
result = 0
while num:
result += num & 1
num >>= 1
return result
@nb.njit
def w_bits_sym_cb_nb(n, p=0.2):
if n > 0:
q = 1 - p
m = 2 ** n
result = np.empty((m, m), dtype=np.float_)
for i in range(m):
for j in range(i + 1):
t = count_set_bits_nb(i ^ j)
s = n - t
result[i, j] = result[j, i] = (p ** t) * (q ** s)
return result
else:
return np.empty((0, 0))
（请注意，count_set_bits() 也已加速）。
或者，可以通过重复 1 位情况的基本概率矩阵来构造逐元素乘法概率矩阵：
0 1
0 q p
1 p q
具有两次重复的幂，例如两个字节：
q p q p q q p p
p q p q X q q p p
q p q p p p q q
p q p q p p q q
这可以再次通过“手动”循环计算：
@nb.njit
def w_bits_sym_lm_nb(n, p=0.2):
if n > 0:
b = 2
m = b ** n
q = 1 - p
base = np.array([[q, p], [p, q]])
result = np.ones((m, m), dtype=base.dtype)
for k in range(n):
bk = (b ** k)
for i in range(m):
for j in range(m):
result[i, j] *= base[i // bk % b, j // bk % b]
return result
else:
return np.empty((0, 0))
但是，一种更快的方法是使用广播乘法执行与重复元素的逐元素矩阵乘法（@PierreD's answer 的完善版本）：
import numpy as np
def bc_mul(a, b):
nm = len(a) * len(b)
return (a[:, None, :, None] * b[None, :, None, :]).reshape(nm, nm)
def w_bits_sym_bm(n, p=0.2):
if n > 0:
base = np.array([[1 - p, p], [p, 1 - p]])
result = base.copy()
for i in range(1, n):
result = bc_mul(base, result)
return result
else:
return np.empty((0, 0))
请注意，由于bc_mul() 是关联的，因此可以将循环内的行写为result = bc_mul(base, result) 或result = bc_mul(result, base)，但性能却截然不同！
最后一种方法也相当快，尤其是对于较大的n 渐近，主要是因为它执行的乘法次数呈指数级减少。
同样可以用 Numba 重写，但性能相似（但性能稍慢）：
@nb.njit
def bc_mul_nb(a, b):
n = len(a)
m = len(b)
nm = n * m
result = np.empty((nm, nm), dtype=a.dtype)
for i in range(n):
for j in range(m):
for k in range(n):
for l in range(m):
result[i * m + j, k * m + l] = a[i, k] * b[j, l]
return result
@nb.njit
def w_bits_sym_bm_nb(n, p=0.2):
if n > 0:
base = np.array([[1 - p, p], [p, 1 - p]])
result = base.copy()
for i in range(1, n):
result = bc_mul_nb(base, result)
return result
else:
return np.empty((0, 0))
更多关于执行速度（包括基准）的信息如下。
值相关/位置无关位转换
一个稍微复杂，更有趣的场景场景是当 0 到 1 和 1 到 0 的概率不同，但仍然独立于位置等时。
两者都可以从base概率矩阵计算：
0 1
0 p00 p01
1 p10 p11
其中p00、p01、p10 和p11 是一位从一种状态转换到另一种状态的概率。
当然：
p00 == 1 - p01p11 == 1 - p10
和以前一样，对于更多位，可以通过乘法组合多个位转换的概率。
这本质上是上述的不对称版本。
从a 到b 的转换概率（其中a 和b 是相同长度的两位配置）取决于转换的数量t00_ab、t01_ab、t10_ab , t11_ab 乘以它们各自的概率（对称情况下使用的符号，t01 和 t10 对应于 t，t00 和 t11 对应于 s）：
p(a, b) == (
(p00 ** t00_ab) *
(p01 ** t01_ab) *
(p10 ** t10_ab) *
(p11 ** t11_ab))
例如，转换：0b00011 到 0b00101 由以下公式给出：
p(0b00011, 0b00101) == (p00 ** 2) * (p01 ** 1) * (p10 ** 1) * (p11 ** 1)
当然，所有这些都可以与上述类似的计算。
设置位计数方法可以直接在~a & b 和a & ~b 上与a & b 一起使用来计数位转换：
@nb.njit
def w_bits_cb_nb(n, p01=0.2, p10=-1):
if n > 0:
p10 = p10 if p10 >= 0 else p01
p00 = 1 - p01
p11 = 1 - p10
m = 2 ** n
result = np.empty((m, m), dtype=np.float_)
for i in range(m):
for j in range(m):
t11 = count_set_bits_nb(i & j)
t01 = count_set_bits_nb(~i & j)
t10 = count_set_bits_nb(i & ~j)
t00 = n - (t11 + t01 + t10)
result[i, j] = \
(p00 ** t00) * (p11 ** t11) * (p01 ** t01) * (p10 ** t10)
return result
else:
return np.empty((0, 0))
或者可以在单个循环中稍微更有效地完成（与@Viglione's current answer 中的类似但更快）：
@nb.njit
def bit_diff_nb(a, b, n):
t11 = t01 = t10 = 0
t00 = n
while a | b:
aa = a & 1
bb = b & 1
t11 += aa & bb
t01 += ~aa & bb
t10 += aa & ~bb
a >>= 1
b >>= 1
t00 = n - (t11 + t01 + t10)
return t00, t11, t01, t10
@nb.njit
def w_bits_bd_nb(n, p01=0.2, p10=-1):
if n > 0:
p10 = p10 if p10 >= 0 else p01
p00 = 1 - p01
p11 = 1 - p10
m = 2 ** n
result = np.empty((m, m), dtype=np.float_)
for i in range(m):
for j in range(m):
t00, t11, t01, t10 = bit_diff_nb(i, j, n)
result[i, j] = \
(p00 ** t00) * (p11 ** t11) * (p01 ** t01) * (p10 ** t10)
return result
else:
return np.empty((0, 0))
另外，所有其他方法都可以轻松扩展到这种情况：
@nb.njit
def w_bits_lm_nb(n, p01=0.2, p10=-1):
if n > 0:
p10 = p10 if p10 >= 0 else p01
b = 2
m = b ** n
base = np.array([[1 - p01, p01], [p10, 1 - p10]])
result = np.ones((m, m), dtype=base.dtype)
for k in range(n):
bk = (b ** k)
for i in range(m):
for j in range(m):
result[i, j] *= base[i // bk % b, j // bk % b]
return result
else:
return np.empty((0, 0))
def w_bits_bm(n, p01=0.1, p10=-1):
if n > 0:
p10 = p10 if p10 >= 0.0 else p01
base = np.array([[1 - p01, p01], [p10, 1 - p10]])
result = base.copy()
for i in range(1, n):
result = bc_mul(base, result)
return result
else:
return np.empty((0, 0))
def w_bits_bmi(n, p01=0.1, p10=-1):
if n > 0:
p10 = p10 if p10 >= 0.0 else p01
base = np.array([[1 - p01, p01], [p10, 1 - p10]])
result = base.copy()
for i in range(1, n):
result = bc_mul(result, base)
return result
else:
return np.empty((0, 0))
结果一致性
为了完整起见，我还包含了currently accepted and top voted answer 方法（类似于w_bits_bd_nb()，但使用二进制字符串且没有加速）和一些桥接代码来获取底层 NumPy 数组：
import pandas as pd
def calc_p(sent, observed, p01, p10):
p00 = 1 - p01
p11 = 1 - p10
r = 1
for i, _ in enumerate(sent):
if sent[i] == "0":
r *= p00 if observed[i] == "0" else p01
else:
r *= p10 if observed[i] == "0" else p11
return r
def generate_error_matrix(n, p01, p10):
labels = [f"{i:0{n}b}" for i in range(0, 2 ** n)]
result = pd.DataFrame(index=labels, columns=labels)
for rowIndex, row in result.iterrows():
for columnIndex, _ in row.items():
result.at[rowIndex, columnIndex] = calc_p(rowIndex, columnIndex, p01, p10)
return result
def w_bits_bs_pd(n, p01=0.2, p10=-1):
p10 = p10 if p10 >= 0.0 else p01
return generate_error_matrix(n, p01, p10).to_numpy().astype(float)
funcs = (
w_bits_bm, w_bits_bmi,
w_bits_cb_nb, w_bits_bd_nb, w_bits_lm_nb,
w_bits_bm_nb, w_bits_bmi_nb,
w_bits_sym_cb_nb, w_bits_sym_bm_nb, w_bits_sym_lm_nb,
w_bits_bs_pd)
n = 2
base = funcs[0](n)
print(f"{'ProbRowsSumTo1:':>27} {np.allclose(np.sum(base, 0), np.ones(2 ** n))}")
x = w_bits_bm(10, 0.2, 0.2)
print(f"{'(p01 == p10) -> Symmetric:':>27} {np.allclose(x, x.T)}")
x = w_bits_bm(10, 0.2, 0.4)
print(f"{'(p01 != p10) -> Asymmetric:':>27} {not np.allclose(x, x.T)}")
print()
for func in funcs:
res = func(n)
print(f"{func.__name__!s:>20} Same: {np.allclose(base, res)}")
print(func(2))
print()
ProbRowsSumTo1: True
(p01 == p10) -> Symmetric: True
(p01 != p10) -> Asymmetric: True
w_bits_bm Same: True
[[0.64 0.16 0.16 0.04]
[0.16 0.64 0.04 0.16]
[0.16 0.04 0.64 0.16]
[0.04 0.16 0.16 0.64]]
w_bits_bmi Same: True
[[0.64 0.16 0.16 0.04]
[0.16 0.64 0.04 0.16]
[0.16 0.04 0.64 0.16]
[0.04 0.16 0.16 0.64]]
...
下面的代码表明：
- 所有函数都给出相同的结果
- 如果
p01 == p10 转移矩阵是对称的
- 如果
p01 != p10 转移矩阵是不对称的
- 所有行加起来为一（单独）
基准测试
由于大多数对称实现与非对称实现非常相似，因此它们已从基准测试中省略。
funcs = (
w_bits_bm, w_bits_bmi,
w_bits_cb_nb, w_bits_bd_nb, w_bits_lm_nb,
w_bits_bm_nb, w_bits_bmi_nb,
w_bits_sym_cb_nb, w_bits_bs_pd)
timings = {}
for n in range(1, 12):
print(f"n = {n}")
timings[n] = []
base = funcs[0](n)
for func in funcs:
res = func(n)
timed = %timeit -r 4 -n 8 -q -o func(n)
timing = timed.best * 1e6
timings[n].append(timing)
print(f"{func.__name__:>24} {np.allclose(base, res)} {timing:10.3f} µs")
要绘制：
import pandas as pd
df = pd.DataFrame(data=timings, index=[func.__name__ for func in funcs]).transpose()
df.plot(marker='o', logy=True, xlabel='Num. bits n / #', ylabel='Best timing / µs')
制作：
这确实表明基于广播乘法的解决方案对于较大的n 而言是渐进的，性能最高，但总体上在所有尺度上都相当出色。
请注意，由于计算复杂度呈指数增长，因此时序已按 y 对数标度绘制。
另请注意，w_bits_bs_pd() 比其他的要慢几个数量级。
更好的输出
像往常一样，在处理表格/矩阵等众所周知的对象时，使用特定的工具会很有好处。
如果想要获得漂亮的输出，可以使用Pandas（类似于@Viglione's answer 中所做的）和Seaborn 以获得更好的可视化效果：
import pandas as pd
import seaborn as sns
def gen_bit_transitions(n, p01=0.2, p10=-1, func=w_bits_bm):
data = func(n, p01, p10)
labels = [f"{i:0{n}b}" for i in range(2**n)]
return pd.DataFrame(data, index=labels, columns=labels)
df = gen_bit_transitions(3, 0.4, 0.2)
sns.set(rc={'figure.figsize': (8, 7)})
sns.heatmap(df, annot=True, vmin=0.0, vmax=1.0)
df = gen_bit_transitions(5, 0.4, 0.2)
sns.set(rc={'figure.figsize': (9, 8)})
sns.heatmap(df, annot=False, vmin=0.0, vmax=1.0)