本文介绍如何利用Python自行生成随机序列,实现了 Whichmann / Hill 生成器。
参考:
基本原理
本文粗略将随机数分为两种:均匀分布以及非均匀分布。均匀分布随机数通过非线性变换可得到非均匀分布的随机数。故而均匀分布随机数更基本。引文[3]中提到了三种生成算法:线性同余法、联合法、反馈位移寄存法。限于时间,简述线性同余法:
线性同余是基于求余计算的方法,对原始序列进行重排。该方法由一个“种子”开始,通过递归得到整个随机序列[3]:
begin{equation}
x_i = ax_{i-1}+c quad left(mod m quad right)
end{equation}
式中 $x_i$、$a$、$c$ 均为正整数,$x_0$ 称为种子,$a$ 称为乘子,$c=0$ 时称为乘同余法,反之称为混合同余法。
查阅引文[3]的书,全英文,篇幅长,但还是翻了翻。翻到47页的时候发现了很清晰的算法表述。该算法在1982年由 Wichmann 和 Hill 发表,通过三个发生器组合成一个同余发生器。因为其易于编程,且随机性优越(书中简述了相关研究),本文选择该算法进行实现[1]:
begin{equation}
begin{cases}
x_i equiv 171x_{i-1}mod 30269 \
y_i equiv 172x_{i-1}mod 30307\
z_i equiv 170x_{i-1}mod 30323\
u_i = left( dfrac{x_i}{30269} + dfrac{y_i}{30307} + dfrac{z_i}{30323} right) mod 1
end{cases}
end{equation}
其中,发生器的种子为 $left(x_0,y_0,z_0right)$,是一个三元矢量。该发生器直接产生 $left(0,1right)$ 区间上的 $u_i$ 。可以由最终表达式看出:分母的值很大。
实现
为了减少代码量,我将 $x,y,z$ 三组数都放到了一个变量下面,名字为 x 。a 与 b 设为 np.array 格式,避免序列之间不支持对应元素的乘除求余问题。
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 28
@author: adgk07
A script to generate noise signal
"""
import numpy as np
from pylab import *
S = [1,3,7] # user defin SEED
L = 1024 # user defin LENGTH
a = np.array([171., 172, 170])
b = np.array([30269., 30307, 30323])
x = np.zeros((L,3))
x[0] = S
y = np.zeros(L)
y[0] = np.sum(S/b)%1
for num in np.arange(L-1):
x[num+1] = (a*x[num])%b
y[num+1] = np.sum(x[num+1]/b)%1
plot(y)
show()
上述例程可改写为函数,不在话下。
END