随机数组合发生器-matlab实现

单个的随机数发生器在周期性等方面存在缺陷，组合发生器是将两个或者多个独立的随机数发生器以某种方式组合起来，使得新组合的随机数发生器具有更长的周期和良好的统计性质。

在组合发生器中比较著名的是M.D.Maclaren和G.Marsaglia于1965年提出的组合同余法，生成规则如下：

1. 采用第一个线性同余发生器M1生成k个随机数，把k个数按序依次存放在某一个向量T中，T=(t1,t2,…,tk)，令n=1；
2. 用第二个线性同余发生器M2生成一个随机整数j，满足1 <= j <= k；
3. 令xn = tj，用M1生成一个新的随机数y替代tj，即令tj = y，置n = n+1；
4. 重复2-3步，得到随机数序列{xn}，即为组合发生器生成的序列。若M1的模为m，令 un = xn/m，则{un}即为由该组合发生器生成的均匀随机数序列。

当然，根据此生成规则以及常用的随机数发生器可以生成多个组合发生器，如平方取中发生器与其自身的组合、与乘积取中、线性同余、常数乘子、斐波那契等随机数发生器的组合(单个随机数发生器的matlab代码可参考文章：五种常用的随机数发生器-matlab实现)。下面用matlab实现其中的一种组合。

注：下面直接在程序里进行了参数调用，也可将整个过程封装成函数然后调用。

% 组合发生器生成随机数：
% 为简便起见，循环之前首先生成了满足规则2和规则3的两组向量，
% 在循环体中从两组向量中随机抽取满足条件的整数j和y，只是其中一种解决方式，仅做参考
clc
close all

% 第一步，第一个随机数发生器生成100个随机数，并按序放至向量T中
% 生成随机向量Y，从Y中随机取出一个y来替代tj
[T,u1] = linear_mod_random(73, 7, 213, 5, 100);    % 参数可更改
Y = linear_mod_random(53, 2, 157, 6, 100);         % 参数可更改

# 将三个随机数序列放在一起可视化
subplot(3,2,1);
plot(T(1:end), '*-');
title('第一个随机数发生器');
hold off;

subplot(3,2,2);
hist(u1(2:end));
title('随机数直方图分布');
hold off;

% 第二步，第二个随机数发生器生成随机的整数j，满足1<=j<=100
% 生成随机向量J，在循环中从向量J中随机抽取一个满足条件的j
[J,u2] = linear_mod_random(47, 2, 157, 6, 100);      % 参数可更改

subplot(3,2,3);
plot(J(1:end), '*-');
title('第二个随机数发生器');
hold off;

subplot(3,2,4);
hist(u2(2:end));
title('随机数直方图分布');
hold off;

J1 = [];
u21 = [];
for i = 1 : 100
    % 第二步，第二个随机数发生器生成随机的整数j，满足1<=j<=100
    j = rem(J(int16(rand * (length(J) - 1) + 1)),100)+1;
    % 第三步，第一个随机数发生器生成一个新的随机数y以替代tj
    y = Y(int16(rand * (length(Y) - 1) + 1));
    % 第三步，将t[j]赋值给x[i],并用y代替t[j]
    x(i) = T(j);
    T(j) = y;
end

u = x / 213;   # 此处求u用的是第一个随机数发生器的m值

subplot(3,2,5);
plot(x(1:end), '*-');
title('组合随机数发生器');
hold off;

subplot(3,2,6);
hist(u(2:end));
title('随机数直方图分布');
hold off;

% 定义linear_mod_random函数
function [x,u] = linear_mod_random(a, c, m, x0, iterMax)
    x = [];
    u = [];
    x(1) = x0;
    u(1) = x(1) / m;
    for i = 2 : iterMax
        x(i) = a * x(i-1) + c;
        x(i) = mod(x(i), m);
        u(i) = x(i) / m;
    end
end

需要注意的是，组合发生器理论上能够提高随机数的性能，但也会受到组合规则、组合所用发生器自身等因素的影响，性能是否提高需要利用随机数序列检验方法进行检验。