Số ngẫu nhiên trong R: rnorm, runif, sample và set.seed

set.seed() đi trước tiên

Các giá trị ngẫu nhiên thay đổi sau mỗi lần chạy thì vô dụng cho việc dạy học, chấm bài, gỡ lỗi hay nghiên cứu khoa học: nếu mô phỏng của bạn hôm nay cho ra 0,146 và ngày mai cho ra 0,153, con số nào sẽ vào báo cáo? set.seed() cố định điểm bắt đầu của bộ sinh số, khiến toàn bộ dòng "ngẫu nhiên" tiếp theo lặp lại chính xác:

set.seed(42)
print(rnorm(3))

set.seed(42)          # same seed...
print(rnorm(3))       # ...same "random" numbers

set.seed(7)           # different seed, different stream
print(rnorm(3))

Cùng seed, cùng giá trị rút ra, mọi lần, trên mọi máy. Những con số này là giả ngẫu nhiên: một dãy tất định được thiết kế để vượt qua mọi phép kiểm định thống kê về tính ngẫu nhiên, với seed quyết định bạn bắt đầu từ đâu trong dãy. Bản thân giá trị seed không mang ý nghĩa gì, 42, 7 hay 20260807 đều được; chỉ cần ghi lại. Thói quen cần xây dựng: một lệnh set.seed() ở đầu mọi script có dùng tính ngẫu nhiên. (Lưu ý rằng mỗi lần rút giá trị đều làm trạng thái tiến lên, nên thứ tự các lệnh gọi cũng ảnh hưởng đến khả năng tái tạo.)

Hệ thống d/p/q/r: một bảng giải mã cả thư viện

R đặt tên mọi hàm phân phối theo dạng tiền tố + họ phân phối, và khi đã thấy bảng này, bạn đọc được toàn bộ thư viện thống kê:

p và q là hàm ngược của nhau, và đây chính là cặp bạn đã gặp trong khoảng tin cậy dưới dạng qt(0.975, df):

print(pnorm(1.96))      # P(Z <= 1.96), about 0.975
print(qnorm(0.975))     # which z has 97.5% below it? About 1.96
print(dnorm(0))         # height of the bell curve at its peak

Thay tên họ phân phối và mọi thứ vẫn giữ nguyên: runif/dunif/punif/qunif, rbinom/..., rpois/..., rt/..., rexp/.... Học bốn tiền tố, dùng được hàng chục phân phối.

rnorm(): rút giá trị từ phân phối chuẩn

rnorm(n, mean, sd) rút giá trị từ một đường cong hình chuông, công cụ chủ lực để mô phỏng dữ liệu kiểu số đo:

set.seed(42)
iq <- rnorm(1000, mean = 100, sd = 15)

print(mean(iq))
print(sd(iq))
print(mean(iq > 130))   # share more than 2 sd above the mean

Trung bình và độ lệch chuẩn của mẫu rơi vào gần 100 và 15, chứ không đúng bằng: khoảng chênh đó chính là nhiễu lấy mẫu, và nó thu hẹp khi n tăng. Dòng cuối là một mẹo đáng học: mean() của một vector logic là tỷ lệ các giá trị TRUE, và tỷ lệ vượt quá 130 ra gần với giá trị lý thuyết 1 - pnorm(130, 100, 15) ≈ 2,3%. Giá trị mặc định là mean = 0, sd = 1 (phân phối chuẩn tắc). Một histogram của iq cho thấy hình chuông quen thuộc.

runif(), rbinom(), rpois()

Ba họ phân phối nữa đáp ứng hầu hết nhu cầu mô phỏng:

set.seed(42)

print(runif(5))                    # uniform on [0, 1]
print(runif(3, min = 10, max = 20))

print(rbinom(8, size = 10, prob = 0.5))   # successes in 10 coin flips, 8 times

print(rpois(8, lambda = 3))        # event counts averaging 3 per interval

runif() rải các giá trị đều trên một khoảng ("uniform", không phải "run if", ai cũng từng đọc nhầm một lần). rbinom(n, size, prob) mô phỏng n thí nghiệm, mỗi thí nghiệm gồm size lần thử, và trả về số lần thành công của mỗi thí nghiệm; vậy mỗi giá trị ở trên là số lần ra mặt ngửa trong 10 lần tung đồng xu. rpois(n, lambda) tạo ra số lần xảy ra của các sự kiện độc lập với tốc độ trung bình đã biết: số yêu cầu hỗ trợ mỗi giờ, số lỗi chính tả mỗi trang.

sample(): lấy mẫu và xáo trộn

Trong khi các hàm r* tạo ra giá trị mới từ một phân phối, sample() rút từ những giá trị bạn đã có sẵn:

set.seed(42)

