Python_随机序列生成

本文介绍如何利用Python自行生成随机序列，实现了 Whichmann / Hill 生成器。

参考：

　　[1]Random Number Generation and Monte Carlo Methods(P.47)

　　[2]简单产生白噪声的算法

　　[3]各种分布白噪声的产生

基本原理

　　本文粗略将随机数分为两种：均匀分布以及非均匀分布。均匀分布随机数通过非线性变换可得到非均匀分布的随机数。故而均匀分布随机数更基本。引文^[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_0\right)$，是一个三元矢量。该发生器直接产生 $\left(0,1\right)$ 区间上的 $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

巴特西

Python_随机序列生成_白噪声

最新文章

热门文章