详细介绍：【概率论】均匀分布的伪随机数

文章目录

• 伪随机数的基本概念
• 线性同余法
• Mersenne 旋转算法
• Xoshiro256++
• PCG64

伪随机数的基本概念

伪随机数（Pseudo-random Number）是通过确定性算法生成的序列，看似随机但可重现。其核心特点是：

• 确定性：相同种子产生相同序列。
• 统计特性：需满足均匀性、独立性等统计检验。
• 周期性：序列最终会重复（周期需足够长）。

与真随机数（依赖物理熵源如热噪声、放射性衰变）不同，伪随机数因高效可控，广泛用于仿真、蒙特卡罗方法等过程辨识场景。

线性同余法

在过程辨识课程中介绍机器生成随机数（特别是伪随机数）时，需从理论基础、生成方法、算法实现到应用场景进行系统阐述。以下为分步骤的讲解框架，重点聚焦于生成**[0,1]均匀分布随机数**的具体方法及其原理。 线性同余生成器： linear congruential generator（LCG）

原理：通过线性递推公式生成整数序列，再映射到[0,1]区间。

• 递推公式：
  $$X_{n+1}=(a{X}_n+c)\mathrm{mod}m$$
  • (X_n)：当前状态（种子）
  • (a)：乘数（如1103515245）
  • (c)：增量（如12345）
  • (m)：模数（取(2^{31})等大整数）。

转换到[0,1]：
$$U_n=\frac{X_n}{m}$$
示例代码（Python）：

def lcg(seed, a=1664525, c=1013904223, m=2**32):
while True:
seed = (a * seed + c) % m
yield seed / m # 映射到[0,1]

直方图测试

import numpy as np
import matplotlib.pyplot as plt
seed=0
iters=1000;
v=np.zeros(iters)
for k in range(iters):
seed= lcg(seed=seed)
v[k]=seed/2**32
plt.hist(v)
plt.show()

直方图测试

优缺点：

• ✅ 实现简单、计算快。
• ❌ 周期短、低维均匀性差（需谨慎选择参数）。

Mersenne 旋转算法

梅森旋转算法（Mersenne Twister，MT. 以MT19937著称）是一种基于线性反馈移位寄存器的伪随机数生成算法，由松本真和西村拓士于1997年提出。其核心优势在于极长周期（2¹⁹⁹³⁷ -1）、高维均匀分布性（623维内均匀），以及高效性。以下是其具体算法步骤及原理详解，结合MT19937（32位版本）为例：

1. 算法基础与核心参数

MT19937-32的关键参数如下表所示：
参数 符号 取值（16进制） 作用
状态数组长度 n 624 存储内部状态
偏移量 m 397 状态混合偏移
位掩码 r 31 (0x7FFFFFFF) 分割高低位
旋转矩阵 a 0x9908B0DF 非线性变换
初始化乘数 f 1812433253 种子扩展
位移参数 u, s, t, l 11, 7, 15, 18 输出调和
掩码参数 b, c, d 0x9D2C5680, 0xEFC60000, 0xFFFFFFFF 输出调和

注：d 通常取全1掩码（0xFFFFFFFF），在32位运算中可省略。

2. 算法详细步骤

(1) 初始化阶段：构建初始状态数组

输入一个32位种子 seed，生成长度为 n=624 的状态数组 MT[0…623]：

def initialize(seed):
MT = [0] * 624
MT[0] = seed
for i in range(1, 624):
# 高位保留30位，低位2位参与运算 
prev = MT[i-1]
temp = f * (prev ^ (prev >>
30)) + i
MT[i] = temp &
0xFFFFFFFF # 取低32位 

关键操作：
• prev ^ (prev >> 30)：保留高位随机性，消除低位规律性。

• 乘数 f=1812433253 扩散随机性，加法 +i 避免重复序列。

(2) 扭转阶段（Twist）：更新状态数组

当所有状态值使用后（即生成624个随机数），需重构整个状态数组：

def twist():
for i in range(624):
# 组合MT[i]的高1位和MT[i+1]的低31位 
x = (MT[i] &
0x80000000) | (MT[(i+1) % 624] &
0x7FFFFFFF)
# 右移1位并可能异或a 
xA = x >>
1
if x % 2 != 0: # 若x为奇数 
xA ^= a # 异或旋转矩阵a 
# 更新状态：异或偏移m处的状态 
MT[i] = MT[(i + m) % 624] ^ xA