print(sample(1:10, 3))                              # 3 values, no repeats
print(sample(c("win", "lose"), 5, replace = TRUE))  # 5 draws, repeats allowed
print(sample(1:5))                                  # no size = shuffle everything

Đối số replace là mấu chốt: FALSE (mặc định) giống như chia bài, mỗi giá trị chỉ xuất hiện một lần, và yêu cầu nhiều hơn số bạn có sẽ báo lỗi; TRUE giống như gieo xúc xắc, mỗi lần rút đều bắt đầu lại. Lấy mẫu có hoàn lại từ chính dữ liệu của bạn là cỗ máy của bootstrap. Khi chỉ truyền một vector, sample(x) trả về một hoán vị ngẫu nhiên: đây là cách viết thường dùng để xáo trộn.

Lấy mẫu các hàng của một data frame thì dùng sample() bên trong phép lập chỉ mục hàng:

set.seed(42)

random_cars <- mtcars[sample(nrow(mtcars), 5), ]
print(random_cars[, c("mpg", "wt", "hp")])

sample(nrow(mtcars), 5) chọn 5 số thứ tự hàng ngẫu nhiên; phép lập chỉ mục lấy ra những hàng đó. Đây là thao tác chuẩn để kiểm tra nhanh một bộ dữ liệu lớn hoặc chia tập train/test.

Một mô phỏng Monte Carlo nhỏ

Đây là thành quả của cả bộ công cụ. Câu hỏi: một quy trình cho ra các số đo có phân phối N(100, 15); bạn lấy trung bình của 10 số đo, xác suất để trung bình đó vượt quá 105 là bao nhiêu? Thay vì suy ra đáp án, hãy mô phỏng nó: làm thí nghiệm 10.000 lần rồi đếm:

set.seed(42)

trial <- function() {
  mean(rnorm(10, mean = 100, sd = 15)) > 105
}

results <- replicate(10000, trial())
print(mean(results))

# The exact answer, for comparison:
print(1 - pnorm(105, mean = 100, sd = 15 / sqrt(10)))

Kết quả mô phỏng chỉ lệch vài phần nghìn so với giá trị chính xác (khoảng 0,146). Đó là Monte Carlo trong một hơi thở: viết một lần thử thành một hàm, replicate() nó hàng nghìn lần, lấy mean() của các lần thành công. Ở đây có đáp án chính xác vì bài toán đơn giản như trong sách giáo khoa; ngay khi câu hỏi trở nên phức tạp (phân phối kỳ lạ, giá trị lớn nhất của các lần rút có tương quan, một trò chơi theo luật), con đường giải tích khép lại còn công thức mô phỏng vẫn chạy nguyên vẹn. Để ý 15 / sqrt(10) trong phép kiểm tra: giá trị trung bình dao động ít hơn từng giá trị riêng lẻ, cũng chính là quy luật căn bậc hai quyết định độ rộng của khoảng tin cậy.

Đủ tốt cho mô phỏng, không đủ cho bảo mật

Một giới hạn cần tôn trọng: bộ sinh số mặc định của R (Mersenne Twister) được xây dựng cho chất lượng thống kê và tốc độ, không phải cho bí mật. Đầu ra của nó là tất định khi biết seed, và trạng thái bên trong có thể được khôi phục từ đầu ra quan sát được: những lỗ hổng chết người với mật khẩu, token hay bất cứ thứ gì liên quan đến bảo mật. Cho mô phỏng, bootstrap và dạy học thì nó rất tuyệt; cho mật mã học, hãy dùng thư viện chuyên dụng (ví dụ package openssl), không bao giờ dùng sample() hay runif().

Những điều cần nhớ

• Một lệnh set.seed() ở đầu script giúp mọi kết quả "ngẫu nhiên" tái tạo được: cùng seed, cùng giá trị.
• Các tiền tố d/p/q/r giải mã cả thư viện phân phối: rút ngẫu nhiên, mật độ, xác suất tích lũy, phân vị.
• rnorm(n, mean, sd), runif(n, min, max), rbinom(n, size, prob), rpois(n, lambda) đáp ứng hầu hết nhu cầu mô phỏng.
• sample() rút từ chính các giá trị của bạn: replace = TRUE cho xúc xắc, mặc định cho chia bài, không có size để xáo trộn; df[sample(nrow(df), k), ] lấy mẫu các hàng.
• Monte Carlo = hàm một lần thử + replicate() + mean(): công thức trả lời những câu hỏi xác suất mà toán học khó tiếp cận thuận tiện.
• Bộ sinh số ngẫu nhiên của R dành cho mô phỏng, không dành cho mật mã học.

Tiếp theo: gỡ lỗi, ý nghĩa thật sự của các thông báo lỗi trong R và cách đọc traceback.