2. Phát Sinh Các Số (giả) Ngẫu Nhiên | Blog Của Chiến

Trang chủ » Số Giả Ngẫu Nhiên Là Gì » 2. Phát Sinh Các Số (giả) Ngẫu Nhiên | Blog Của Chiến

Có thể bạn quan tâm

§2 Phát sinh các số (giả) ngẫu nhiên

Trở về Mục lục cuốn sách

Mã nguồn người dịch mới upload: randomness.zip

Nhiều trường hợp mô phỏng trong khoa học, kinh tế hay khoa học xã hội cần đến biến ngẫu nhiên. Thường thì mô hình tự nó bộc lộ những tham số ngẫu nhiên mà vẫn được giữ cố định trong suốt quá trình mô phỏng; ta nói đến quenched disorder (một dạng phi trật tự). Một ví dụ nổi tiếng trong lĩnh vực vật lý thể đặc là thủy tinh spin, vốn là hợp kim trộn ngẫu nhiên các vật liệu từ tính và phi từ tính. Trong trường hợp này, khi thực hiện những mô phỏng đối với hệ thống nhỏ, để thu được đại lượng vật lý cần thiết, ta phải tiến hành lấy trung bình các biểu hiện nhiễu loạn khoác nhau. Mỗi biểu hiện nhiễu loạn bao gồm các vị trí của hạt từ tính và phi từ tính được chọn một cách ngẫu nhiên. Để phát sinh ra biểu hiện nhiễu loạn phục vụ mô phỏng, ta cần có các số ngẫu nhiên. Nhưng ngay cả khi hệ được mô phỏng bản thân không có tính ngẫu nhiên thì thuật toán dùng để tính lại thường yêu cầu số ngẫu nhiên, chẳng hạn để lập nên một tập hợp thống kê (ensemble) chứa những nhiệt độ hữu hạn, hoặc khi dùng đến thuật toán ngẫu nhiên. Tóm lại, ứng dụng của số ngẫu nhiên trong mô phỏng điện toán rất phổ biến. Trong mục này, chúng tôi trình bày sự phát sinh số ngẫu nhiên. Trước hết là phần giải thích cách phát sinh chúng bằng mọi cách trên máy tính. Sau đó, các phương pháp khác nhau sẽ được trình bày nhằm mục đích thu được số tuân theo một dạng phân bố mong muốn: phương pháp nghịch đảo, phương pháp loại bỏ, và phương pháp Box-Müller. Các thông tin đầy đủ về những phương pháp này và tương tự có thể được tìm thấy trong các tài liệu tham khảo: [Morgan (1984); Devroye (1986); Press và nnk. (1995)]. Trong mục này tôi coi rằng bạn đã quen thuộc với những khái niệm cơ bản về lý thuyết xác suất và thống kê, như đã trình bày trong Mục 1.

2.1 Các số (giả) ngẫu nhiên phân bố đều

Trước hết, cần chỉ ra rằng máy tính thông dụng không có tính tất định. Do vậy, nó không có khả năng trực tiếp phát sinh những số ngẫu nhiên thực sự. Song ta có thể tạo tương tác với người dùng. Chẳng hạn, ta có thể đo độ dài khoảng thời gian giữa hai lần gõ phím liên tiếp, về bản chất thì thời gian đó phân bố ngẫu nhiên. Nhưng khoảng thời gina thu được lại p hụ thuộc nhiều vào người dùng hiện tại và có nghĩa rằng các thuộc tính thống kê không thể kiểm soát được. Mặt khác, có những thiết bị ngoại vi được cài đặt những quá trình ngẫu nhiên vật lý, chúng có thể được gắn vào máy tính [Quantis; Westphal] hoặc sử dụng thông qua internet [Hotbits]. Dù sao, vì những con số này thực sự ngẫu nhiên, chúng không cho phép thực hiện mô phỏng ngẫu nhiên (stochastic simulation) theo một cách kiểm soát được và tái hiện được. Điều này rất quan trọng đối với khoa học, vì những kết quả đẹp mắt hoặc bất ngờ đều thường được cố gắng tái hiện bởi những nhóm nhà khoa học khác. Ngoài ra, có những lỗi chương trình có thể xuất hiện chỉ với những số ngẫu nhiên nhất định. Vì vậy, với mục tiêu gỡ lỗi điều quan trọng là phải chạy lại được đúng chương trình mà ta đã mô phỏng từ trước.

