【概率论】均匀分布的伪随机数

伪随机数的基本概念

伪随机数（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++ 不适用于密码学（非加密安全），但适合高性能模拟。

五、应用场景建议

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()接口，兼具易用性与扩展性，尤其适合需高并发和高性能的场景。