1. 程式人生 > >【編程珠璣】【第一章】生成隨機數、隨機取樣的問題

【編程珠璣】【第一章】生成隨機數、隨機取樣的問題

當前 rand 可用 生成 奇數 sel 浪費 print 運行時

一、利用隨機數函數生成隨機數

問題1(《編程珠璣》習題12.1後半段):

給定一個rand(),可以產生從0RAND_MAX的隨機數,其中RAND_MAX很大(常見值:16int能表示的最大整數32767),寫出利用rand()生成[a,b]中任意整數的函數,其中a>=0, b<=RAND_MAX,且b-a<<RAND_MAX。

分析:

這是在編程工作最常見的隨機函數的應用,在這裏做一個起點再合適不過。把隨機數區間的起點從0變為a,同時把一共RAND_MAX+1個數的區間縮小至只含有b-a+1個數的區間,寫為 a + rand()%(b-a+1),此時顯然最大值是a+(b-a)= b

進一步地,這個b-a<<RAND_MAX的條件雖然看上去不起眼,其實很重要。

附加思考:

如果b-aRAND_MAX很接近會發生什麽情況?讀者不妨先做思考,問題2的分析會做出解答。這個rand(),其實相當於《編程珠璣》提到的bigrand()

問題2(筆試題變形):

給定一個隨機數函數rand7(),它能以等概率生成1~77個整數。請根據rand7()寫出類似的rand5()

分析:

如果直接像問題1中一樣,把1+rand7()%5作為rand5()會有什麽情況發生?這時確實能產生1~5的隨機數沒錯,可是各個數的概率相等嗎?

對於隨機數2,既有可能來自於1+1%5,也有可能來自於

1+6%5,顯然其概率是2/7,而不是1/5,不滿足rand5()等概率產生各隨機數這一隱含要求。不同於問題1,問題1中一個很大的區間收縮成較小區間時,各個元素映射後的新元素概率雖然概率可能不完全一樣,但卻是近似相同的。

為了滿足等概率的要求,可以這麽做:

int rand5() {
    int res;
    do {
        res = rand7();
    } while(t>=6);
    return res;
}

雖然保證了1~5的概率都變成了1/5,但是有一個無法避免的缺點是,每當產生了6或者7都要拋棄,相當於這一次運行是“空轉”,浪費了時間。如果對1/5這個概率不明白,可以有兩種理解:每次產生

67就被拋棄,剩下數的概率相等,必然為1/5;或者用更嚴密的推理:產生1~5的隨機數,最終得到某一個的概率為:1/7+(1/7)*(2/7)+(1/7)*(2/7)2+...,無限項求和,結果是1/5

問題3(筆試題原題)

給定一個隨機數函數rand7(),它能以等概率生成1~77個整數。請根據rand7()寫出類似的rand10()

分析:

有了問題2的概率基礎,把rand7()變成rand10()僅僅需要一點點思考了。

int rand10() {
    int t1,t2,res;
    do {
        t1 = rand7();      //t1以等概率1/5成為{1,2,3,4,5}中的一個數字
    } while(t1>=6);
   
    do {
        t2 = rand7();    //t2以1/2概率成為奇數或偶數,因為{1,2,3,4,5,6]中3奇3偶
    } while(t2==7); 
    res = t1+5*(t2%2);       //res是1~10中的任意一個數的概率都是1/10,註意到%和*具有相同的優先級,這裏去掉括號結果是錯的    
return res;
}

解釋:

t1的概率是{1~5}中每個元素概率為1/5t21/2概率為奇或者偶,所以(t2%2)是以1/2概率選擇0或者1,因此5*(t2%2)是以1/2概率為0或者5。這樣,選擇1的概率是1/5 * 1/2=1/10...,選擇6的概率是1/5 * 1/21/10....由此可見。

問題2和問題3是對《編程珠璣》上使用範圍不大的randint(i,j)生成其他範圍隨機數方法的解答。在掌握了問題2和問題3的解法後,你已經學會隨機數區間的收縮和擴張,類似的問題迎刃而解。

問題4(《編程珠璣》習題12.1前半段)

C庫函數rand()常返回15個隨機位,寫出bigrand(),能夠返回30個隨機位。

分析:

其實問題4和問題3有點像,但是不同之處在於,這次我們的視角是從位出發的,把rand()看做了將15個位每一個位以1/2概率設為01,從而生成0~RAND_MAX。從某種意義上來說,按這種理解方式來解這個問題更輕松一些,但是僅局限於2的冪減1這樣的數值的區間,比如從011...11。在此基礎之上把兩部分的位拼接起來即可。