Một điểm khác là, với những số ngẫu nhiên thực sự, hoặc là tốc độ phát sinh số ngẫu nhiên sẽ bị hạn chế nếu việc phát sinh số ngẫu nhiên thực sự không tốn tài nguyên máy tính lắm, còn nếu không thì các giải thuật phát sinh sẽ rất tốn kém. Đó là lý do tại sao các số giả ngẫu nhiên thường được dùng. Chúng được phát sinh từ những quy tắc tất định. Cơ chế này tạo ra hàm rand() để phát sinh số ngẫu nhiên tuân theo phân bố chuẩn. Mỗi lần rand() được gọi đến, một số (giả) ngẫu nhiên mới được trả về. (Để cho tiện, từ giờ chữ “giả” sẽ được bỏ bớt.) Những số ngẫu nhiên này phải “trông giống” số ngẫu nhiên thực sự và có nhiều thuộc tính của chúng. Ta nói chúng phải là số ngẫu nhiên “tốt”. Bây giờ chúng tôi sẽ nói cụ thể “trông giống” và “tốt” nghĩa là thế nào. Ta muốn có thuật toán sinh số ngẫu nhiên sao cho mỗi số nếu có thể được tạo ra thì sẽ xuất hiện với cùng xác suất. Ngoài ra, nếu hai số được tạo ra là ri, rk chỉ hơi khác nhau thì các số ngẫu nhiên ri + 1, rk + 1 được trả về từ lần gọi hàm tiếp theo sẽ khác nhau rõ rệt, vì vậy các số liên tiếp phải có mức tương quan yếu. Có nhiều cách để chỉ định mối tương quan, vì vậy không có một tiêu chuẩn duy nhất nào. Dưới đây, cách đơn giản nhất sẽ được thảo luận. Những phương pháp đơn giản nhất để phát sinh số giả ngẫu nhiên là phát sinh đồng dư tuyến tính. Theo đó một dãy số nguyên x1, x2, … có giá trị nằm giữa 0 và m − 1 được tạo ra bằng công thức đệ quy:

xn + 1 = (axn + c)mod m .      (47)

Giá trị ban đầu x0 được gọi là nhân. Sau đây chúng tôi trình bày một đoạn mã lệnh C, lin_con(), để thực hiện phương pháp trên. Nó lưu lại số hiện thời vào biến địa phương x vốn được khai báo là static, để đảm bảo được nhớ lại, ngay cả khi hàm đã kết thúc. Có hai đối số. Đôi số thứ nhất set_seed chỉ định rằng liệu ta có muốn đặt giá trị cho nhân không. Nếu có, thì giá trị mới của nhân cần được truyền vào đối số thứ hai, còn nếu không thì giá trị của đối số thứ hai sẽ được bỏ qua. Hàm này trả về nhân nếu nó bị thay đổi, hoặc số ngẫu nhiên mới. Lưu ý các hằng số a và c được định nghĩa bên trong hàm, còn số chia M được lập qua một macro RNG_MODULUS để giúp nó có thể nhìn thấy từ bên ngoài lin_con():

#define RNG_MODULUS 32768 /* modulus */ int lin_con(int set_seed, int seed) { static int x = 1000; /* so ngau nhien hien thoi */ const int a = 12351; /* he so */ const int c = 1; /* dich chuyen */ if(set_seed) /* nhan moi ? */ x = seed; else /* so ngau nhien moi ? */ x = (a*x+c) % RNG_MODULUS; return(x); }

Mã nguồn: file rng.c

Nếu máy bạn chạy hệ Windows, đừng dùng Notepad để xem mã nguồn vì sẽ không hiện đúng các dấu xuống dòng của hệ UNIX. Hãy dùng Wordpad hoặc các trình biên tập cao cấp khác.

Nếu bạn chỉ muốn thu được số ngẫu nhiên tiếp theo, bạn không cần quan tâm đến nhân. Vì vậy, để cho tiện, ta sẽ dùng rn_lin_con() để gọi lin_con() với đối số thứ nhất bằng 0:

int rand_lin_con() { return(lin_con(0,0)); }

Nếu ta muốn đặt giá trị cho nhân, để cho tiện ta cũng dùng một hàm nhỏ được tạo riêng seed_lin_con():

void srand_lin_con(int seed) { lin_con(1, seed); }

