現在,在計算機,用來產生隨機數的算法是“線性同余”法。所謂線性同余,其實就是下面兩個式子。假設I就是壹個隨機數的序列,Ij+1與Ij的關系如下:
Ij+1 =Ij * a+c (mod m)
或是Ij+1 =Ij *a (mod m),
其中,不妨取a=16807,m=2147483647,以為壹常數。寫個簡單的程序就是:
long r;
void scand( long v)//初始化隨機種子數
{
r = v;
}
long rand()//產生隨機數
{
r = (r*a + c)%m;//a,c,m為常數
return r;
}
再看壹下稍復雜壹點的:(Random () 的 Borland 的實現)
long long RandSeed = #### ;
unsigned long Random(long max)
{
long long x ;
double i ;
unsigned long final ;
x = 0xffffffff;
x += 1 ;
RandSeed *= ((long long)134775813);
RandSeed += 1 ;
RandSeed = RandSeed % x ;
i = ((double)RandSeed) / (double)0xffffffff ;
final = (long) (max * i) ;
return (unsigned long)final;
}
二.計算機產生的隨機數不是真的隨機數
[引:]我們建立了真正調用偽隨機數生成器的 random()。但什麽是偽隨機數生成器?假定需要生成介於 1 和 10 之間的隨機數,每壹個數出現的幾率都是壹樣的。理想情況下,應生成 0 到 1 之間的壹個值,不考慮以前值,這個範圍中的每壹個值出現的幾率都是壹樣的,然後再將該值乘以 10。請註意,在 0 和 1 之間有無窮多個值,而計算機不能提供這樣的精度。
為了編寫代碼來實現類似於前面提到的算法,常見情況下,偽隨機數生成器生成 0 到 N 之間的壹個整數,返回的整數再除以 N。得出的數字總是處於 0 和 1 之間。對生成器隨後的調用采用第壹次運行產生的整數,並將它傳給壹個函數,以生成 0 到 N 之間的壹個新整數,然後再將新整數除以 N 返回。這意味著,由任何偽隨機數生成器返回的數目會受到 0 到 N 之間整數數目的限制。
在大多數的常見隨機數發生器中,N 是 232? (大約等於 40 億),對於 32 位數字來說,這是最大的值。換句話說,我們經常碰到的這類生成器能夠至多生成 40 億個可能值。而這 40 億個數根本不算大,只是指尖這麽大。
偽隨機數生成器將作為“種子”的數當作初始整數傳給函數。這粒種子會使這個球(生成偽隨機數)壹直滾下去。偽隨機數生成器的結果僅僅是不可預測。由偽隨機數生成器返回的每壹個值完全由它返回的前壹個值所決定(最終,該種子決定了壹切)。如果知道用於計算任何壹個值的那個整數,那麽就可以算出從這個生成器返回的下壹個值。
結果,偽隨機數生成器是壹個生成完全可預料的數列(稱為流)的確定性程序。壹個編寫得很好的的 PRNG 可以創建壹個序列,而這個序列的屬性與許多真正隨機數的序列的屬性是壹樣的。例如:
PRNG 可以以相同幾率在壹個範圍內生成任何數字。
PRNG 可以生成帶任何統計分布的流。
由 PRNG 生成的數字流不具備可辨別的模式。
PRNG 所不能做的是不可預測。如果知道種子和算法,就可以很容易地推算出這個序列。
計算機產生的隨機數壹般都只是壹個周期很長的數列,不是真的隨機數。也就是說,隨機數壹般是偽隨機數,每個隨機數都是由隨機種子開始的壹個已定的數列(周期很長)。壹般地,為了隨機數更真壹點,隨機種子在系統中通常是參照系統時鐘生成的。
看看下面這個例子,從中,讀者應能有所體會。
main()
{
int i;
scand(time(NULL)); /*可向計算機讀取其時鐘值,並把值自動設為隨機數種子*/
for(i=0;i<10;i++){
printf("%10d",1+(rand()%6));/*這裏1是移動值,他等於所需的連續整數*/
} /*數值範圍的第壹個數;b是比例因子;他*/
return(0); /*等於所需的連續整數值的範圍寬度;*/
}
從數學上講,為了得到壹個周期性更長的隨機序列,我們可以使用如下方法:(這是我在壹個書上看到的,詳細的情況大家可以查壹查,我自己也記不清了,呵呵)
float rand(long* idum)
{
int j; long k; static long iy=0; static long iv[NTAB];
float temp;
if(*idum<=0||!iy)
{
if(-(*idum)<1) *idum=1;
else *idum=-(*idum);
for(j=NTAB+7;j>=0;j--){
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-k*IR;
if(*idum<0)
*idum+=IM;
if(j<NTAB)
iv[j]=*idum;
} iy=iv[0]; }
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-k*IR;
if(*idum<0) *idum+=IM;
j=iy/NDIV;
iy=iv[j];
iv[j]=*idum;
if((temp=float(AM*iy))>RNMX) return float(RNMX);
else return temp;
}
[註:]有文講 ,根據Randomize的工作原理,利用函數Timer得到從午夜開始到現在經過的秒數,然後再根據要得到的隨機數值大小對該數值進行"“衰減”處理,這樣得到的數值則可稱得上是真正意義的隨機數值,我認為,這也是人為的方法,仍有它的確定性和周期性,仍稱不上是真正的隨機數,單純改變偽隨機數的生成邏輯計算方法並不能達到目的,最有效的辦法只能是改變rand_seed,就是種子。而且,改變後的rand_seed不應該是人為的。(註:目前,在 Intel 的PIII處理器中內置了壹個與CPU溫度相關的隨機數生成器,算是壹個比較有效的種子生成器。)更好的辦法是根據“隨機事件”生成隨機數,如鍵盤和鼠標輸入值、中斷、磁盤讀取等等。然而,許多服務器沒有鍵盤和鼠標,許多“黑盒”產品也不帶有硬盤,因此很難找到好的隨機數源,當然,通訊密鑰也就壹樣不安全。而如網絡狀態等也不能很好地保證隨機數的“隨機性”。電器噪聲和聲音頻率也許是很好的隨機數源,但大多數人恐怕並不願意在計算機中增加這種設備,而且也可能出現設備失靈和外部操縱(影響)等問題。對於要處理大量連接的網關服務器,是必須要考慮的問題。如果可以通過,精確檢測機器cpu的通電電流強度,來作為隨機數種子,或是其他壹些沒有人為因素的幹擾的,且瞬間變化快的方法獲得種子,必將能產生符和要求的隨機數。
三.幾種隨機數的獲取辦法
1.產生壹個範圍內的隨機數
壹般地,我們可用j=1+(int)(n*rand()/(RAND_MAX+1.0))來生成壹個0到n之間的隨機數。
若用int x = rand() % 100;來生成 0 到 100 之間的隨機數這種方法是不可取的,比較好的做法是: j=(int)(100.0*rand()/(RAND_MAX+1.0)) ,當然,我們也可是使用random(100)。下面的例子都是用了random(n).
2、篩選型隨機數 如希望取0-99的隨機數,但不能是6。
解決方法:
x = random(100);
while (x==6) {
x = random(100);
}
又如希望取0-99的隨機數,但不要5的倍數 解決方法:
x = random(100);
while ((x % 5)==0) {
x = random(100);
}
3、從連續的壹段範圍內取隨機數。
如從40--50的範圍內取隨機數。 解決方法: x=random(11)+40
4、從壹組亂數中取隨機數。 如:從 67, 87, 34, 78, 12, 5, 9, 108, 999, 378十個數中隨機取數。 解決方法:可以用數組將些十個數存貯,然後把0--9中取出的隨機數作為序號,實現隨機取數。
a = new Array(67, 87, 34, 78, 12, 5, 9, 108, 999, 378);
j = random(10);
x = a[j];
四.產生具有壹定分布的隨機數
關於這點,有壹篇文章寫得很好,請看:
非均勻分布隨機數的產生及其在計算機模擬研究中的應用
作者:胡性本 劉向明 方積乾
單位:胡性本(湖北職工醫學院 荊州434000); 劉向明(湖北職工醫學院 荊州434000 現為華中理工大學生物工程系博士研究生); 方積乾(湖北職工醫學院 荊州434000 中山醫科大學)
關鍵詞:非均勻分布隨機數;計算機模擬;程序設計
摘 要:非均勻分布隨機數在進行計算機模擬研究中有重要作用,但計算機高級語言通常都只提供產生均勻分布隨機數的函數,給研究工作帶來困難。該文提出的方法,較好地解決了這壹問題,有很強的適用性。
中圖分類號:O 242.1
文章編號:1004-4337(2000)01-0059-02▲
1 問題的提出
常用的計算機高級程序設計語言,大多提供了產生在〔0,1〕區間內連續均勻分布的獨立隨機數r的函數。若將產生的隨機數作簡單的變換X=a+(b-a)r,就能得到在區間〔a,b〕上均勻分布的隨機數X。如果與取整或舍入函數結合運用,還可得到離散均勻分布的隨機數,給計算機模擬提供了方便。但在很多情況下進行計算機模擬需要用到按某些特定規律分布的非均勻隨機數。例如,在離子通道門控模型的模擬研究中,就經常要用到服從指數分布、對數分布、正態分布、普畦松分布等非均勻分布的隨機數。這時可以在源程序中編寫壹個函數或過程來產生。我們在近年來所進行的計算機模擬研究工作〔1,2〕中,大量使用了非均們分布的隨機數,獲得成功。由於離散隨機數可以由連續隨機數通過取整或舍入等途徑轉換而得,本文只側重介紹非均勻分布連續隨機數的產生原理和應用實例,並給出了用類PASCAL語言〔3〕編寫的相應的函數和過程。掌握了任何壹種計算機高級語言編程技術和數據結構知識的人,都能很方便地轉換成能上機運行的程序,有很大的靈活性與很強的實用性。
2 產生原理
假定連續隨機變量X是連續隨機變量R的函數
x= (r) (1)
由概率統計知識可知,X與R在對應點的分布函數的數值應當相等,即產生R的反變換
F(x)=F(r) (2)
其中,r為(1)式的解,即r= -1(x)=Ψ(x)
在特殊情況下,若R是在〔0,1〕區間內的均勻分布的連續隨機變量,那麽R<r的概率P(R<r)=r。此時,(2)式可以簡化為
F(x)=F(r)=r (3)
其反函數
x=F-1(r) (4)
具有分布函數F,現予以證明:因為分布函數是單調遞增函數,其反函數存在。由概率知識可知,對任意實數X,有
P(X<x)=P(F-1(R)<x)=P(R<F(x))
因為R是在〔0,1〕區間上的均勻分布的連續隨機變量,故P(R<F(x))=F(x),即F-1(r)具有分布函數F。
(4)式表明,要產生壹個服從某種分布的隨機數,可以先求出其分布函數的反函數的解析式,再將壹個在〔0,1〕區間內的均勻隨機數的值代入其中計算就可以了。
3 應用實例
例1 產生密度函數滿足指數分布
的隨機數,其中α>0。
因為 ,故分布函數為:
由(4)式可得
(5)
當r為〔0,1〕區間內的均勻隨機數時,1-r也是〔0,1〕區間內的均勻隨機數,故(5)式還可以簡化為:
(6)
假定計算機高級語言已提供了產生在〔0,1〕區間內均勻分布的隨機數的函數RND(0),則根據(6)式可以寫出產生服從指數分布的隨機數的類PASCAL語言的函數如下:
FUNC get_exp_random(alpha:real):real;
{產生服從指數分布的隨機數,值參alpha為比例參數}
get_exp_random=-ln(RND(0))/alpha;
ENDF; {get_exp_random}
例2 產生服從正態分布的隨機數
正態分布X~N(a,σ)都可以由標準的正態分布X′~N(0,1)經過變換X=a+σX′得到,在此只討論服從標準正態分布的隨機數的產生方法。標準正態分布的密度函數為:
(7)
但其分布函數不能用有限的解析形式來表示,給轉換帶來困難,但可以用以下的變換來產生隨機數。
假定r1與r2是〔0,1〕區間的兩個獨立的均勻分布隨機數,現將其作如下的變換,令
(8)
再將其反變換為:
(9)
(9)式中的C為常數。由此可導出其密度函數為:
(10)
比較(10)式與(7)式,可知X1與X2是兩個相互獨立的服從正態分布的隨機變量,因而(8)式是與(4)相當的據以產生隨機數的表達式。由此可以寫出產生服從標準正態分布隨機數的類PASCAL語言的過程如下:
PROC gauss(VAR x1,x2:real);
{產生服從正態分布的隨機數對,用變量參數x1與x2傳遞隨機數}
pi:=3.14159265; {賦π值}
r1:=RND(0); r2:=RND(0);
s:=SQRT(-2*ln(r1)); {中間變量賦值}
x1:=s*cos(2*pi*r2);
x2:=s*sin(2*pi*r2);
{正態分布隨機數賦給變量參數}
ENDP; {gauss}
此過程每調用壹次能產生兩個相互正交的標準正態分布隨機數。