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/*.h,编译时只需要在gcc命令中添加相应参数即可,可以参考下面两个命令:
gcc `pkg-config --cflags --libs gsl` example.c -o example
gcc -lgsl -lgslcblas -lm example.c -o example
在Windows中使用GSL库要麻烦一些,如果你使用MingW的话,可以到这里下载GSL的库文件。GSL为MingW和Cygwin等环境提供了.a格式的库,使用起来相对还要容易一些,但是如果你使用Virtual C++或者Virtual Studio的话,就要麻烦一些。这里展示一个使用Virtual C++ 6.0的可行的解决方案,如下:
到http://www6.in.tum.de/~kiss/WinGsl.htm下载WinGsl-Lib-1.4.02.zip,解压到某个目录下面,例如C:\WinGsl;
打开Virtual C++,点击菜单Tools/Options...,在弹出的对话框里选择Directories选项卡;
在Show directories for标签下选择Include files,然后在下面的Directories里添加刚刚解压的目录,这里是C:\WinGsl;
在Show directories for标签下选择Library files,然后在下面的Directories里添加目录C:\WinGsl\Lib;
复制C:\WinGsl\Bin下的两个文件WinGsl.dll和WinGslD.dll到Virtual C++安装目录里的相应目录,例如C:\Program Files\Microsoft Visual Studio\VC98\Bin;
建立工程,打开菜单Project/Settings...,在弹出的对话框选择Link选项卡,工程需要用到哪些库里的函数,就把要用到的库文件名添加到Object/library modules下面的文本框里,通常添加文件WinGslLib_s.lib即可。
Erlang分布
Erlang分布是一种连续型概率分布,与指数分布一样多用来表示独立随机事件发生的间隔,Erlang分布不具有马尔科夫性。其概率密度函数如下:
Erlang分布有两个参数,一个是阶数(stage)n,一个是均值[tex]\mu[/tex](或用[tex]\lambda=\frac{1}{\mu}[/tex])。阶数n=k的Erlang分布称作K-阶Erlang分布。
遵循Erlang分布的随机变量可以分解成多个同参数指数分布的随机变量之和,当阶数n=1时,Erlang分布就退化为指数分布。
此事可以直观理解,柏松过程的两个相邻事件到达时间间隔服从指数分布,而从0时刻直到第n个事件发生的到达时间则服从Erlang分布。
设一柏松过程的到达时间间隔序列为[tex]\{X_n, n\geq 1\}[/tex],则若[tex]X_n[/tex]服从参数为[tex]\lambda[/tex]的指数分布,到达时间[tex]S_n = \sum_{i=1}^nX_i[/tex]就服从参数为[tex](n, \lambda)[/tex]的Erlang分布。
Erlang分布是Gamma分布的特殊情况,Gamma分布的概率密度函数为:
其中Gamma函数[tex]\Gamma(\alpha)[/tex]可以看作是阶乘的扩展,它有下面的特征(来源):
当[tex]\alpha[/tex]为整数时,Gamma分布就被称为Erlang分布。后面GSL库中Erlang分布的计算其实就是利用Gamma分布的计算函数进行的。
GSL库产生Erlang分布的随机数
GSL的源代码文件erlang.c中提供了两个关于Erlang分布的函数:
double gsl_ran_erlang (const gsl_rng * r, const double a, const double n);
double gsl_ran_erlang_pdf (const double x, const double a, const double n);
上述两个函数的参数r为随机数生成器指针,a代表Erlang分布概率密度函数中的[tex]\frac{1}{\lambda}=\mu[/tex],n就是阶数。函数gsl_ran_erlang产生服从Erlang分布的随机数,而gsl_ran_erlang_pdf计算参数为(n,a)的Erlang分布在x处的概率密度。
产生服从Erlang分布的随机数
函数gsl_ran_erlang产生随机数,服从参数为(n, a)的Erlang分布,这个函数产生随机数的基础是r指向的随机数生成器生成的Uniform分布的随机数。在实现上,这个函数将传给它的三个参数原样传给下面的Gamma分布计算函数计算:
double gsl_ran_gamma (const gsl_rng * r, const double a, const double n);
gsl_ran_gamma首先确保[tex]a\geq1[/tex],如果a<1,则先生成一个(0,1)区间内服从Uniform分布的随机数u,然后返回:
gsl_ran_gamma (r, 1.0 + a, n) * pow (u, 1.0 / a);
当[tex]a\geq1[/tex]时,gsl_ran_gamma首先循环调用函数gsl_ran_gaussian_ziggurat,使用Marsaglia-Tsang ziggurat方法产生方差[tex]\sigma=1.0[/tex],服从Gaussian分布的随机数x,然后计算[tex]v=1+\frac{x}{3\sqrt{a-\frac{1}{3}}}[/tex],直到[tex]v\geq0[/tex]为止。
得到满足条件的v以后,计算[tex]v=v^3[/tex],然后再产生一个(0,1)区间内服从Uniform分布的随机数u,如果v满足下面两个条件中的任意一个,就接受这个v,否则重复上述生成v的整个过程,直到生成合格的v为止:
最后,函数返回[tex]nv(a-\frac{1}{3})[/tex]。
计算概率密度
函数gsl_ran_erlang_pdf计算x处Erlang分布的概率密度,其计算公式为:
其中lngamma是调用函数gsl_sf_lngamma计算的[tex]ln \Gamma(n)[/tex]。其实上式就是由上面Gamma函数的概率密度函数得到的,进行一下简单推导不难验证。
附代码样例:使用GSL库生成Erlang分布的随机数
#include
#include
#include
#include
#include
#include
#define MAXRNDNUM 100
int main(int argc, char * argv[])
{
const gsl_rng_type *T; // 随机数生成器类型指针
gsl_rng *r; // 随机数生成器指针
int i;
double u; // 随机数变量
const double erlang_a=0.6; // lambda=0.6
const double erlang_n=2.0; // n=2.0
// gsl_rng_default (gsl_rng_mt19937)
// gsl_rng_ranlxs0, gsl_rng_ranlxs1, gsl_rng_ranlxs2
// gsl_rng_ranlxd1, gsl_rng_ranlxd2
T = gsl_rng_ranlxs0;
gsl_rng_default_seed = ((unsigned long)(time(NULL))); // 取当前时间作为种子
r = gsl_rng_alloc(T); // 创建随机数生成器实例
for (i=0; i
{
/** Functions:
* double gsl_ran_erlang (const gsl_rng * r, const double a, const double n)
* double gsl_ran_erlang_pdf (const double x, const double a, const double n)
*/
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
c语言 算法库 设计 gsl,GSL科学计算库、随机变量的Erlang分布与Weibull分布相关推荐
- python的科学计算库有哪些_Python科学计算库-Numpy
NumPy 是 Python 语言的一个扩充程序库.支持高级大量的维度数组与矩阵运算,此外也针对数组运算提供大量的数学函数库,也是学习 python 必学的一个库. 1. 读取文件 numpy.gen ...
- 【Python 标准库学习】数据科学计算库 — math
欢迎加入 Python 官方文档翻译团队:https://www.transifex.com/python-doc/ math 模块官方文档:https://docs.python.org/3/lib ...
- GSL科学计算库——计算高斯-勒让德积分
相关文章: Windows系统Qt5配置GSL科学计算库 Windows系统下Cygwin+Dev C ++ 配置GSL科学计算库 假设计算下列积分: ∫0πexcos(x)dx\int_0^\pi ...
- 高精度计算器_尝试设计个.NET Core高精度科学计算库
猪年快到了,还记得在我在鸡年立了个flag,计划在狗年推出<纸书科学计算器>的2.0版本(HowardZ:近来的一些修修补补).然而现实很骨感,我因为时间太紧而咕咕了--就此陷入了人类三大 ...
- python科学计算库-数值计算库与科学计算库
BLAS 接口 BLAS , LAPACK , ATLAS 这些数值计算库的名字很类似,他们之间有什么关系呢?BLAS是一组线性代数运算接口,目前是事实上的标准,很多数值计算/科学计算都实现了这套接口 ...
- python科学计算主要学什么_以下哪些是python常用的科学计算库?_学小易找答案
[单选题]17-51. 在 Windows 中,若要终止未响应的应用程序,可使用( ) [单选题]19-55.在 Windows控制面板中,下列无法实现的操作是 [单选题]witness的元素属性(比 ...
- 【Python基础】科学计算库Scipy简易入门
0.导语 Scipy是一个用于数学.科学.工程领域的常用软件包,可以处理插值.积分.优化.图像处理.常微分方程数值解的求解.信号处理等问题.它用于有效计算Numpy矩阵,使Numpy和Scipy协同工 ...
- 初识 Python 科学计算库之 NumPy(创建多维数组对象)
文章目录 参考 描述 NumPy 特点 获取 导入 多维数组对象 np.array() np.asarray() 范围 随机 概览 np.random.randn() np.random.normal ...
- numpy不用科学记数发 python_Python科学计算库Numpy常用的函数使用
林小森博客: Python科学计算库Numpy常用的函数使用 - 林小森www.linxiaosen.com Numpy具有强大的计算功能,本文介绍Numpy常用的函数,可以有效的提高工作效率. 首 ...
最新文章
- OCH\OMS\OTS\MSP\RS\SPI解释
- 51nod 1022 石子合并v2
- 【Java中级】(三)IO
- 前端学习(2064):vue的生命周期函数有什么
- 应用虑镜特效时遇到浏览器权限问题
- python处理二进制文件_使用Python进行二进制文件读写的简单方法(推荐)
- DesignPatterns-装饰器模式
- 物联网技术的基站能耗监控解决方案
- 微信小程序报错 40125 已解决
- HTML5游戏实战 1 50行代码实现正面跑酷游戏
- 实现发送xml格式的请求
- 六轴机器人运动学正解
- java通过framer生成word_DSO Framer Control Object 操作word文件
- 2-6轴点胶控制系统方案设计
- iwlwifi(AC9260)移植总结
- 方舟服务器金币系统,方舟指令快速获取金币方法大揭晓 你也可以变土豪
- 软件安全之代码注入技术 向目标 PE 文件注入 DLL notepad lpk.dll 远程线程函数 提权函数 OpenProcess VirtualAllocEx
- 【转】GPS定位基本原理浅析
- 如何提升微信小程序排名的技巧?
- 关于discuz论坛邮箱配置
热门文章
- 计算机组装与维护 游戏设计,《计算机组装与维护课程设计报告.doc
- WIN10快捷键探索
- 江南爱软装十大品牌 软装的风格你都知道多少?
- android适配华为m5,2019-05-29 Android悬浮窗适配全机型,包含8.0,小米魅族华为悬浮窗权限适配demo看这一篇就够了...
- 利用PCL库从点云数据生成深度图像及关键点提取
- python调用origin画图_Python科学绘图
- 【笔记】机器学习所涉及到的“积分学”知识
- 合肥先进光源束测后台的初步设计
- postgresql_internals-14 学习笔记(五)Buffer Cache
- 视频抽帧:多视频、可视化、手动旋转