Để phát sinh các số ngẫu nhiên r phân bố trong khoảng [0, 1) ta phải chia số ngẫu nhiên hiện thời cho số chia m. Tốt nhất là nhận được các kết quả phân bố đều trên khoảng. Các số ngẫu nhiên được phát sinh từ phân bố này có thể được dùng làm đầu vào để phát sinh ra những số ngẫu nhiên tuân theo dạng phân bố khác, về cơ bản, bất kì phân bố nào cũng được. Sau đây, bạn sẽ thấy cách phát sinh số ngẫu nhiên tuân theo các dạng phân bố khác. Hàm C đơn giản sau đây phát sinh ra những số ngẫu nhiên trong khoảng [0, 1) dùng macro RNG_MODULUS định nghĩa ở trên:

double drand_lin_con() { return( (double) lin_con(0,0) / RNG_MODULUS); }

Ta phải chọn các tham số a, c, m sao cho thu được các số ngẫu nhiên “tốt”, ở đây “tốt” có nghĩa là “có ít tương quan”. Lưu ý rằng trong quá khứ đã có những mô phỏng mà kết quả bị chứng minh là sai vì việc dùng phép phát sinh số ngẫu nhiên xấu [Ferrenberg và nnk. (1992); Vattulainen và nnk.(1994)].

Ví dụ 1. Để thấy được rằng “phép phát sinh xấu” nghĩa là thế nào, hãy xét ví dụ các tham số a = 12351, c = 1, m = 215 và giá trị nhân I0 = 1000. 10000 số ngẫu nhiên được phát sinh bằng cách chia mỗi số nhận được cho m. Chúng được phân bố đều trong khoảng [0, 1). Hình 5 cho thấy sự phân bố của số ngẫu nhiên.

Hình 5. Phân bố của số ngẫu nhiên trong khoảng [0, 1) nhận được bằng cách chuyển đổi một biểu đồ tần suất sang một hàm phân bố xác suất, xem Mục 3.3. Các số ngẫu nhiên được phát sinh bằng phép đồng dư tuyến tính với các tham số a = 12351, c = 1, m = 215.

Phân bố thu được có vẻ khá phẳng, nhưng nhìn kĩ hơn ta thấy được có đôi chỗ “quy luật” nhất định. Những đặc điểm có tính quy luật này có thể được xem xét bằng cách ghi lại bộ k số ngẫu nhiên liên tiếp (xi, xi + 1, …, xi + k − 1). Nếu là phép sinh số ngẫu nhiên tốt, không bộc lộ tương quan, sẽ làm đầy không gian k một cách đều đặn. Không may là với phép phát sinh đồng dư tuyến tính, thay vì những điểm nằm trên mặt phẳng (k − 1)-chiều, có thể cho thấy rằng cùng lắm chỉ có cỡ m1 / k mặt phẳng như vậy. Một phép phát sinh xấu còn có ít mặt phẳng hơn nữa. Đây là trường hợp ứng với ví dụ trước, xem phần trên của Hình 6.

Hình 6. Tương quan hai điểm xi + 1(xi) giữa các số ngẫu nhiên liên tiếp xi, xi + 1. Trường hợp trên được phát sinh bằng phép đồng dư tuyến tính với các tham số a = 12351, c = 1, m = 215, trường hợp dưới có a = 12349.

Kết quả với a = 123450 còn tệ hơn nữa: chỉ có 15 số “ngẫu nhiên” khác nhau được phát sinh (với nhân là 1000), đến đó phép tính lặp đạt đến một điểm cố định (không được chỉ ra trên hình). Nếu thay vào đó ta chọn a = 12349 thì tương quan hai điểm có dạng như nửa dưới Hình 6. Hiển nhiên là biểu hiện này ít quy luật hơn rõ rệt, nhưng mối tương quan yếu có thể trở nên dễ thấy với các bộ k lớn. Một phép phát sinh đạt được yêu cẩu của vài phép thử kinh nghiệm là a = 75 = 16807, m = 231 − 1, c = 0. Khi lập trình phát sinh này bạn phải cẩn thận vì trong quá trình tính toán, những con số sinh ra sẽ không vừa phạm vi 32 bit. Một cách lập trình khéo léo được trình bày ở [Press và nnk. (1995)]. Sau cùng, cần nhấn mạnh rằng phép phát sinh này, cũng giống như mọi phép phát sinh đồng dư tuyến tính khác, có các bit thấp kém ngẫu nhiên nhiều so với các bit cao. Vì vậy, khi bạn muốn phát sinh số nguyên trong khoảng [1,N], bạn nên viết

r = 1+(int) (N*x_n/m);

