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

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

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

 

Code Snippet
  1. double prob_fitness(double x)
  2. {
  3.     if(x < 0.0 ) return prob_fitness(-0.5 * x);
  4.     else if(x < 10.0) return 10.0 ;
  5.     else if(x < 20.0) return x - 10.0;
  6.     else if(x < 30.0) return 2 * x;
  7.     else if(x < 40.0) return 0.5 * x * x - 20.0 * x + 14;
  8.     else if(x < 50.0) return sqrt(x) + x + 1.2;
  9.     else return prob_fitness( x * 0.000001);
  10. }

 

只要能將機率模型表達出來即可。另較佳的方式應為  反函式法  (也不只這方式),有興趣可自行閱讀。這裡講的方法較像是蒙地卡羅測試,筆者不確定到底有沒有正式的名稱 ( 看起來是有點像 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
  1. double RandomGenerate(
  2.     double x_max, double x_min,
  3.     double y_max, double y_min,
  4.     double (*fx)(double))
  5. {
  6.     double a, b, xa, yb;
  7.     do{
  8.         a = (double)rand() / RAND_MAX;
  9.         b = (double)rand() / RAND_MAX;
  10.         xa = a * (x_max - x_min) + x_min;
  11.         yb = b * (y_max - y_min) + y_min;
  12.     }while( yb > fx(a));
  13.     return xa;
  14. }

 

 

 

edisonx 發表在 痞客邦 PIXNET 留言(0) 人氣()