核心目的：
• 通过非线性移位（异或 a）和状态混合（异或 MT[i+m]）打破线性相关性。

(3) 输出处理阶段（Tempering）：生成最终随机数

从当前状态 MT[index] 提取值，经4步位操作调和：

def extract_number():
y = MT[index]
y ^= (y >> u) # 消除高位局部相关性 
y ^= (y << s) & b # 消除低位局部相关性 
y ^= (y << t) & c
y ^= (y >> l)
return y &
0xFFFFFFFF # 输出32位整数 

调和步骤的作用：
• 通过位移和掩码（b, c）将状态位的线性关系转化为混沌输出，增强统计均匀性。

3. 代码实现（Python精简版）

class MersenneTwister
:
def __init__(self, seed):
self.n, self.m = 624, 397
self.MT = [0] * self.n
self.index = 0
self.init_genrand(seed)
def init_genrand(self, seed):
self.MT[0] = seed
for i in range(1, self.n):
temp = 1812433253 * (self.MT[i-1] ^ (self.MT[i-1] >>
30)) + i
self.MT[i] = temp &
0xFFFFFFFF
def twist(self):
for i in range(self.n):
x = (self.MT[i] &
0x80000000) | (self.MT[(i+1) % self.n] &
0x7FFFFFFF)
xA = x >>
1
if x % 2 != 0:
xA ^= 0x9908B0DF
self.MT[i] = self.MT[(i + self.m) % self.n] ^ xA
self.index = 0
def genrand_int32(self):
if self.index >= self.n:
self.twist()
y = self.MT[self.index]
y ^= y >>
11
y ^= (y <<
7) &
0x9D2C5680
y ^= (y <<
15) &
0xEFC60000
y ^= y >>
18
self.index += 1
return y &
0xFFFFFFFF

直方图测试

import numpy as np
import matplotlib.pyplot as plt
iters=1000
v=np.zeros(iters)
N=0
for k in range(iters):
mt = MersenneTwister(N)
rand_uniform = mt.genrand_int32() / (2**32) # 转换为[0,1)区间 
v[k]=rand_uniform
N=mt.genrand_int32()
plt.hist(v)
plt.show()

直方图测试

4. 性能优化与变体

5. 空间优化：
   • TinyMT：状态数组仅127位，适用于嵌入式系统。

6. 速度优化：
   • SFMT（SIMD-Friendly MT）：利用CPU向量指令并行生成随机数，提速200%。

7. 安全增强：
   • CryptMT：通过过滤函数输出，满足密码学安全要求。

8. 算法验证与应用场景

• 统计检验：

通过Diehard、TestU01等测试验证623维内均匀性。
• 典型应用：

• 蒙特卡洛模拟（如金融风险评估）

• 游戏引擎（如地图生成、NPC行为）

• 机器学习（如参数初始化、数据增强）

关键注意事项

• 种子选择：避免全0种子（导致状态全0失效），推荐使用时间戳。

• 非密码安全：攻击者可反向推导状态（需340个连续输出），加密场景需换用安全算法。

梅森旋转算法以其平衡的性能与随机性质量，成为Python（random模块）、NumPy等库的默认生成器。理解其状态更新机制（Twist）和输出调和（Tempering）是掌握现代PRNG设计的关键。

Xoshiro256++

Xoshiro256++ 是由 David Blackman 和 Sebastiano Vigna 设计的伪随机数生成器（PRNG），属于 xoshiro（xor/shift/rotate）系列算法。其核心设计目标是：

• 高速度：利用位运算（移位、异或、旋转）实现极低开销的随机数生成。
• 长周期：周期为 (2^{256}-1)，远高于传统算法（如梅森旋转算法的 (2^{19937}-1)）。
• 统计质量优异：通过严苛的随机性测试（如 TestU01 的 BigCrush 套件）。

1. 状态维护
算法维护一个包含 4 个 64 位整数的状态数组 S = [s0, s1, s2, s3]，初始状态需非全零。

2. 状态更新步骤
每生成一个随机数，状态按以下步骤更新：