thay vì dùng toán tử modulo như với r=1+(x_n % N). Trong C tiêu chuẩn, có một hàm phát sinh số ngẫu nhiên lập sẵn đơn giản có tên rand() (xem tài liệu đi kèm), với số chia m = 215, nghĩa là rất tồi. Đa số các hệ điều hành cũng có drand48(), hàm này dựa trên m = 248 (a = , c = 11) và cũng yêu cầu phép tính số học riêng. Nó đủ dùng cho những mô phỏng không cần quá nhiều số ngẫu nhiên và không yêu cầu chất lượng ngẫu nhiên cao nhất. Vào những năm gần đây, một vài phép phát sinh số ngẫu nhiên đã được phát triển. Có những chương trình phát sinh rất tốt được phát hành tự do trong thư viện khoa học GNU. Vì vậy, bạn không cần phải tự lập trình chúng. Cho đến giờ, ta đã thấy cách phát sinh số ngẫu nhiên phân bố đều trên khoảng [0, 1). Nói chung, ta thường quan tâm đến những số ngẫu nhiên được phân bố theo một dạng phân bố xác suất cho trước với mật độ p(x) nào đó. Trong những mục tiếp theo, một số kĩ thuật phục vụ việc này sẽ được trình bày.

2.2 Các biến ngẫu nhiên rời rạc

Trong trường hợp phân bố rời rạc với số hữu hạn những khả năng xảy ra, ta có thể tạo một bảng các kết quả khả dĩ cùng các xác suất pX(xi) của chúng (i = 1, …, imax), với giả thiết rằng xi được xếp theo chiều tăng dần. Để rút ra một số, ta phải rút một số ngãu nhiên u với phân bố đều trong khoảng [0, 1) rồi lấy chỉ số j của bảng sao cho tổng sj ≡ ∑ i = 1jpX(xi) của các xác suất lớn hơn hơn u, nhưng sj − 1 ≡ ∑ i = 1j − 1pX(i) < u. Lưu ý rằng ta có thể nhanh chóng tìm kiếm trên mảng bằng cách tìm kiếm nhị phân: Mảng được liên tiếp chia đôi và mỗi lần việc tìm kiếm được tiến hành trên nửa mảng có chứa chỉ số j. Bằng cách này, việc phát sinh một số ngẫu nhiên có độ phức tạp thời gian chỉ tăng với cấp logarit theo số kết quả khả dĩ, imax. Cách này sẽ phát huy tác dụng nếu số kết quả khả dĩ là rất lớn.

Ở Bài tập (1) bạn được yêu cầu viết một hàm để lấy mẫu từ phân bố xác suất của một biến rời rạc, cụ thể là phân bố Poisson. Trong mục tiếp theo ta sẽ tập trung vào kĩ thuật phát sinh các biến ngẫu nhiên liên tục.

2.3 Phương pháp nghịch đảo

Cho trước một hàm sinh số ngẫu nhiên drand(); giả sử rằng nó có thể phát sinh số ngẫu nhiên U phân bố đều trong khoảng [0, 1). Mục đích của ta là phát sinh số ngẫu nhiên Z có theo mật độ xác suất pZ(z). Hàm phân bố tương ứng là

F_Z(z)\equiv \text{Pr}(Z\le z) \equiv \int_{-\infty}^z dz^{\prime} p_Z(z^{\prime}) \label{eq:pra:transform}     (48)

Mục tiêu là tìm một hàm g(u), sao cho sau phép đổi biến Z = g(U) các kết quả của Z sẽ được phân bố tuân theo (48). Giả sử rằng g có thể nghịch đảo và có tính đơn điệu tăng. Khi đó ta có

F_Z(z)=\text{Pr}(Z\le z)=\text{Pr}(g(U)\le z) = \text{Pr}(U\le g^{-1}(z))     (49)

Vì hàm phân bố F_U(u)=\text{Pr}(U\le u) của một biến phân bố đều chỉ là FU(u) = u (u ∈ [0, 1]) nên ta thu được FZ(z) = g − 1(z). Như vậy, ta phải chọn g(z) = FZ − 1(z) cho hàm biến đổi để có thể thu được các số ngẫu nhiên tuân theo phân bố xác suất FZ(z). Dĩ nhiên, cách này chỉ hoạt động khi FZ nghịch đảo được. Nếu không, bạn có thể dùng những phương pháp sẽ được trình bày trong các chương tiếp sau, hoặc bạn có thể tự phát sinh một bảng giá trị của hàm phân bố, vốn thực chất là xấp xỉ dạng rời rạc cho hàm phân bố, rồi dùng các phương pháp phát sinh số ngẫu nhiên như trong Mục 2.2. Thậm chí cách này có thể được tinh chỉnh bằng việc dùng một xấp xỉ tuyến tính cho hàm phân bố. Ở đây, ta không đi sau vào chi tiết, mà trình bày một ví dụ trong đó hàm phân bố có thể được nghịch đảo.

