确定性：从均匀到一切
1. sampling 有哪些
信号处理中的采样定理
音乐制作的采样
统计学里大量的非概率性采样：比如分层采样，系统采样
2. 均匀分布
通过时间函数求模就可以获得
如果有分布函数，采样是简单的，因为分布函数实际 是一个“限制”很强的函数，比如 CDF 一定是单调递增， 并且所有概率求和或积分为 1, 利用这两个性质，尤其是 单调递增性质保证了 CDF 的逆的存在，这才能够
随着函数形式的放松，包括：
- blackbox methods
- non-standar distributions
- 基于能量？
连续分布如何采样？
3. Inverse Transform Sampling
Inverse Transform Sampling : Data Science Concepts - YouTube 这个算法就是介绍随机库使用的一种算法，将均匀分布转换成各类其他分布。 算法目的就是给定一个数，将其转成另外一个数 但更全局的性质约束是，给定的数符合 0-1 之间的均匀分布（例如或许可通过计算机计时器求模得到均匀分布），输出的数符合另外一个分布。 即如果有一个均匀分布采样器，uni(), 构造另一个函数 uni2exp，那么 uni2exp(uni()) 能够从指数分布中采样
因此从数学的设计上，inverse_trans 是一个工厂函数，根据需求返回不同的分布采样函数 但是从计算机的实现上，我们通过 random() 函数获得一个数，然后把这个数输入到另外一个函数中， 得到的结果是满足另外一个分布的。注意均匀分布只在 0-1 之间，但映射到别的分布的话，其变量取值可能是无限。 如何从 0-1 到无限呢？关键在于，所有的分布的概率和是 1, 因此任何分布的累积分布函数都在 0 到 1 之间， 这是能够用均匀分布来模拟的核心。
3.1. ritvik 的讲解
听 ritvik 讲完真的觉得一点也不难，基本就是一个数学变换，讲课的优点在于：
- 说清楚这个算法的作用，这是 behind the scene 的原理性部分，我们可能经常用 random 函数，但是不知道他们是如何 sample 出的，不需要完全里理解细节，但作为冰山下 7/8 的支撑
这里用到了 transform 的概念，而不是函数， 线性代数中也叫 transform, 更像是对一个结构进行转换。 例如这里把分布进行转换，以指数分布为例，给定一个均匀分布中的数， T(U)=X 得到一个符合指数分布的数。
如何思考的？ 假设 U 服从均匀分布， X 服从指数分布
- 我们只能关注 U 和 X 的结构，这实际是一个概率的周作业级别的题目(也许是一天的课后题） X 符合指数分布，因此 U~T^-1(X) 符合均匀分布（0-1），但这思路就有问题了。。
假设指数分布的 CDF 是 \(F_X(x)=P(X<=x)\), 此时 X 可以写成 T(U), 因此 CDF 就是 \(P(T(U) \le x)\),
这里的思考实际是模拟性质的， P(X<=x) 表示的是从 X 分布中采样出的结果小于 x 的概率，而采样是通过 T(u) 得到的（例如通过一个复杂的机器得到），现在我们要构造这个 T. P(T(u)<=x) 的概率，继续转换得到 \(P(u<=T^{-1}(x))\) 而由于 u 是 0 到 1 之间的均匀分布，这个概率就等于 \(T^{-1}(x)\) 本身： \(F_X(x)=P(X<=x)=P(u<=T^{-1}(x))=T^{-1}(x)\) 那么 \(T(x)=F^{-1}_X(x)\) 因此给定任何分布的 CDF 函数，只要能求逆，并且逆变换能够保证不等号不改变方向，就能构造出一个 T
3.2. 采样指数分布的例子
例如指数分布的 cdf 是 \(F(x)=1-e^{-\theta x}, x\geq 0\) 其反函数是 \(-\frac{1}{\theta}\log (1-y)\)
于是通过以下就可以从指数分布中采样
import numpy as np
import matplotlib.pyplot as plt
theta=0.5
x = -1/2*np.log(1-np.random.uniform(size=5000))
plt.hist(x, bins=500)
plt.show()
注意 CDF 是描述性质的，它类似判别性质的，给定一个区间，得到概率（一个数，甚至可以用这个概率来分类） 但通过该算法，得到了一个生成性质的函数，给定一个均匀分布的样本，可以构造出另外的样本。
3.3. 采样离散分布的例子 – 幸运大转盘
实际上对于离散的分布，逆变换采样就是更常见一点的模拟大转盘游戏了，假设有一个转盘划分成三个区，每个区大小分别是 1,2,3 给定 random() 函数，模拟大转盘结果，可以在 prefix sum 上做二分搜索来实现：
import random
import bisect
p = [1/6,2/6,3/6]
p_sum = [1/6, 1/2, 1]
x = [bisect.bisect_left(p_sum, random.random()) for x in range(100)]
plt.hist(x, bins=3)
plt.show()
而更加朴素的方法是，通过 0-1 分布生成 1-10 的均匀分布，那么直接转换成 int(random.random() *10) + 1
x = [int(random.random() * 10) + 1 for x in range(10000)]
plt.hist(x)
plt.show()
3.4. 危险的扩充？
这与 VAE 做的事情是一样的，给定一个噪声，最后采样出一个高维的特征
局限是：
- 很多时候就没有 CDF, 例如只有数据，不知道分布，因此要从数据（现象层）里直接采样
- 即便有一个确定的分布（模型层），CDF 求逆也很难，尤其高维度
- 比较适合生成连续概率分布，如果生成离散概率分布，需要做 prefix sum 搜索，参考
- 只能通过连续的0-1均匀分布来生成
如何通过离散均匀分布获得其他离散均匀分布？
4. Accept-reject sampling
4.1. ritvik 解释
也称为 rejection sampling 概率 P(x) =f(x)/NC, 这里 f(x) 一般是知道的，但分母需要归一化，而归一化代价很高，在这种情况下如何进行采样？
比如，假设班上有 1000 个同学，每个人的身高都是已知的，现在刮风，身高越高的人越有概率被风里的叶子击中。 设计一种方式模拟每次被叶子击中的人的身高
这里采用一种很经典的近似的思想，假设 f(x) 知道，但它不是分布，于是我们没法直接采样，那么可以选择一个现有的标准分布 g, 然后从 g 里采样，显然 g 如果和 P(x) 非常接近，那么采样结果直接就是 P(x) 的结果，然而我们不清楚 f(x) 的形状，因此 f(x) 和 g(x) 可能很不一样，这时候对于采样的结果就要进行判断，是接受它作为 P(x) 的结果还是拒绝。
而判断标准是，从从 0-1 中均匀随机一个值，如果这个值小于 f(x)/Mg(x) ，那么就接受，否则拒绝。这里 x 是当前从 g 中采出的样本。
M 是使得 f(x)<Mg(x) 的一个参数，使得 f(x)/Mg(x) 是小于 1 的。 注意实际 M 可以选择足够大，但这会导致采样效率很低，选择与 f 接近 g 是保证采样的高效而不是算法的正确性。
因为可以证明，不论是什么样的 g(x), 这个方法采样出来被接受的数据 D 都是符合分布 P(x) 的，证明手段是 bayes 法则 P(D|A) 表示的是 accept 的数据的分布： \[ P(D|A)=\frac{P(D)P(A|D)}{P(A)} = \frac{g(x)\frac{f(x)}{Mg(x)}}{P(A)} = \frac{\frac{f(x)}{M}}{\int f(x)/M dx} = \frac{f(x)}{\int f(x) dx} = P(x) \]
g 最后都被消除了，因此可以是任意分布。
概念层的理解在于： 用 f(x)/Mg(x) 作为接受的概率，如果这个比值比较大，说明 f(x) 相对 g(x) 较大，意味着在原始 P(x) 中有更高的概率出现，而在近似分布 g 中出现概率低，此时样本 x 显得很珍贵，因为是原始分布里的高概率点，而在 g(x) 里 又不常出现；如果 f(x)/Mg(x) 小，说明在原始分布中不常出现，而在近似分布里概率高，因此不是很重要。
因此拒绝采样的缺点就是对 g 的选择，类似超参数，另外对于现实中高维度的 f(x) 图形会很不均匀，导致 M 很大，采样低效
4.2. 两个离散均匀分布的拒绝采样
一个有趣的例子是，假设给定一个 rand10 表示 1 到 10 之间的均匀分布，如何构造 rand7 ? 或者如何用 randX 构造 randY, 这里 X>Y, 由于二者都是均匀分布，他们的差别仅仅在于范围大小， 因此我们永远可以构造出一个 g(x) 使得 Mg(x) 在 randY 的变量范围内与 f(x) 完全一致，例如 rand10 在 1-7 范围内概率是 1/10, 而 rand7 的概率是 1/7, 因此 M 可以取 10/7 此时在 8-10 的范围内， Mg(x) 恒大于 0, 而 f(X)=0, 因此 f(x)/Mg(x)=0, 也就是始终拒绝。
这里提供的一个思路是，如果有两个均匀分布，第一个范围比第二个大，那么始终拒绝第二个分布里的无效分布即可
在 470. 用 Rand7() 实现 Rand10() - 力扣（Leetcode） 中有一题是要把 rand7 生成 rand10 这里先要用 trick 对 rand7 扩充 这是一个非常好的对拒绝采样理解的题目
4.3. 拒绝采样的局限
由于需要完整的 q(x) 来使得 Mq(x) >= p(x), 因此对于高维度变量（例如图模型），我们也要用一个比较理想的 同维度的高维分布来采样，使得难度大大增加（因为很难找到高维的与原分布接近的分布，导致采样效率非常低，大部分都拒绝），因此拒绝采样也只是用来生成低维度样本，作为 MCMC 的 building block
5. 重要性采样
根据 Importance Sampling - YouTube 介绍，重要性采样是一种蒙特卡洛算法，核心就是当要对 f(x) 求 p(x) 的期望时，如果 p(x) 容易采样（比如是一些标准分布），那么直接用 motecalor 就够了，当如果 p(x) 采样非常困难的时候，就可以引入一个容易采样的分布 q, 只要保证 q(x) 大于 p(x)f(x) ，根据大数定理，采样出的方差就会很小，甚至比采样 p(x) 结果还小，但问题就是如何找到？这是非常困难的，一种纯经验尝试。
\[ \mathbb{E}_{p(x)}[f(x)] \approx \frac{1}{N} \sum_{i=1}^{N} \frac{p(x_i)}{q(x_i)} f(x_i), \quad x_i \sim q(x) \]
ChatGPT背后的数学原理[2]: 重要度采样, 强化学习中的移花接木之术。Importance Sampling_哔哩哔哩_bilibili
p(x)/q(x) 就是样本权重，或者说重要性
这是拒绝采样的扩展，并且，这里的推理用到了与 EM, 变分推断等类似的 trick, 是它们的入门版本 如果我们要推理出某个变量的
6. MCMC sampling 框架
rejection sampling 每次采样是独立的，然而假设采样到一个 x 使得 f(x)/Mg(x) 很大，这说明 x 在原分布概率是比较高的，此时或许可以多从它周围去采样，获得更多概率高的样本。
MCMC 就是利用这种一阶性质来更高效的采样，它目的在于构造一个转移概率，这是一个条件概率 T(x|y) 使得从任意一个 样本出发，遵循 T(x|y) 的转移概率，采样到一定程度后，markov 链就进入均衡状态 detailed balance 也就是之后每次采样得到的分布都是 p(x)
从另一个角度来说，对于满足 p(x) 的任意两个样本 x, y, 进入 detailed balance 后， p(x)T(y|x)=p(y)T(x|y)
MCMC 是一种采样的思路，把采用过程中前后样本的局部性考虑进来，用 Markov 来刻画这种局部性。
拒绝采样和重要性采样都可以作为 MCMC 的一个 building block, 正如 ITS 可以作为生成基本随机函数的手段
7. MCMC 之 Gibbs sampling
Gibbs Sampling : Data Science Concepts - YouTube
背景在于： 高维度分布很难采样，但 P(x_i|x_j) 更简单 用迭代的思路，随机初始化 x0, y0, 根据 P(X|y0) 采样出 x1，然后根据 P(Y|x0) 采样出 y1, 这就得到一个样本点 x1,y1. 依次不断迭代下去。
首先需要能够写出条件概率 P(x|y), 或者通过别的采样方式能够得出条件概率，然后 gibbs 的变化方式：
- 如果有 n 个维度，可以随机选择下一次要变化的维度
- 或者按 block 的方式去采样，例如转移是 p(x1x2|x3)
概念的理解：
- 这是一种 MCMC 采样，本次采样的样本和上一次的有关，并且是一阶的
证明：
- 假设分布是多变量的，两个维度 x 和 y, 那么要证明迭代最够多次后 p(x,y)p(x|y)p(y|x)=p(x,y)
假设当前点是 <x0,y0> 其概率是 p(x0,y0), 那么 <x1, y0> 中 x1 是基于 y0 采样出来的， 因此它遵循概率是 \[ \int p(x_0,y_0)p(x_1|y_0)dx_0 = \int p(y_0)p(x_0|y_0)p(x_1|y_0)dx_0 = p(y_0)p(x_1|y_0) = p(x_1, y_0) \]
积分原因在于 x1 和 x0 是无关的，因此要考虑 x0 的所有范围，第一个等号来自条件概率公式，第二等号来是因为 第二个式子中只有中间项有 x0, 而 p(x0|y0) 对所有 x0 求和就是 1, 于是最后得到下一采样得到的样本的概率 就是均衡后的概率，这符合 detailed balance。
待定： 这里的思路和 EM 非常相似，同样需要证明迭代的结果是符合原分布的，同样也与优化 SVM 的 SMO 算法一样。 甚至这种迭代和随机游走也是类似的 以上EM ，SMO 我自己也不是很清楚，不能妄下结论。
我可以很快写出这个算法的伪代码
def gibbs(dist, k):
p_x_y = get_cond_dist(dist, y)
p_y_x = get_cond_dist(dist, x)
x0, y0 = random(2)
samples = [(x0, y0)]
for i in range(k):
x0, x1 = samples[-1]
x1 = sample(p_x_y(y0))
y1 = sample(p_y_x(x1))
samples.append((x1, y1))
return samples
具体的实现见 sampling
gibbs 的局限：
- 如果 p(x|y) 会使得某些 x 概率为 0 或 p(x|y) 是一个很聚集的概率分布, 那么这些 x 值不能马上采样到，需要很多步才能跳转到（如图中右下角 issue 中对角矩阵的形式）
8. MCMC 之 Metropolis-Hasting
Metropolis - Hastings : Data Science Concepts - YouTube 以上视频讲的一如既往的好，控制了概念的复杂度，主要讲解 Metropolis 部分，讲清楚了缘由和使用方式
Coding MCMC : Data Science Code - YouTube 这个视频里对比了 rejections 采样不同参数的效果和 Metroplis sampling 的差异，代码见 YouTubeVideoCode/MCMC Experiments.ipynb at main · ritvikmath/YouTubeVideoCode
这里有一个类比，从 Gibbs sampling 到 Metropolis sampling 类比于从 Inverse Transform 到 Accept-reject
原因在于 Gibbs 中我们需要完整的分布，例如至少要知道条件分布 p(x|y), Inverse Transform 则需要知道目标分布 p(x) 而 Accept-reject 只需要知道分布的 logits 函数 f(x), 然后通过在另外一个简单的分布 g 里采样，依据 f(x)/Mg(x) 的比值来做接受或拒绝。
Metropolis 是在 Markov 的假设下来实现类似的效果，只知道 f(x) 以及上个采样的值 a, 根据这个采样 a 来构造 一个分布如 g(a), 比如常见的就是以 a 为均值的正态分布, 从 g 中采样出新的样本 b, 如果 f(b)/f(a) 大于 1, 意味着我们采样到一个比之前概率更高的样本，这是比较珍贵的，因为毕竟概率越高采样的点应该越多，那么保留这个样本，如果小于 1, 意味着找到一个比之前概率更低的样本，那么以 f(b)/f(a) 的概率来接收它。 这是直觉上的理解
而证明也是直接的，根据 markov 性质，要证明给定一个点 a, 根据采样策略得到另一个点后，这个点的分布还是满足 p 关键就在于根据 a 写出转移后的结果 b 的概率 b ~ p(a)g(b;a)A(a->b) a ~ p(b)g(a;b)A(b->a) 要使得以上两个分布相等： A(a->b)/A(b->a)=(p(b)g(b;a))/(p(a)g(a;b))=(f(b)g(b;a))/(f(a)g(a;b)) 设 pq = (f(b)g(b;a))/(f(a)g(a;b)) 那么 pq > 1 时，可以令 A(a->b) =1, 直接接收，A(b->a) 就是 1/pq 但 pq < 1 时，可以令 A(a->b) = pq, 以 pq 概率接收 b, A(b->a) 就是 1 这里可以看到 pq 非常灵活
因此始终保证 detailed balance 成立
Metropolis-Hasting 实际包括两类：
- Metropolis sampling: 选择的分布 g(a)是对称的，这时候 g(b)(a) = g(a)(b), 所以只要考虑 f(a)/f(b)
- Metropolis-Hasting sampling: 选择的分布不对称时
9. Tompson sampling
YouTubeVideoCode/Thompson Sampling.ipynb at main · ritvikmath/YouTubeVideoCode 这里是用 Bayes 的思路处理 multi-bandit 问题，takeaway 的点是，可以用方差很大的高斯分布来模拟均匀分布， 以此作为先验概率， 但这个例子非常理想，是一个 toy model，即方差固定 1 ，需要估计的是三个饭店好评的均值，假设他们都是高斯分布， 然后开始时几乎随机选择一个，然后根据这 1 个样本用 Bayes 去更新当前饭店的好评的分布，然后再从三个分布里继续采样， 选择得分最高的去访问，依次迭代，思路是很简单的
multi-bandit 不是一个混合概率，我们知道有 n 个概率分布，每次可以选择从哪个取出，它的目的是
10. 要多大的样本集
YouTubeVideoCode/Sample Size.ipynb at main · ritvikmath/YouTubeVideoCode 中展现了一个这样的例子，假设公司的产品 A 已经卖了很多年，评分是 4 现在要对该产品做一点改进，变成产品 B，为了检查改进带来的评分变化，需要先找一部分人来体验这个产品并且打分， 现在问题是，要找多少人来测评才能比较反映出新产品真实的评价
显然，理想状态是把产品 B 给用过产品 A 的人所有用一遍，但这就直接产品了。
但我们不需要真正分析产品 B, 而是研究采样大小对结果的影响。 于是手动构造一个分布，这个分布的均值 u 和方差 s 都是研究对象，因此当前关注的变量是 u, s 和 n
为了简化，可以只关注 u 和 n. 如果从理论研究，那么采样出越多 n, 统计均值就趋近于真实均值，导向大数定律的证明。当这里并不需要这种证明 从分布中采样出 n 个样本，统计其均值并且 plot 出图进行分析。 最后的结论显然符合直觉：
- 如果新产品的真实的 u 与原产品 A 的评分很接近，那么需要更多的样本才能确定产品 B 相对 A 是好还是坏
这个问题更多反映一种元思路：
- 因为我们分析的是决定一次采样的超参数 n，所以又要得到更多不同 n 的采样
11. 随机的开始
给定数字 2 ，你可以说它是整数，是自然数，是偶数，但没法判断它是否是随机数。
随机不是针对一个独立的数而言的，而是数之间的关系，因此我们更多谈论随机序列
如果一个数和跟在它后面的那个数之间的真实关系没有明显的现实意义（无客观意义），它就是
最初的随机数来自于物理机械装置，比如掷色子，计算机发明之后，可以把掷色子的结果搜集成表格做成要给随机数表保存在内存里，而如果要让机器能够直接动态获得随机数，可以依赖于电器装置，比如电流噪声，
1946 年前后，冯·诺伊曼开始用普通算术操作来产生随机数，比如平方取中法：
11.1. 看似随机即可
def sqrt_mid_random(seed, length):
def _gen():
nonlocal seed
while 1:
seed_str = str(seed).zfill(length*2)
start = (length + 1) // 2
print(seed_str)
seed = int(str(int(seed_str)**2)[start: start + length])
yield seed
return _gen
gen = sqrt_mid_random(33, 4)()
for i in range(10):
print(next(gen))
00000033 89 00000089 21 00000021 1 00000001
ValueErrorTraceback (most recent call last)
<ipython-input-17-6bc17fafb6e2> in <module>
13
14 for i in range(10):
---> 15 print(next(gen))
<ipython-input-17-6bc17fafb6e2> in _gen()
6 start = (length + 1) // 2
7 print(seed_str)
----> 8 seed = int(str(int(seed_str)**2)[start: start + length])
9 yield seed
10 return _gen
ValueError: invalid literal for int() with base 10: ''
N. Metropolis 对平方取中的随机性进行了广泛分析
11.2. 随意设计一个算法获得随机的概率？
Kunth 原本觉得设计一个足够复杂的算术函数来产生伪随机数序列并不是很难的事情，但实际上，大多数设计出来的随机数 生成器都不太随机，比如冯·诺伊曼的平方取中法：
这使得伪随机的产生也需要一套理论：
- 如何判断随机程度
- 哪种方法能够更随机
如何验证序列的随机度
习题 6,7,9,10 验证随机性
计算机程序设计艺术太底层了，应该梳理完上层应用再回到随机本质