天天看點

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