Ví dụ 2. Ta hãy xét phân bố lũy thừa với tham số μ, và hàm phân bố FZ(z) = 1 − exp( − z / μ). Từ đó, ta có thể nhận được các số ngẫu nhiên Z tuân theo phân bố lũy thừa bằng cách phát sinh các số ngẫu nhiên u theo phân bố đều rồi chọn z =  − μln(1 − u).

Hình 7. Biểu đồ tần số hàm mật độ xác suất của các số ngẫu nhiên được phát sinh theo phân bố lũy thừa (μ = 1) so với hàm mật độ xác suất (đường thẳng) trong biểu đồ có thang logarit.

Mã nguồn: file expo.c. Chương trình C đơn giản sau phát sinh một số ngẫu nhiên tuân theo phân bố lũy thừa. Tham số μ của phân bố được truyền vào dưới dạng đối số. Lưu ý rằng ở dòng 4 ta dùng hàm đơn giản drand48() để phát sinh số ngẫu nhiên; hàm này đã có sẵn trong thư viện C chuẩn và chạy tốt với những ứng dụng yêu cầu vừa phải về chất lượng tính toán thống kê. Với những hàm được thiết kế tinh vi hơn, bạn đọc có thể xem thư viện khoa học GNU chẳng hạn.

Hình 7 cho thấy một biểu đồ tần số hàm mật độ xác suất với 105 số ngẫu nhiên được phát sinh bằng cách này và hàm xác suất của phân bố lũy thừa với μ = 1 trên trục y chia theo thang logarit. Chỉ ở những giá trị lớn ta mới nhìn thấy sự khác biệt giữa hai đường. Sự khác biệt được phát sinh do những dao động thống kê vì ở đây pZ(z) rất nhỏ.

2.4 Phương pháp bác bỏ

Như đã đề cập từ trước, phương pháp nghịch đảo chỉ dùng được khi hàm phân bố P có thể nghịch đảo bằng cách giải tích. Với những phân bố không thỏa mãn điều kiện này, đôi khi ta có thể giải quyết vấn đề bằng cách rút ra vài số ngẫu nhiên và kết hợp chúng lại một cách khéo léo.

Hình 8. Phương pháp bác bỏ: Các điểm (x, y) rải rác đều trên một hình chữ nhật bị giới hạn. Xác suất để y ≤ p(x) sẽ tỉ lệ với p(x).

Phương pháp bác bỏ có tác dụng với các số ngẫu nhiên trong đó hàm mật độ xác suất p(x) đặt vừa trong khung [x0, x1) × [0, ymax), nghĩa là, p(x) = 0 với x ∉ [x0, x1] và p(x) ≤ ymax. Ý tưởng cơ bản của việc phát sinh một số ngẫu nhiên tuân theo phân bố p(x) là để phát sinh những cặp số ngẫu nhiên (x, y), được phân bố đều trong [x0, x1] × [0, ymax] và chỉ chấp nhận những số x nào thỏa mãn y ≤ p(x), nghĩa là những cặp số ứng với điểm phía dưới p(x), xem Hình 8. Do vậy, xác suất rút được x thì thỉ lệ với p(x), như mong muốn.

Mã nguồn: file reject.c.

Hàm C sau thực hiện phương pháp bác bỏ, viết cho một hàm mật độ xác suất bất kì. Nó nhận vào các đối số gồm biên giới của khung y_max, x0 và x1 cùng với một con trỏ pdf đến hàm chỉ định mật độ xác suất. Nếu cần giải thích về con trỏ, bạn có thể tự đọc tài liệu C.

double reject(double y_max, double x0, double x1, double (* pdf)(double)) { int found; /* co danh dau khi tim duoc so hop le */ double x,y; /* diem ngau nhien trong [x0,x1]x[0,p_max] */ found = 0; while(!found) /* lap den khi so duoc phat sinh */ { x = x0 + (x1-x0)*drand48(); /* dong deu tren [x0,x1] */ y = y_max *drand48(); /* dong deu tren [0,p_max] */ if(y <= pdf(x)) /* chap nhan ? */ found = 1; } return(x); }

