[C&++] 任意亂數分佈產生器 (?)－Edison.X. Blog

在學數值分析時，會講到如何以均勻分佈，轉換到各種分佈裡，常見的大概會有波松分佈、指數分佈、常態分佈。較有名的方式 Miller-Box Transform 是專用在轉常態分佈 ( 因常態分佈是算常出現的機率模型 )。

一些常見的機率分佈亂數，C++ 都已包到 TR1 裡去，有興趣可自行參考。

這裡提的轉換方式，當然要有一份均勻的亂數產生器 ( 統計用的均勻就夠了，如 rand() )，另外假設讀者也知道該亂數分佈的機率函式 f(x) 為何，並可用程式語言表達出該 f(x) ，其 f(x) 長得奇型怪狀也沒差，如

Code Snippet
double prob_fitness(double x)
{
    if(x < 0.0 ) return prob_fitness(-0.5 * x);
    else if(x < 10.0) return 10.0 ;
    else if(x < 20.0) return x - 10.0;
    else if(x < 30.0) return 2 * x;
    else if(x < 40.0) return 0.5 * x * x - 20.0 * x + 14;
    else if(x < 50.0) return sqrt(x) + x + 1.2;
    else return prob_fitness( x * 0.000001);
}

只要能將機率模型表達出來即可。另較佳的方式應為反函式法 (也不只這方式)，有興趣可自行閱讀。這裡講的方法較像是蒙地卡羅測試，筆者不確定到底有沒有正式的名稱 ( 看起來是有點像 Acceptance/Rejection Method )。

首先，必須先確定，定義域之範圍 [x_min , x_max]，值域範圍 [y_min, y_max]，再令機率模型為 f(x)，步驟如下。

(1) 產生 [0,1] 浮點數亂數 a

(2) 產生 [0,1] 浮點數亂數 b

(3) xa = a * (x_max - x_min) + x_min;

(4) yb = b * (y_max - y_min) + y_min;

(5) yb <= f(xa) ，成立，取出 xa；不成立，回到步驟 1 。

原理可再細思，與蒙地卡羅法求 PI 有幾分相似，程式碼略如下述。

Code Snippet
double RandomGenerate(
    double x_max, double x_min,
    double y_max, double y_min,
    double (*fx)(double))
{
    double a, b, xa, yb;
    do{
        a = (double)rand() / RAND_MAX;
        b = (double)rand() / RAND_MAX;
        xa = a * (x_max - x_min) + x_min;
        yb = b * (y_max - y_min) + y_min;
    }while( yb > fx(a));
    return xa;
}