def next_state(s0, s1, s2, s3):
t = (s0 + s3) &
0xFFFFFFFFFFFFFFFF # 64位模加
s0_new = (s0 ^ (s1 <<
17)) &
0xFFFFFFFFFFFFFFFF
s1_new = s1 ^ s2
s2_new = s2 ^ s3
s3_new = s3 ^ ((s3 <<
45) | (s3 >>
19)) # 循环左移45位
return s0_new, s1_new, s2_new, s3_new, t

关键操作说明：

• 异或移位（s0 ^ (s1 << 17)）：破坏线性相关性，增强混沌性。
• 循环移位（(s3 << 45) | (s3 >> 19)）：避免状态退化，维持长周期。
• 模加法（t = s0 + s3）：作为输出前的非线性变换，改善低维均匀性。

3. 输出生成
最终输出的随机数为 t（64 位整数），可通过除法映射到 [0, 1) 区间：

output = t / (2**64) # 转换为[0,1)浮点数

三、Python 完整实现

class Xoshiro256PP
:
def __init__(self, seed=0):
# 初始化状态：使用SplitMix64生成初始状态（确保非零）
self.state = [0] * 4
temp_seed = seed
for i in range(4):
temp_seed = self._splitmix64(temp_seed)
self.state[i] = temp_seed
def _splitmix64(self, seed):
# SplitMix64 用于种子扩展
seed = (seed + 0x9E3779B97F4A7C15) &
0xFFFFFFFFFFFFFFFF
z = seed
z = (z ^ (z >>
30)) * 0xBF58476D1CE4E5B9 &
0xFFFFFFFFFFFFFFFF
z = (z ^ (z >>
27)) * 0x94D049BB133111EB &
0xFFFFFFFFFFFFFFFF
return z ^ (z >>
31)
def next(self):
s0, s1, s2, s3 = self.state
# 计算临时输出值（非线性变换）
t = (s0 + s3) &
0xFFFFFFFFFFFFFFFF
# 更新状态
s0_new = (s0 ^ (s1 <<
17)) &
0xFFFFFFFFFFFFFFFF
s1_new = s1 ^ s2
s2_new = s2 ^ s3
s3_new = (s3 ^ ((s3 <<
45) | (s3 >>
19))) &
0xFFFFFFFFFFFFFFFF # 循环左移45位
self.state = [s0_new, s1_new, s2_new, s3_new]
return t
def random(self):
# 生成[0,1)区间的浮点数
return self.next() / (2**64)

直方图测试

import numpy as np
import matplotlib.pyplot as plt
iters=1000
v=np.zeros(iters)
N=0
for k in range(iters):
rng = Xoshiro256PP(seed=k)
rand_uniform = rng.random() # 转换为[0,1)区间 
v[k]=rand_uniform
plt.hist(v)
plt.show()

直方图测试

关键代码解析：

1. _splitmix64：用于将种子扩展为 4 个 64 位状态值，避免初始状态弱熵问题。
2. 循环移位：(s3 << 45) | (s3 >> 19) 等效于循环左移 45 位（64 位环境下）。
3. 位掩码：& 0xFFFFFFFFFFFFFFFF 确保结果在 64 位范围内（Python 无整数溢出需手动处理）。

四、算法特性与对比

特性	Xoshiro256++	梅森旋转（MT19937）
周期	(2^{256}-1)	(2^{19937}-1)
状态空间	32 字节	2.5 KB
速度	⭐⭐⭐⭐ (比 MT 快 3-5 倍)	⭐⭐
统计质量	通过 TestU01 BigCrush	通过，但高维有缺陷
适用场景	游戏、模拟、轻量加密	通用场景

注：Xoshiro256++ 不适用于密码学（非加密安全），但适合高性能模拟。

五、应用场景建议

1. 蒙特卡洛模拟：高吞吐量需求下的随机数生成。
2. 游戏开发：NPC 行为、地图生成等高频随机事件。
3. 机器学习：数据增强、参数初始化（需结合种子管理）。

六、扩展优化方向

1. 并行化：每个线程独立维护状态数组，避免锁竞争。
2. 跳跃函数：预计算 (2^{128}) 步状态，支持分布式 RNG 同步。

如需加密安全方案，可替换为 SHA-256 或 ChaCha20 等算法。

PCG64

PCG64（Permuted Congruential Generator）是一种高性能的伪随机数生成器（PRNG），由Melissa O’Neill于2014年提出。它结合了线性同余生成器（LCG）的高效性和输出变换函数的统计优化，适用于科学计算、模拟和高并发场景。以下从算法原理、关键特性及Python实现三方面展开详解。

