本節主要受到《程式設計珠璣》第12章隨機取樣問題的啟發，但不僅僅限於隨機取樣問題，進一步地，研究討論了一些在筆試面試中常見的和隨機函式以及概率相關的問題。

　　閱讀本文所需的知識：

　　　　1.對C語言中或其他語言中等價的rand()、srand()有所瞭解。本文不討論種子的設定和偽隨機數的問題；

　　　　2.中學或以上水平的概率基本概念。

　　目錄

利用隨機數函式生成隨機數

問題1（《程式設計珠璣》習題12.1後半段）：

給定一個rand()，可以產生從0到RAND_MAX的隨機數，其中RAND_MAX很大（常見值:16位int能表示的最大整數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-a和RAND_MAX很接近會發生什麼情況？讀者不妨先做思考，問題2的分析會做出解答。

這個rand()，其實相當於《程式設計珠璣》提到的bigrand()。

問題2（筆試題變形）：

　　給定一個隨機數函式rand7()，它能以等概率生成1~7這7個整數。請根據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這個概率不明白，可以有兩種理解：每次產生6或7就被拋棄，剩下數的概率相等，必然為1/5；或者用更嚴密的推理：產生1~5的隨機數，最終得到某一個的概率為：1/7+(1/7)*(2/7)+(1/7)*(2/7)²+...，無限項求和，結果是1/5。

問題3（筆試題原題）

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

分析：

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

int rand10() {
    int t1,t2,res;
    do {
        t1 = rand7();
    } while(t1>=6);
    //t1以等概率成為1/5

    do {
        t2 = rand7();
    } while(t2==7); 
    //t2以1/2概率成為奇數或偶數

    res = t1+5*(t2%2);
    //res是1~10中的任意一個數的概率都是1/10
    //注意到%和*具有相同的優先順序，這裡去掉括號結果是錯的    
    return res;
}

等概率的rand10()

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

問題4（《程式設計珠璣》習題12.1前半段）

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

分析：

其實問題4和問題3有點像，但是不同之處在於，這次我們的視角是從位出發的，把rand()看做了將15個位每一個位以1/2概率設為0或1，從而生成0~RAND_MAX。從某種意義上來說，按這種理解方式來解這個問題更輕鬆一些，但是僅侷限於2的冪減1這樣的數值的區間，比如從0到11...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，並把不符合的部分拋棄。從這裡可以看出這種方法其實還是有缺陷的：如果用rand7()生成rand50()，那這50個數可是填滿這個表格後還有填不下的。

利用隨機數函式產生隨機事件

　　用隨機事件表示隨機事件？經過上面的區間收縮的思考，看上去並不難。把這個問題具體化為：使用rand()表示以M/N的概率發生的隨機事件，M<=N，並可用作：if(事件A發生)  ，其中P(A) = M/N，那麼表示為：

if(rand()%N< M)
    ...

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

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

從概率角度出發

考慮整數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;
}

　　其概率證明可見於Knuth的The Art of Computer Programming第2卷。進一步地可以優化為：

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;
}

 從集合插入出發

由於集合元素不重複，如果按等概率選擇一個隨機數，不在集合中就把它插入，反之直接拋棄，直到集合元素個數達到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=100而m=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;
}

從“打亂順序”出發

這是個來源於實際的想法：將所有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節）

從一般到特殊

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

　　1.n = 10⁶而m=n-10，這時可以先生成10個元素，然後輸出其餘的元素。進行這種處理的額外程式碼可以提高演算法的平均速度。

　　2.n=2³¹而m = 10⁷，m<<n，這時可以先生成1.1*10⁷個元素，排序後去掉重複的（由於n很大，m中出現重複元素的概率很低），得到 10⁷個元素的樣本。

附：

　　（《程式設計珠璣（續）》習題13.5）Doug McIlory的求N個元素中取M種的第G種情況的組合的演算法，原書這個演算法我沒有理解，也沒有看到比較滿意的解釋，可能用到了某些我不熟悉的組合數性質。原書中沒有對其進行進一步解釋，並且似乎只是作者的題外話而已。如果看著有困難，這部分程式碼完全可以跳過。介紹這個演算法的原因是用它能在獲得隨機數G後，直接獲得第G種N個元素取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++)
    //i<n代表某種終止條件，n未知    
        if(rand()%i ==0)
            num = i;
    return num;
}

 　　對於這個問題的常見應用：馬爾科夫文字生成器裡選擇雜湊表中某項對應連結串列的任意一個結點。

概率問題選編

1.（習題12.4）對於集合插入問題，每次產生隨機數後都需要在集合中檢測該元素是否已經存在，這個測試次數和呼叫隨機數函式的次數相同。對於從n個元素中選擇m個元素的問題，平均需要測試多少次才能保證選擇了m個元素？

這裡先處理特殊的問題，即“贈券收集問題”：必須收集多少張棒球卡才能保證擁有所有的n種卡？

記P_i為從收集了i-1種到i種的概率，P_i=(n-i+1)/n

此時需要收集的卡片數目n_i = 1* P_i + 2*(1-P_i)*P_i + 3*(1-P_i)²*P_{i +} ... + n_*(1-P_i)ⁿ*P_{i + ...}，這個無窮級數可以求解為n_i=1/P_i

那麼所需要總的卡片數目sum(n_i)  = n₁+n₂ +... +n_n約為nlnn + γn + O(1)。相關維基百科

對於從n個元素中選擇m個，相應地sum(n_i)  = n₁+n₂ +... +n_m，約為nlnm+γn。（根據調和數的推導）

2.（習題12.11）每個玩家有一張包含16個覆蓋點的紙牌，覆蓋點下面隱藏著1~16的隨機排列。玩家每次刮開一個點，如果出現3，則判玩家負；如果出現1或2，則判玩家勝。那麼，隨機選擇覆蓋點刮開，獲勝的概率是多少？

解答：

其實勝負只對應兩種排列：1和2都出現在3前則勝，否則負。前者的即*1*2*3*或*2*1*3*，後者為其餘四種排列。顯然獲勝的概率為2/3。

不要把問題複雜化，應該儘量簡化。4~16的數字都是沒有用的。