Ở các dòng 9–10, điểm ngẫu nhiên được phát sinh theo phân bố đều trong hộp. Các dòng 11–12 chứa lệnh kiểm tra xem liệu có tìm được điểm nằm phía dưới đường mật độ xác suất không. Việc tìm kiếm trong vòng lặp (các dòng 7–13) tiếp tục đến tận khi một số ngẫu nhiên được chấp nhận; nó sẽ được trả lại ở dòng 14. Phương pháp bác bỏ được áp dụng cho một hàm mật độ xác suất, có mật độ bằng 1 trong đoạn [0, 0. 5) và tăng từ 0 đến 4 trong đoạn [1, 1. 5). Ở mọi chỗ khác nó bằng không. Hàm mật độ này được viết cụ thể bằng hàm sau trong C:

double pdf(double x) { if( (x<0)|| ((x>=0.5)&&(x<1))|| (x>1.5) ) return(0.0); else if((x>=0)&&(x<0.5)) return(1.0); else return(4.0*(x-1)); }

Hàm mật độ xác suất thực nghiệm thu được như trên Hình 9.

Hình 9. Biểu đồ hàm mật độ xác suất của 105 số ngẫu nhiên phát sinh bằng phương pháp bác bỏ từ một hàm mật độ xác suất nhân tạo.

Phương pháp bác bỏ luôn áp dụng được nếu mật độ xác suất bị chặn, nhưng nó có nhược điểm là cần phát sinh nhiều số ngẫu nhiên hơn số thực tế được dùng đến: Nếu A = (x1 − x0)ymax là diện tích của khung thì trung bình ta phải phát sinh 2A số ngẫu nhiên phụ trợ để nhận được chỉ một số ngẫu nhiên trong phân bố mong muốn. Nếu điều này dẫn đến hiệu quả tính toán thấp, bạn có thể xem xét việc dùng nhiều khung cho các phần khác nhau của hàm mật độ xác suất.

2.5 Phân bố Gauss

Trong trường hợp hàm phân bố không thể nghịch đảo, và xác suất cũng không thể vừa trong một khung, thì ta cần áp dụng các phương pháp riêng khác. Xét ví dụ, một phương pháp phát sinh số ngẫu nhiên phân bố theo dạng Gauss. Các phương pháp khác và ví dụ về cách kết hợp những kĩ thuật khác nhau được tập hợp trong [Morgan (1984)]. Hàm mật độ xác suất của phân bố Gauss với trị trung bình μ và phương sai σ2 được cho bởi PT (36), xem thêm Hình 10. Ngoài phân bố đều, phân bố Gauss là dạng thường gặp nhất trong mô phỏng.

Hình 10. Phân bố Gauss với trị trung bình bằng không và bề rộng bằng đơn vị. Các điểm tròn biểu diễn hàm mật độ xác suất thu được từ 104 số rút ra bằng phương pháp Box-Müller.

Ở đây, trường hợp phân bố Gauss chuẩn hóa (μ = 0,  σ = 1) được xét đến. Nếu bạn muốn lập trường hợp tổng quát, bạn phải rút một số ngẫu nhiên z tuân theo phân bố Gauss chuẩn hóa rồi dùng biến đổi σz + μ để thu được phân bố mong muốn. Vì phân bố Gauss trải dài trên khoảng vô hạn và vì hàm phân bố không thể nghịch đảo, nên những phương pháp đã trình bày đều không áp dụng được. Các đơn giản nhất để phát sinh số ngẫu nhiên tuân theo phân bố Gauss được dựa trên định lý giới hạn trung tâm. Định lý này phát biểu rằng bất kì tổng của K biến ngẫu nhiên phân bố độc lập Ui (với trị trung bình μ và phương sai v) sẽ hội tụ về một phân bố Gauss có trị trung bình Kμ và phương sai Kv. Nếu một lần nữa, Ui được rút từ phân bố đều trên khoảng [0, 1) (có trị trung bình μ = 0. 5 và phương sai v = 1 / 12), thì ta có thể chọn K = 12, khi đó biến ngẫu nhiên Z = ∑ i = 1KUi − 6 sẽ phân bố xấp xỉ với phân bố Gauss chuẩn hóa. Nhược điểm của phương pháp này là cần có 12 số ngẫu nhiên để phát sinh ra một số ngẫu nhiên kết quả, và những số lớn hơn 6 hoặc nhỏ hơn  − 6 sẽ không bao giờ xuất hiện. Đối lập với cách làm này, phương pháp Box-Müller lại có tính chính xác. Bạn cần hai biến ngẫu nhiên U1, U2 phân bố đều trên [0, 1) để phát sinh ra hai biến Gauss độc lập, N1, N2. Điều này có thể đạt được bằng cách phát sinh u1, u2 từ U1, U2 và gán

