本章主要介绍NS-3当中的(伪)随机数生成器。NS-3提供了简单方便的随机数类型,可以生成各种分布的随机数。随机数生成器也属于NS-3对象框架的类,适用一切对象框架中的特性。此外,随机数可以作为属性值,并且方便地从字符串创建。为了生成互不相交的随机数序列,NS-3使用了随机数流的概念,保证不同的随机数生成器的实例当中,生成的随机数尽量不存在相互交叉的情况。
这里
)。它能够生成的随机数序列的周期为$3.7\times 10^{57}$。这个周期是什么概念?夸张一点,假设我们可以一纳秒生成一个随机数,那么要达到这个周期需要约$1.189\times 10^{41}$年。可以对比的是根据大爆炸理论,现在宇宙的年龄为$1.382\times 10^{10}$……到宇宙毁灭都用不完这么多随机数啊。因此,我们可以将这么长的序列分成多段,形成不同的子流。L’Ecuyer等人将随机数流划分为$1.8\times 10^{19}$个相互独立的流,而每个流又划分为$2.3\times 10^{15}$,每个子流的长度为$7.6\times 10^{22}$。需要$2.443\times 10^{6}$年,即约20万年可以用完这个随机数子流。
因此,对于上面提到的三个衡量准则,该算法:
可以选择不同的流来获得随机性,同时,可以固定一个流来获得确定性。
对于一个特定的流,不同的随机数生成器实例可以选择不同的子流来保证相互的独立性。
随机数生成器的周期长到可以为所欲为。
如果你足够细心的话可能还记得我们在配置路径一章提到过NS-3提供了5个全局属性。其中和随机数生成的属性一个有两个:RngSeed和RngRun。其中RngSeed即表示随机数种子。随机数的种子设置不同的值的话,可以保证随机数生成不同的序列,这可以用以保证随机性。但如果需要重复相同的程序结果,可以保持随机数种子不变,以保证确定性。除此之外,NS-3还提供了另外一个属性RngRun,通过修改这个参数,也可以改变随机数序列的随机性。那么这两个属性究竟有什么联系和区别呢?
使用RngSeed和RunRun的一种典型的情况是对某个特定的仿真过程进行多次独立试验。以便根据大量的独立运行结果来计算统计信息。用户可以改变全局种子并再次运行该模拟过程,也可以推进随机数的流的状态,这被称为增加运行次数(run number)。
那么设置新的种子和推进子流的状态,哪一个更好?当改变随机数种子的时候,我们无法保证流产生的两个随机种子不是互相重叠的,即两个不同的种子生成的随机数之间可能有一部分是相互重叠的。保证两个流不互相重叠的唯一办法是使用由生成算法实现提供的子流功能。使用子流功能来实现同一个模拟过程的多个彼此独立的运行。从统计意义上来说,设置独立的重复实验更严格的方法应该是使用固定的种子并改变运行次数(RngRun)。使用子流时,独立的重复实验的最大值为$2.3×10^{15}$(因为每个流只有这么多个子流)。
总结:不同的随机数生成器实例对应不同的流,因此,各个生成器实例之间的随机数序列是相互独立的。不同的运行次数之间,某个生成器实例又使用不同的子流,因此,可以保证每个生成器实例在每次运行之间生成的随机数序列是独立的。因此更加推荐在重复的多次仿真当中使用运行次数的推进。
除了设置全局属性之外,NS-3提供了另外一种方法来设置随机数的种子和运行次数:
1 2
RngSeedManager::SetSeed (3 ); RngSeedManager::SetRun (7 );
这和设置全局属性的效果是一样的。当然,设置随机数种子和运行次数需要在程序刚开始创建任何随机变量之前进行。
当然这种方式每次修改种子和运行次数之后都会引起程序重新编译。由于RngSeed和RngRun都是全局属性,因此,我们也可以使用命令行一章介绍的方法来通过命令行参数修改他们的值:
1
$ ./waf --command-template ="%s --RngSeed=3 --RngRun=7" --run program-name
这种方式无论多少次修改随机数种子或者运行次数都不会引起程序重新编译。
马欢的博客
):
try-uniform-random.cc
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
#include "ns3/core-module.h" #include <iostream> using namespace ns3;NS_LOG_COMPONENT_DEFINE ("TryUniformRandom" ); int main (int argc, char *argv[]) { LogComponentEnable("TryUniformRandom" , LOG_LEVEL_INFO); const double MIN = 0.0 ; const double MAX = 10.0 ; const uint32_t MAX_NUMBERS = 100000 ; Ptr<UniformRandomVariable> x = CreateObject<UniformRandomVariable> (); x->SetAttribute ("Min" , DoubleValue (MIN)); x->SetAttribute ("Max" , DoubleValue (MAX)); std ::ofstream outfile ("uniform.txt" ) ; for (uint32_t i = 0 ; i < MAX_NUMBERS; i++) { double value = x->GetValue (); outfile << value << "\n" ; } outfile.flush(); outfile.close(); Simulator::Run (); Simulator::Destroy (); }
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
set term pngcairo lw 2 set output "uniform.png" n=50 #number of intervals max=0.0 #max value min=10.0 #min value width=(max-min)/n #interval width #function used to map a value to the intervals hist(x,width)=width*floor(x/width)+width/2.0 set boxwidth width*0.9 set style fill solid 0.5 # fill style #count and plot plot "uniform.txt" u (hist($1,width)):(1.0) smooth freq w boxes lc rgb"green" notitle
当然生成的数字越多,越符合大数定理,越能展示其分布。最终直方图如下:
博客
。
对于离散型的随机变量,可以使用通用方法生成随机变量,我们只需要知道其分布函数即可。例如,我们需要生成几何分布的随机变量,其分布函数为:$Pr(k) = (1-p)^(k-1)p$,即随机变量的值是k的概率。那么每个取值都对应一个概率。因此,我们可以计算一张表,例如p=0.5的时候,有:
因此,我们只需要生成一个(0,1]分布的随机变量值,然后写一个循环,查找这个值究竟在哪个范围内即可确定对应的值是多少。然而这种方法太慢了,当变量的取值范围非常大的时候,需要循环很多次才能找到正确的区间。
不同的离散型随机变量生成方法大多都有不同的优化方法。例如,我们要生成几何分布的随机变量。首先,我们需要先明白几何分布的随机变量的意义才能更好地进行优化。这里,几何分布表示反复进行贝努利实验,第k次才实验成功(前k-1次都失败)的概率。
什么是贝努利实验呢?其实可以简单的理解,这个实验只有两种结果,并且其中一种的概率用p来表示。例如,抛硬币就是一个典型的贝努利实验,只有正面和反面两种结果,我们可以将结果为正面的概论定为p,那么结果为背面的概率自然就是1-p。
因此几何分布中,第k次才成功的概率为:$P(k)=(1-p)^{k-1}p$。其中$(1-p)^{k-1}$表示前k-1次均失败的概率,而p表示第k次成功的概率。
$\sum_{i=1}^{j-1}P(i)$表示第1次成功的概率一直累加到第j-1次成功的概率,其中$j\ge 1$。因此有$\sum_{i=1}^{j-1}P(i) = 1 - \sum_{i=j}^{+\infty}P(i)$。而$\sum_{i=j}^{+\infty}P(i)$表示前$j-1$次都未成功的概率,根据贝努利实验的原理可知$\sum_{i=j}^{+\infty}P(i)=(1-p)^{j-1}$。因此最终有$\sum_{i=1}^{j-1}P(i)=1-(1-p)^{j-1}$。这就将概率的累加用公式表示出来了,因此,我们编程的时候就无需循环去累加概率了。和通用方法类似的,我们只需要去生成一个(0,1]分布的随机变量u满足如下公式:
$$1-(1-p)^{j-1}\le u < 1-(1-p)^{j}$$
那么对应的j就是我们要求的值。简化后得到:
$$(1-p)^{j}< 1- u \le (1-p)^{j-1}$$,即
$$X=\min\{j:(1-p)^{j}<(1-u)\}=\min\{j:ln(1-p)^{j}<ln(1-u)\}=\min\{j:jln(1-p)<ln(1-u)\}$$
这里需要注意的是,由于(1-p)一定小于1,则ln(1-p)一个是个负数。因此有:
$$X=\min\{j:j>\frac{ln(1-u)}{ln(1-p)}\}=\lfloor\frac{ln(1-u)}{ln(1-p)}\rfloor+1$$
由于u是均匀分布的随机变量,那么生成u和生成1-u的概率是一样的,因此,有
$$X=\lfloor\frac{ln(u)}{ln(1-p)}\rfloor+1$$。
这里需要注意的是p不能等于1,程序当中需要特殊处理一下。不过从意义上来理解,如果p等于1了,那说明成功的概率为1,肯定实验第一次就成功了呀。所以……你懂的
这样生成的几何分布的均值为$E(X)=1/p$,方差为$D(X) = (1-p)/p^2$。
有了以上的公式,我们就可以生成几何分布的随机变量了。只要写一个类继承RandomVairableStream,然后按公式生成随机数即可。程序如下:
try-geometric-random.cc
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
#include "ns3/core-module.h" #include <iostream> using namespace ns3;NS_LOG_COMPONENT_DEFINE ("TryGeometricRandom" ); class GeometricRandomVariable : public RandomVariableStream {public : static TypeId GetTypeId (void ) ; GeometricRandomVariable (); double GetProbability (void ) const ; virtual double GetValue (void ) ; virtual uint32_t GetInteger (void ) ; private : double m_probability; }; NS_OBJECT_ENSURE_REGISTERED(GeometricRandomVariable); TypeId GeometricRandomVariable::GetTypeId (void ) { static TypeId tid = TypeId ("ns3::GeometricRandomVariable" ) .SetParent<RandomVariableStream>() .AddConstructor<GeometricRandomVariable> () .AddAttribute("Probability" , "The lower bound on the values returned by this RNG stream." , DoubleValue(0.5 ), MakeDoubleAccessor(&GeometricRandomVariable::m_probability), MakeDoubleChecker<double >(0 , 1 )) ; return tid; } GeometricRandomVariable::GeometricRandomVariable () { NS_LOG_FUNCTION (this ); } double GeometricRandomVariable::GetProbability () const { return m_probability; } double GeometricRandomVariable::GetValue (void ){ return GetInteger(); } uint32_t GeometricRandomVariable::GetInteger (void ){ if (m_probability == 1 ) { return 1 ; } double u = Peek()->RandU01(); return (int )(log (u) / log (1 - m_probability)) + 1 ; } int main (int argc, char *argv[]) { LogComponentEnable("TryGeometricRandom" , LOG_LEVEL_INFO); const double PROBABILITY= 0.5 ; const uint32_t MAX_NUMBERS = 10 ; Ptr<GeometricRandomVariable> x = CreateObject<GeometricRandomVariable> (); x->SetAttribute ("Probability" , DoubleValue (PROBABILITY)); for (uint32_t i = 0 ; i < MAX_NUMBERS; i++) { double value = x->GetValue (); NS_LOG_INFO(value); } Simulator::Run (); Simulator::Destroy (); }
运行程序,得到如下结果:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Waf: Entering directory `/home/rainsia/Applications/ns-allinone-3.29/ns-3.29/build' [1980/2049] Compiling scratch/try-geometric-random.cc [2009/2049] Linking build/scratch/try-geometric-random Waf: Leaving directory `/home/rainsia/Applications/ns-allinone-3.29/ns-3.29/build' Build commands will be stored in build/compile_commands.json 'build' finished successfully (7.878s) 1 1 2 2 2 1 4 1 1 1
同样的,我们也生成10万个符合几何分布的随机数,然后画出其直方图:
从图中可以看到,大部分时候,实验都是一两次就成功了。因为默认的成功概率p为0.5。假如我们认为成功是指抛硬币抛到正面,那么,这就解释了为什么我们抛硬币可以在短短几次之内就能抛到正面。