產生任意機率分布的亂數的公式推導和驗證?

Inverse distribution transform:

設 F(x) 為一連續型單變量 distribution function, 其 support 為

一區間 [a,b]. 則在此區間 F 有 inverse, 以 G 表示之. 即:

對任意 x in [a,b], G(F(x)) = x;

對任意 p in [0,1], F(G(p)) = p.

今若 U 是 [0,1] 間 uniform distributed 的隨機變數, 則 Y = G(U)

具有 distribution function F.

G 為 F 之 inverse, 也就是 p = F(x) iff. x = G(p), 或說是 p 和

x 之間具有 p = F(x) = ∫_a^x dF(t) = ∫_a^x f(t) dt 的關係,其

中 f 為 F 之 density function.

因此, 欲產生 distribution function 為 F 之 pseudo random number,

可考慮先產生 [0,1) 之 uniform PRN, 再做 Y = G(U) 之轉換. 沒有

inverse function G 之 explicit form 時, 就是解方程式 F(x) = U.

用 p.d.f. 表示, 就是解 ∫_a^x f(t) dt = U. (x 為方程式之 unknown.)

但有時我們只有分布的 "形" 而不是 distribution function 或 p.d.f.,

也就是 ∫_R f(x) dx ≠ 1. 但基於機率的基本特性, 在

∫_R f(x) dx 為 finite 的條件下, 我們知道 f(x)/∫_R f(t) dt

就是所要的 p.d.f.. 所以, 上列方程式就成為

∫_a^x f(t) dt / ∫_a^b f(t) dt = p

也就是

∫_a^x f(t) dt = p * ∫_a^b f(t) dt

這也就是所引 BBS 文章中寫

integral[a,p]f(x)dx=u*integral[a,b]f(x)dx

的緣故.

如 p.d.f. 是 f(x) = 2(x-a)/(b-a)^2, 則 d.f. 是

F(x) = (x-a)^2/(b-a)^2 in a < x < b

而 F(x) = p 即 x = a + √[p(b-a)^2].

若 p.d.f. 是 f(x) = 2x/(b^2-a^2), 則 d.f. 是

F(x) = (x^2-a^2)/(b^2-a^2) in 0 ≦ a < x < b.

此時解 F(x) = p 得到的是 x = √[a^2+p(b^2-a^2)].

以 p.d.f. 的 "形" 來說, 前者是 f*(x) = x-a in a < x < b,

後者是 f*(x) = x in 0 ≦ a < x < b. 除非 a = 0, 否則二者有一

個垂直差距 a 的平移. 或者說, 前者有 x-intersection a, 而後者

卻通過原點.

合二者為一 general form:

f*(x) = x - h, h ≦ a < x < b,

則 ∫_a^b f*(x) dx = [(b-h)^2-(a-h)^2]/2, 故 d.f. 為

F(x) = [(x-h)^2-(a-h)^2]/[(b-h)^2-(a-h)^2]/2

而 p.d.f. 為

f(x) = 2(x-h)/[(b-h)^2-(a-h)^2].

方程式 F(x) = p 之解為

x = h + √{ (a-h)^2 + p [(b-h)^2-(a-h)^2]}

而如果 U 是一個 [0,1) uniform distributed 的 random variable,

則 Y = h + √{ (a-h)^2 + U [(b-h)^2-(a-h)^2]} 就是一個具有上

述 distribution function F 或 p.d.f. f 之 r.v..