\begin{aligned} n_1 & = & \sqrt{-2 \log (1-u_1)} \cos(2\pi u_2) \\ n_2 & = & \sqrt{-2 \log (1-u_1)} \sin(2\pi u_2) \end{aligned}

Một cách chứng minh rằng n1 và n2 đúng là phân bố theo (36) có thể thấy ở [Press và nnk. (1995); Morgan (1984)]; tài liệu đó còn đề cập đến các phương pháp khác để phát sinh số ngẫu nhiên Gauss, kể cả những cách hiệu quả hơn. Một phương pháp dựa trên mô phỏng các hạt (chất điểm) trong hộp được giải thích trong [Fernandez và Criado (1999)]. Hình 10 là một đồ thị hàm mật độ xác suất của 104 số ngẫu nhiên rút bằng phương pháp Box-Müller. Lưu ý rằng bạn có thể tìm thấy một cách lập trình phương pháp Box-Müller trong lời giải của Bài tập (3).

Bài tập

(1) Lấy mẫu từ phân bố rời rạc Hãy thiết kế, lập trình và kiểm tra một hàm nhằm trả lại một số ngẫu nhiên tuân theo một hàm phân bố rời rạc nào đó, được lưu trong một mảng F, như đề cập đến ở Mục 2.2. Mẫu của hàm này như sau:

/******************** rand_discrete() *****************/ /** Tra lai so ngau nhien phan bo **/ /** theo dang roi rac cho boi **/ /** ham phan bo trong mang 'F' **/ /** Su dung cach tim kiem tren array de phat sinh so **/ /** PARAMETERS: (*)= tham so tra lai **/ /** n: so phan tu trong mang array **/ /** F: array voi ham phan bo. **/ /** RETURNS: **/ /** so ngau nhien **/ /******************************************************/ int rand_discrete(int n, double *F)

Để đơn giản, bạn có thể dùng hàm drand48() từ thư viện C chuẩn cho việc phát sinh các số ngẫu nhiên phân bố theo U(0, 1). Ngoài ra, hãy thiết kế, lập trình và kiểm tra một hàm để cấp phát bộ nhớ và khởi tạo cho mảng F từ một phân bố Poisson với tham số μ, xem PT (27) để biết về hàm khối xác suất. Hàm này phải tự động xác định được rằng cần có bao nhiêu phần tử của F, tùy theo tham số μ. Mẫu hàm này như sau:

/********************* init_poisson() *****************/ /** Tra lai mang voi ham phan bo **/ /** cho dang phan bo Poisson voi tri trung bin mu: **/ /** p(k)=mu^k*exp(-mu)/x! **/ /** Kich co mang duoc tu dong dieu chinh. **/ /** PARAMETERS: (*)= tham so tra lai **/ /** (*) n_p: con tro toi so phan tu trong bang. **/ /** mu: tham so dang phan bo. **/ /** RETURNS: **/ /** con tro den mang bieu dien ham phan bo. **/ /******************************************************/ double *init_poisson(int *n_p, double mu)

Gợi ý: Để xác định kích thước của mảng, đầu tiên bạn có thể lặp qua những xác suất rồi lấy giá trị đầu k_0 với p(k_0) = 0 với mức độ chính xác phù hợp. Giá trị k_0 này đóng vai trò làm kích thước của mảng. Hoặc, bạn có thể bắt đầu với một kích thước nào đó rồi mở rộng mảng bằng cách gấp đôi kích thước hiện có mỗi khi cần. Để kiểm tra, bạn có thể phát sinh nhiều số, tính trị trung bình rồi so sánh kết quả với μ. Hoặc, bạn có thể ghi lại một histogram rồi so sánh kết quả với PT (27).

Lời giải: poisson.c Hãy cố gắng tự làm bài tập trước khi xem lời giải!

(2) Phương pháp nghịch đảo với phân bố Fisher-Tippett Hãy thiết kế, lập trình và kiểm tra một hàm nhằm trả lại một số ngẫu nhiên tuân theo phân bố Fisher-Tippett, PT (43) với tham số λ. Hãy dùng phương pháp nghịch đảo Mẫu của hàm như sau:

