天天看点

c语言 算法库 设计 gsl,GSL科学计算库、随机变量的Erlang分布与Weibull分布

更新:本文升级版“GSL科学计算库、随机变量的Erlang分布与Weibull分布”已经迁移至我的新博客http://gnailuy.com/。新文章对已有内容做了修改,并新增关于Weibull分布的内容。欢迎猛击这里:

=====旧文:=====

本文介绍了GSL库和Erlang分布的一般知识,以及使用GSL库生成服从Erlang分布的随机数的方法。对计算机生成的随机数,可以有许多算法对它们进行测试,例如一个常见的测试是Kolmogorov-Smirnov测试,也称作K-S Test。云师姐的这篇文章有使用Matlab进行KS-test的示例以及代码。

GSL及生成随机数的一般方法

GSL(GNU Scientific Library)是GNU组织的数值计算C/C++函数库。它是自由软件,依从GPL协议发布。GSL提供了大量关于数学计算的函数库,当然也包括本文用到的随机数生成函数。更多关于GSL的信息可以到GSL的主页去了解。

计算机中产生服从各种分布的随机数,其基础是产生服从均匀分布的随机数。得到服从均匀分布的随机数以后,可以通过许多不同的算法产生服从其他分布的随机数,例如较常见的使用Polar(Box-Mueller)方法(gsl-1.9/randlist/gauss.c中函数gsl_rand_gaussian)或者使用Ziggurat方法(gsl-1.9/randlist/gausszig.c中函数gsl_rand_gaussian_ziggurat)产生Gaussian分布的随机数等(参考William H.Press等人的著作《C数值算法》)。

服从均匀分布的随机数亦可由许多不同的随机数生成器来产生,不同的随机数生成器生成随机数的速度、随机性等均有差别。GSL库提供了12种随机数生成器(来源)。其中速度最快的是taus、gfsr4和mt19937(default)这三个生成器,而随机性最好的则是ranlux系列算法,也就是GSL的ranlxs系列生成器(来源)。ranlxs系列生成器中,ranlxs0、ranlxs1和ranlxs2产生24位单精度随机数,ranlxd1和ranlxd2产生48位双精度随机数。这五个生成器名字后面的数字代表luxury的程度不同,较高luxury程度的生成器产生的样本数据之间相关程度较低。值得一提的是,计算机中这种使用确定算法产生的所谓随机数,都是伪随机数(参考Knuth的《计算机程序设计艺术》卷二)。然而上述产生伪随机数的生成器由于具其生成的数据具有一定的随机性而得到了广泛的应用。

使用GSL库生成随机数的一般方法如下:

一、创建一个随机数生成器实例

r = gsl_rng_alloc(T);

这里的T是gsl_rng_type类型的指针,它可以是gsl_rng_default(即gsl_rng_mt19937)、gsl_rng_ranlxs0或者gsl_rng_ranlxd1等,用于指定使用不同的随机数生成器。

二、生成种子

gsl_rng_default_seed = ((unsigned long)(time(NULL)));

同一般我们使用C语言中的srand和rand函数一样,通常我们采用当前时间作为种子,以提高随机性。

三、生成服从指定分布的随机数

for (i=0; i

{

u = gsl_rng_uniform(r);

}

for (i=0; i

{

u = gsl_ran_erlang(r, erlang_a, erlang_n);

}

这里可以使用GSL提供的一系列函数生成服从各种不同分布的随机数,例如上述代码生成n个服从Uniform分布的随机数和n个服从Erlang分布的随机数,其中参数r就是上述生成的随机数生成器实例的指针,其他参数将在后面解释。

四、释放随机数生成器

gsl_rng_free(r);

类似于C中的malloc和free函数,这里的随机数生成器实例也需要进行释放,避免内存泄漏。

在Linux下使用GSL库很容易,对应的头文件通常是gsl

u = gsl_ran_erlang(r, erlang_a, erlang_n); // 生成服从 Erlang 分布的随机数

printf("%.5f\n", u);

}

// 打印 0.0, 1.0 处的概率密度值

printf("erlang_pdf(0.0)=%.5f\n", gsl_ran_erlang_pdf(0.0, erlang_a, erlang_n));

printf("erlang_pdf(1.0)=%.5f\n", gsl_ran_erlang_pdf(1.0, erlang_a, erlang_n));

gsl_rng_free(r); // 释放随机数生成器

return 0;

}

使用下面命令编译运行即可:

gcc `pkg-config --cflags --libs gsl` gslErlang.c -o gslErlang

./gslErlang