一、算法原理
PCG64的核心结构分为两部分：状态更新（LCG）和输出变换（非线性位操作）。
1. 状态更新（LCG）
采用128位线性同余生成器：
$$s_{n+1}=(a\times s_n+c)\mathrm{mod}{2}^{128}$$

• $s_n$：当前状态（128位整数）
• $a$$：乘数（固定常数，需满足 (a \equiv 5 \mod 8) 以保证周期）
• $c$：增量（奇数，确保满周期 (2^{128}))

2. 输出变换
通过位操作打破LCG的序列相关性，生成64位随机数：

u = s_n >>
(64 - p) # 取高64位，p为位移参数（通常为32）
output = u ^ (u >> p) # 右移后异或，增强随机性

• 位移参数：p 控制混淆强度（默认32）
• 位旋转：部分变体引入循环移位（如PCG64DXSM），进一步改善统计质量。

---

⚡ 二、关键特性
1. 统计质量

• 通过 TestU01的BigCrush测试，高维均匀性优于梅森旋转（MT19937）。
• 改进版PCG64DXSM：解决原始版本在大规模并行流（>百万级）中的统计弱点，采用DXSM输出函数（xorshift-multiply哈希构造），雪崩效应更优。

2. 性能与资源

• 速度：比MT19937快2-3倍（避免大状态数组操作）。
• 内存：仅需16字节（128位状态），MT19937需2.5KB。
• 周期：(2^{128})（可扩展）。

3. 并行能力

• 种子派生：通过SeedSequence.spawn()生成独立子流，避免序列重叠。
• 状态跳跃：jumped(n) 跳过 (2^n) 步状态，快速创建并行流。

三、Python实现
1. 基础用法（NumPy）
NumPy 1.17+ 默认使用PCG64（或PCG64DXSM）：

import numpy as np
# 创建生成器（默认PCG64）
rng = np.random.default_rng(seed=42)
# 生成随机数
uniform = rng.random(size=5) # [0,1) 均匀分布
normal = rng.normal(0, 1, size=5) # 标准正态分布
integers = rng.integers(0, 10, size=5) # [0,10) 整数

2. 并行流管理

# 生成10个独立子流
parent_rng = np.random.default_rng(42)
child_streams = parent_rng.spawn(10)
# 从子流生成随机数
for stream in child_streams:
print(stream.random())

3. 自定义位参数（高级）
若需实现原生PCG64（非NumPy封装）：

class PCG64
:
def __init__(self, seed, a=0x5851F42D4C957F2D, c=0x14057B7EF767814F, p=32):
self.state = seed
self.a = a # 乘数（固定常数）
self.c = c # 增量（奇数）
self.p = p # 位移参数
def next(self):
# 状态更新
self.state = (self.a * self.state + self.c) &
0xFFFFFFFFFFFFFFFF
# 输出变换（取高64位）
u = self.state >>
(64 - self.p)
output = u ^ (u >> self.p)
return output &
0xFFFFFFFF # 返回32位整数
# 使用示例
pcg = PCG64(seed=42)
print(pcg.next()) # 输出随机整数

直方图测试

import numpy as np
import matplotlib.pyplot as plt
seed=0
iters=1000;
v=np.zeros(iters)
for k in range(iters):
pcg= PCG64(seed=seed)
seed=pcg.next()
v[k]=seed/2**32
plt.hist(v)
plt.show()

直方图测试

四、应用场景与最佳实践

1. 蒙特卡洛模拟

   • 金融衍生品定价（百万次路径并行计算）。
   • 物理粒子运动仿真（需高吞吐量随机流）。

2. 机器学习

   • 数据增强（图像/文本扰动）。
   • 模型参数初始化（可复现性要求高）。

3. 游戏开发

   • 地图生成、NPC行为（避免序列重复）。

最佳实践：

• 加密场景：PCG64不适用（非加密安全），改用secrets模块。
• 大规模并行：优先使用PCG64DXSM（NumPy 1.17+默认）。
• 种子管理：用SeedSequence派生子流，而非手动设置种子偏移。

PCG64以简洁的设计平衡了速度、内存与统计质量，成为现代科学计算的首选PRNG。其Python实现依托NumPy的default_rng()接口，兼具易用性与扩展性，尤其适合需高并发和高性能的场景。