//《編程珠璣》的答案
//int bigrand() {return RAND_MAX *rand() +rand();}
//最大值不是30個1,在另一個筆記中有敘述,懷疑有錯我的答案,把先生成的部分左移15位作為高15位
long bigrand(){return ((long)RAND_MAX+1)*rand() +rand();}

擴展:拒絕采樣

這是一個簡單粗暴有效的方法,使用這個方法你可以不用考慮復雜的區間擴展和收縮的問題。以rand7()產生rand10()為例,構造一個二維數組並把它填成:

技術分享圖片

然後使用兩次rand7()分別生成行號和列號,如果對應元素是0,拋棄重來;否則即是結果。其實這種方法依然用到了區間擴展:把7擴展成7*7,並把不符合的部分拋棄。這樣橫向生成一個1~7的隨機數,縱坐標生成1~7的隨機數,這樣兩次隨機數便能定位一個數組元素,若該元素為1~10之間的數,則被選中,否則重新選擇。這能保證1~10之間的每個數被選中的概率都是1/10從這裏可以看出這種方法其實還是有缺陷的:如果用rand7()生成rand50(),那這50個數可是填滿這個表格後還有填不下的。

二、利用隨機數函數產生隨機事件

用隨機事件表示隨機事件?這個問題具體化為:使用rand()表示以M/N的概率發生的隨機事件,其中M<=N,並可用作:if(事件A發生) ,其中P(A) = M/N,那麽表示為:

if(rand()%N< M)

...

通過這種方式,我們就可以做出讓程序“以M/N的概率執行某個命令”這樣的設計了。

三、取樣問題:n個元素中隨機選取m

1.從概率角度出發

考慮整數0,1,2,...,n-1,可以用上節的方法以概率m/n選取0(推導方式略)。但是對於1,必須考慮之前0是否被選取而以(m-1)/(n-1)m/(n-1)的概率選取1,後續就更加麻煩。好在叠代是計算機的長項,只需要把這個是否選擇當前數的隨機事件稍作修改即可,使之變成從r個其余的數中選擇s個:

int gen(int m,int n){
    int i, select = m,remaining = n;
    for(i=0;i<n;i++) {
        if(rand() % remaining <select) {
            printf("%d\n",i);
            select--;
        }
        remaining--;
    }
    return 0;
}

其概率證明可見於KnuthThe Art of Computer Programming2卷。進一步地可以優化為

int genknuth(int m,int n){
    int i;
    for(i=0;i<n;i++)
        if(rand()%(n -i) < m) {
            printf("%d\n",i);
            m--;
        }
    return 0;
}

《編程珠璣》提示,這種算法是“所有n的所有m元子集被選中的概率相等”,這個條件強於“所有元素被選中的概率相同”。下面是習題12.2中提到的“以等概率選擇有元素,但是有些m元子集被選中的概率比其他子集大”的算法:直接選擇1個數,則這個m元集合為它本身即後續的一共m個數,可能包括回繞。

對於這種方法,總要產生n次隨機數。進一步可以寫為for(i=0;i<n&&m>0;i++),但程序的平均運行時間是否變得更快需要權衡。對應地,習題12.7提供了一種稍微減少隨機數產生次數的遞歸函數:

int randselect(int m,int n) {
    int r;
    //assert(m<=n && m>=0);
    if(m>0) {
        r = rand()%n;
        if(r < m) {
            printf("%d\n",n-1);
            randselect(m-1,n-1);
        }
        else
            randselect(m,n-1);
    }
    return 0;
}

2.從集合插入出發

由於集合元素不重復,如果按等概率選擇一個隨機數,不在集合中就把它插入,反之直接拋棄,直到集合元素個數達到m個,同樣可以滿足要求,並且用C++STL很容易實現:

void gensets(int m,int n) {
    set<int> S;
    while(S.size() < m)
        S.insert(rand()%n);
    set<int>::iterator i;
    for(i = S.begin();i!=S.end();++i)
        cout<<*i<<"\n";
}

這個算法的主要問題是,如果拋棄已存在的元素的次數過多,相當於多次產生隨機數並進行集合操作,性能將明顯下降。比如當n=100m=99,取第99個元素時,算法“閉著眼睛亂猜整數,直到偶然碰上正確的那個為止”(《編程珠璣(續)》,13.1節)。雖然這種情況會在“從一般到特殊”提供解決方案,但下面的Floyd算法明顯規避了產生隨機數超過m次的問題。