/******************** rand_fisher_tippett() ***********/ /** Tra lai so ngau nhien von duoc phan bo theo **/ /** dang phan bo Fisher-Tippett. **/ /** PARAMETERS: (*)= tham so tra lai **/ /** lambda: tham so cua dang phan bo **/ /** RETURNS: **/ /** so ngau nhien **/ /******************************************************/ double rand_fisher_tippett(double lambda)

Nhận xét: Để đơn giản, bạn có thể dùng hàm drand48() từ thư viện C chuẩn cho việc phát sinh các số ngẫu nhiên phân bố theo U(0, 1). Để kiểm tra hàm này, bạn có thể tính trị trung bình của các số vừa phát sinh chẳng hạn, rồi so sánh kết quả với giá trị kì vọng  ∼ 0. 57721 / λ.

Lời giải: fischer_tippett.c

(3) phương sai mẫu Thiết kế, lập trình và kiểm tra một hàm để tính phương sai s2 của mẫu các điểm số liệu. Dùng trực tiếp PT (51), nghĩa là không dùng dạng tương đương của PT (21), vì dạng này dễ mắc phải sai số làm tròn. Mẫu hàm này như sau:

/********************** variance() ********************/ /** Tinh ra phuong sai cua n diem du lieu. **/ /** PARAMETERS: (*)= tham so tra lai **/ /** n: so diem du lieu. **/ /** x: mang du lieu **/ /** RETURNS: **/ /** phuong sai **/ /******************************************************/ double variance(int n, double *x)

Nhận xét: Thuộc toán có tên duyệt hai lần được hiệu chỉnh (corrected double-pass algorithm) [Chan và nnk. (1983)] nhằm mục đích giảm sai số xuống nữa. Nó được dựa theo phương trình s^2 = \frac{1}{n} \left[ \sum_{i=0}^{n-1} (x-\overline{x})^2 - \frac{1}{n} \left( \sum_{i=0}^{n-1} (x_i-\overline{x}) \right)^2 \right]\,. Thành phần thứ hai sẽ bằng 0 với những tính toán chính xác, và nếu có sai số thì sẽ dồn về số hạng thứ hai này. Công thức sẽ quan trọng hơn trong trường hợp giá trị kì vọng lớn. Hãy thực hiện những phép thử để phát sinh các số tuân theo phân bố Gauss với σ2 = 1 và μ = 1014, cả trường hợp có lẫn không có hiệu chỉnh.

Lời giải: variance.c

Tài liệu được trích dẫn

Chan, T. F., Golub, G.H., and LeVeque, R. J. (1983). Algorithm for Computing the Sample Variance: Analysis and Recommendations, Amer. Statist. 37, pp. 242-247

Devroye, L. (1986). Non-Uniform Random Variate Generation, (Springer, London).

Fernandez, J. F. and Criado, C. (1999). Algorithm for normal random numbers, Phys. Rev. E 60, pp. 3361-3365.

Ferrenberg, A. M., Landau, D. P. and Wong, Y. J. (1992). Monte Carlo Simulations: Hidden Errors from “Good” Random Number Generators, Phys. Rev. Lett. 69, pp. 3382-3384.

HotBits (trang web): ở đây bạn có thể đặt mua những file chứa số ngẫu nhiên được phát sinh từ sự phân rã chất phóng xạ, xem http://www.fourmilab.ch/hotbits/.

Morgan, B. J. T. (1984). Elements of Simulation, (Taylor & Francis).

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (1995). Numerical Recipes in C, (Cambridge University Press, Cambridge).

Quantis: Thiết bị phần cứng phát sinh số ngẫu nhiên, dựa trên một quá trình cơ lượng tử, đó là tán xạ photon. Có thể nối thiết bị này với máy tính qua cổng USB hoặc khe cắm PCI. Bạn có thể tìm thêm thông tin ở http://www.idquantique.com/products/quantis.htm.

Vattulainen, I., Ala-Nissila, T. and Kankaala, K. (1994). Physica Test for Random Numbers in Simulations, Phys. Rev. Lett. 73, pp. 2513-2516.

Westphal Electronic (http://www.westphal-electronic.com/) có bán những thiết bị tạo số ngẫu nhiên thực sự, dựa trên nhiễu nhiệt trong các đi-ốt Z. Các thiết bị này có thể được nối vào máy tính qua cổng USB hoặc bluetooth.

Đánh giá:

Chia sẻ:

  • Facebook
  • Email
  • Twitter
  • PrintFriendly
  • RSS Feed
Thích Đang tải...

Có liên quan

Từ khóa » Số Giả Ngẫu Nhiên Là Gì