習題12.9提供了一種基於STL集合的隨機數取樣方法,可以在最壞情況下也只產生m個隨機數:限定當前從中取值的區間的大小,每當產生重復的隨機數,就把這一次叠代時不會產生的第一個隨機數拿來替換。

int genfloyd(int m,int n){
    set<int> S;
    set<int>::iterator i;
    for(int j = n-m; j<n;j++) {
        int t = rand()%(j+1);
        if(S.find(t) == S.end())
            S.insert(t);
        else
            S.insert(j);
    }
    for(i=S.begin();i!=S.end();++i)
        cout<<*i<<"\n";
}

不必基於集合而實現,我自己寫了一個原理類似的算法,思想是把“用過”的元素“扔”到下一次叠代的隨機數取樣區間之外:

int gen(int m,int n) {
    int *array;
    int i,j;
    array = (int *)malloc(sizeof(int) * n);
    for(i=0;i<n;i++)
        array[i] = i;
    while(m>=1) {
        j = rand()%n;
        printf("%d\n",array[j]);
        if(j<n-1)
            array[j] = array[n-1];
        m--;
        n--;
    }
    return 0;
} 

3“打亂順序”出發

這是個來源於實際的想法:將所有n個元素打亂,取出前m個。更快的做法是,打亂前m個即可。對應的C++代碼如下:

int genshuf(int m,int n){
    int i,j;
    int *x = new int[n];
    for(i = 0;i<n;i++)
        x[i] = i;
    for(i = 0;i<m;i++) {
        j = randint(i,n-1);  //randint產生i到n-1之間的隨機數
        int t = x[i];x[i] = x[j];x[j] = t;
    }
    //sort(x,x+m);            //sort是為了按序輸出 
    for(i=0;i<m;i++)
        cout<<x[i]<<"\n";
}       

sort註釋掉的這段代碼,可以作為隨機不重復序列產生器。類似的還有Floyd的算法P。(《編程珠璣(續)》,13.3節)

4從一般到特殊

以上討論的幾種方式都不限定mn的取值,只需m<=n即可。對於特殊的取值,有特殊的解決方案,以下是編程珠璣上的兩例:

  1.n = 106m=n-10,這時m很大較為接近n可以先生成10個元素,然後輸出其余的元素。進行這種處理的額外代碼可以提高算法的平均速度。

  2.n=231m = 107m<<n,這時可以先生成1.1*107個元素,排序後去掉重復的(由於n很大,m中出現重復元素的概率很低),得到 107個元素的樣本。

附:

(《編程珠璣(續)》習題13.5Doug McIlory的求N個元素中取M種的第G種情況的組合的算法,原書這個算法我沒有理解,也沒有看到比較滿意的解釋,可能用到了某些我不熟悉的組合數性質。原書中沒有對其進行進一步解釋,並且似乎只是作者的題外話而已。如果看著有困難,這部分代碼完全可以跳過。介紹這個算法的原因是用它能在獲得隨機數G後,直接獲得第GN個元素取M個的取法,相當於只產生了一次隨機數。

int comb(int n,int m,int g,int* array){
    int d = 1,t;
    while(m>0) {
        t = combination(n-d,m-1);
        //combination(n,m)是計算n中取m個的組合數
        if(g-t<0) {
            m--;
            array[m] = d;
            printf("%d\n",d);
        }
        else
            g-=t;
        d++;
    }
    return 0;
}

三、取樣問題:從未知總數的元素中選擇一個

從事先未知總數為n的對象中隨機選擇一個的方法。有兩種常見的具體問題:

  1.讀取一個未知行數的文件,隨機輸出其中的的一行,同時最多只能緩沖一行的內容(《編程珠璣(續)》習題15.3利用了這種形式);

  2.對鏈表進行一趟遍歷,隨機輸出其中的一個結點的元素,只能使用一個臨時指針。

解法是,以概率1選擇第一個元素存入緩存,以概率1/2用第二個元素替換掉緩存,...,直至遍歷完所有元素,最後輸出緩存的內容。可以分析出此時所有元素留在緩存的概率均為1/n,比如1,是1*(1/2)*(2/3)*...*((n-1)/n)

int random_select(void){
    int i,num=1;
    for(i=1;i<n;i++)
       if(rand()%i ==0)     //等於0或者任何一個0~i-1之間的數,表示概率1/i。
           num = i;         //以1/i的概率替換num的緩存值。
    return num;
}

對於這個問題的常見應用:馬爾科夫文本生成器裏選擇哈希表中某項對應鏈表的任意一個結點。

【編程珠璣】【第一章】生成隨機數、隨機取樣的問題