Chương 2 Lý thuyết xác suất

Section này sẽ ôn lại một vài khái niệm cơ bản của Lý thuyết xác suất và cách áp dụng chúng trong R.

Hầu hết các chức năng thống kê trong phiên bản R tiêu chuẩn đều có trong gói lệnh stats. Gói lệnh này cung cấp ngoài những hàm toán cơ bản dùng để tính toán các phép đo mô tả dữ liệu còn cung cấp một lượng lớn các hàm phân phối xác suất. Gói lệnh này được cài mặc định khi dùng R và do vậy, ta không cần thiết phải chạy lệnh: install.packages('stats').

2.1 Biến ngẫu nhiên

Một số khái niệm cơ bản trong lý thuyết xác suất:

  • Những kết quả loại trừ lẫn nhau của một quá trình ngẫu nhiên gọi là outcomes. Loại trừ lẫn nhau có nghĩa là chỉ một trong các kết quả có thể có được quan sát.

  • Xác suất của một outcome là phần trăm outcome này xảy ra trong dài hạn, nghĩa là khi thí nghiệm được lặp lại rất thường xuyên.

  • Bộ các outcome có thể có của một biến ngẫu nhiên được gọi là không gian mẫu.

  • Một event là bộ nhỏ của không gian mẫu và chứa một outcome hoặc nhiều hơn.

Tất cả những ý này đều được gói gọn trong khái niệm biến ngẫu nhiên, là tổng hợp số học của các outcome ngẫu nhiên. Biến ngẫu nhiên có thể rời rạc hoặc liên tục.

2.2 Phân phối xác suất rời rạc

Một ví dụ điển hình cho biến ngẫu nhiên rời rạc \(D\) là kết quả của việc lăn xúc sắc. Ở đây không gian mẫu là \(\{1,2,3,4,5,6\}\) và ta có thể nghĩ về rất nhiều event, chẳng hạn, như những outcome nằm giữa \(2\)\(5\).

Một hàm cơ bản để lấy mẫu ngẫu nhiên từ bộ thành tố là hàm sample(). Ta có thể sử dụng nó để giả lập outcome ngẫu nhiên từ một lần lăn xúc sắc.

sample(1:6,1)
## [1] 1

Xác suất biến ngẫu nhiên rời rạc là bộ các giá trị có thể có của biến và tổng xác suất phải bằng 1. Hàm phân phối tích luỹ cho kết quả xác suất mà biến ngẫu nhiên nhỏ hơn hoặc bằng một giá trị nào đó.

Đối với lăn 1 xúc sắc, phân phối xác suất và phân phối tích luỹ được tổng kết trong bảng sau.

X Freq Cumul Relative CumulativeRelative
1 1 1 0.1666667 0.1666667
2 1 2 0.1666667 0.3333333
3 1 3 0.1666667 0.5000000
4 1 4 0.1666667 0.6666667
5 1 5 0.1666667 0.8333333
6 1 6 0.1666667 1.0000000

2.2.1 Phân phối Bernoulli

Trong trường hợp bộ thành tố mà hàm sample() chọn ra không chứa các con số. Ta có thể giả lập việc tung đồng xu với hai outcome là \(H\) (đầu) và \(T\) (đuôi).

sample(c("H","T"),1)
## [1] "H"

Kết quả của một lần tung đồng xu là biến ngẫu nhiên tuân theo phân phối Bernoulli, nghĩa là một biến nhị phân.

Giả sử tung 10 lần đồng xu này liên tục, và câu hỏi đặt ra là: Xác suất được \(5 H\) là bao nhiêu? Đây là ví dụ cho một thí nghiệm Bernoulli trong đó \(n=10\) lần thử Bernoulli và ta đang quan sát \(k=5\) lần tung thành công ra \(H\). Với xác suất cho một lần \(H\)\(p=0.5\), ta biết rằng số lần thành công \(k\) tuân theo phân phối nhị phân:

\[ k \sim B(n,p) \]

Xác suất quan sát thấy \(k\) lần thành công với thí nghiệm \(B(n,p)\) sẽ là

\[ f(k) = P(k) = {n \choose k} p^k (1-p)^{n-k} = \frac{n!}{k!(n-k)!}p^k(1-p)^{n-k} \]

trong đó \({n \choose k}\) là hệ số nhị phân.

Trong R, ta có thể giải quyết bài toán trên bằng công cụ là hàm dbinom():

dbinom(x=5, size = 10, prob = 0.5)
## [1] 0.2460938

Ta kết luận rằng xác suất quan sát 5 lần thành công trong 10 lần tung xu là khoảng 24.6%.

Bài tập: Viết code cho \(P(4\le k \le 7)\):

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjZGJpbm9tIiwic29sdXRpb24iOiJkYmlub20oeD00LCBzaXplID0gMTAsIHByb2IgPSAwLjUpICsgZGJpbm9tKHg9NSwgc2l6ZSA9IDEwLCBwcm9iID0gMC41KSArIGRiaW5vbSh4PTYsIHNpemUgPSAxMCwgcHJvYiA9IDAuNSkgKyBkYmlub20oeD03LCBzaXplID0gMTAsIHByb2IgPSAwLjUpIn0=

Đáp án đúng cho trường hợp này là \(0.7734375\) (~ \(77.34\%\)). Có nhiều cách để ra được kết quả này, phổ biến nhất là tính riêng từng trường hợp đối với $ k = 4,5,6,7$ sau đó cộng tổng lại. Cách này dễ nhất nhưng việc code rất rườm rà. Để rút gọn ta có thể dùng cấu trúc dbinom() và hàm sum() như sau

sum(dbinom(x=4:7, size = 10, prob = 0.5))
## [1] 0.7734375

Một cách đơn giản khác là sử dụng định nghĩa Hàm phân phối luỹ kế \(F(x) = P(X \le x)\) để tính toán. Khi đó: \[ P(a \le X \le b) = F(b) - F(a)\] với \(F\) là hàm phân phối luỹ kế.

pbinom(size = 10, prob = 0.5, q = 7) - pbinom(size = 10, prob = 0.5, q = 3)
## [1] 0.7734375

Để hình hoạ hoá phân phối của thí nghiệm này, ta dùng commands sau:

k <- 0:10
probability <- dbinom(x = k, size = 10, prob = 0.5)
plot(x = k , y = probability, main = "Probability Distribution Function", type = "h")

Viết code hình hoạ hoá phân phối luỹ kế của thí nghiệm trên.

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjVmlzdWFsaXphdGlvbiIsInNvbHV0aW9uIjoicGJpbm9tIn0=

2.3 Giá trị kỳ vọng, Trung bình

Với biến ngẫu nhiên \(Y\)\(k\) giá trị có thể có \(y_1,...,y_k\) trong đó \(y_1\) là giá trị đầu tiên, \(y_2\) là giá trị thứ hai,… và xác suất \(Y\) lấy giá trị \(y_1\)\(p_1\), xác suất \(Y\) lấy giá trị \(y_2\)\(p_2\),… thì giá trị kỳ vọng của biến \(Y\), \(E(Y)\), được định nghĩa như sau:

\[E(Y) = \sum_{i=1}^k y_ip_i \]

Giá trị kỳ vọng của biến \(Y\) còn được gọi là giá trị trung bình của nó, và được ký hiệu là \(\mu_Y\).

Trong ví dụ xúc sắc, biến ngẫu nhiên, \(D\) nhận \(6\) giá trị có thể có \(d_1 = 1, d_2 = 2, ..., d_6 = 6\). Với một xúc sắc cân đối, mỗi trường hợp đều nhận xác suất là \(1/6\). Do đó dễ dàng tính tay được giá trị kỳ vọng của biến \(D\): \[E(D) = 1/6 \sum_1^6 d_1 = 3.5\]

Với R ta có thể dễ dàng tính giá trị này bằng hàm mean()

mean(1:6)
## [1] 3.5

Khi tiến hành lấy mẫu có thay thế bằng cách tung xúc sắc 3 lần liên tục, ta có thể thực hiện trên R như sau

set.seed(123)
sample(1:6, 3, replace = TRUE)
## [1] 2 5 3

Trong thực tế ta thường tiến hành số lần thí nghiệm nhiều hơn 3 lần, có thể lên đến \(10000\) phép thử:

set.seed(123)
mean(sample(1:6,10000, replace = TRUE))
## [1] 3.4718

Ta thấy rằng giá trị trung bình rất gần giá trị kỳ vọng. Điều này sẽ được bàn kỹ hơn trong các section tiếp theo.

2.4 Phương sai

Một phép đo khác thường xuyên sử dụng đến khi khảo sát xác suất biến ngẫu nhiên chính là phương sai và độ lệch chuẩn. Cả hai đều đo lường độ phân tán của biến ngẫu nhiên.

Phương sai của biến ngẫu nhiên \(Y\), ký hiệu là \(\sigma^2_Y\) được định nghĩa: \[ \sigma^2_Y = \text{Var}(Y) = E\left[(Y-\mu_Y)^2 \right] = \sum_{i=1}^k (y_i - \mu_Y)^2 p_i \]

Độ lệch chuẩn của \(Y\)\(\sigma_Y\), là căn bậc hai của phương sai. Đơn vị tính của độ lệch chuẩn giống như đơn vị của biến.

Đối với R việc tính toán phương sai theo công thức định nghĩa là dành cho tổng thể thì R không thực hiện được. Thay vào đó, R tính toán phương sai (cũng như độ lệch chuẩn) theo mẫu: \[ s_Y^2 = \frac{1}{n-1} \sum_{i=1}^n (y_i - \bar{y})^2 \]

Nhớ rằng phương sai mẫu \(s^2_Y\) khác với phương sai tổng thể \(\sigma^2_Y\): \[ \sigma^2_Y = \frac{1}{N} \sum_{i=1}^N (y_i - \mu_Y)^2\]

Điều này sẽ rõ ràng hơn khi thí nghiệm với xúc sắc và biến \(D\). Phương sai tổng thể: \[ \text{Var}(D) = 1/6 \sum_{i=1}^6 (d_i - 3.5)^2 = 2.92 \]

trong khi đó R tính toán theo hàm var() như sau

var(1:6)
## [1] 3.5

Phương pháp tính phương sai mẫu là một estimator của phương sai tổng thể.

2.5 Phân phối xác suất liên tục

Biến ngẫu nhiên liên tục có vùng giá trị liên tục, chẳng hạn như \([0,1]\). Do vậy ta không thể dùng lối tư duy biến ngẫu nhiên rời rạc lắp vào biến liên tục được. Thay vào đó, phân phối xác suất liên tục được tổng hợp trong hàm mật độ xác suất probability density function (PDF). Với \(f_Y(y)\) là PDF của biến liên tục \(Y\), xác suất \(Y\) nằm giữa \(a\)\(b\) với \(a<b\) là: \[ P(a \le Y \le b) = \int_a^b f_Y(y) dy \]

Hàm phân phối xác suất luỹ kế ** cumulative probability distribution function (CDF) ** được định nghĩa tương tự như trường hợp rời rạc.

Giá trị kỳ vọng của biến \(Y\) được tính như sau \[ E(Y) = \mu_Y = \int yf_Y(y) \]

Phương sai được tính theo công thức sau \[ \text{Var}(Y) = \sigma^2_Y = \int (y-\mu_Y)^2f_Y(y) \]

Chẳng hạn, với biến liên tục \(X\) có hàm mật độ phân phối như sau: \[ f_X(x) = \frac{3}{x^4}, x>1 \]

Ta dễ dàng chứng minh được với vùng giá trị đã cho thì tích phân \(f_X\) bằng 1, ứng với xác suất \(P(X>1) = 100\%\): \[ \int_1 f_X(x)dx = \int_1 \frac{3}{x^4} dx = \lim_{t \rightarrow \infty} \int_1^t \frac{3}{x^4} dx = \lim_{t \rightarrow \infty} -x^{-3} \Big|_{x=1}^t = -\left(\lim_{t \rightarrow \infty} \frac{1}{t^3} - 1 \right) = 1 \] Kỳ vọng \(X\) được tính như sau:

\[\int_1 xf_X(x)dx = \int_1 x.\frac{3}{x^4} dx = \lim_{t \rightarrow \infty} \int_1^t \frac{3}{x^3} dx = \lim_{t \rightarrow \infty} -\frac{3}{2}x^{-2} \Big|_{x=1}^t = -\frac{3}{2}\left(\lim_{t \rightarrow \infty} \frac{1}{t^2} - 1 \right) = \frac{3}{2}\]

Chú ý rằng \(\text{Var}(X) = E(X^2)-E(X)^2\), nên để tính \(\text{Var}(X)\) ta chỉ cần tính \(E(X^2)\):

\[\int_1 x^2f_X(x)dx = \int_1 x^2.\frac{3}{x^4} dx = \lim_{t \rightarrow \infty} \int_1^t \frac{3}{x^2} dx = \lim_{t \rightarrow \infty} -3x^{-1} \Big|_{x=1}^t = -3\left(\lim_{t \rightarrow \infty} \frac{1}{t} - 1 \right) = 3\]

Do đó: \[\text{Var}(X) = 3 - \frac{9}{4} = \frac{3}{4}\]

Tuy nhiên, đối với những PDF không tồn tại dạng tích phân đúng, phương pháp này không thể áp dụng được. May mắn là, R cho phép ta dễ dàng tính toán kết quả trên với công cụ dùng là hàm integrate(). Phương pháp làm như sau

f <- function(x) 3/x^4
g <- function(x) x*f(x)
h <- function(x) x^2*f(x)
area <-  integrate(f, lower = 1, upper = Inf)$value
EX <- integrate(g, lower = 1, upper = Inf)$value
VarX <- integrate(h, lower = 1, upper = Inf)$value
cat('Vùng xác suất = ', area)
## Vùng xác suất =  1
cat('Giá trị kỳ vọng =', EX)
## Giá trị kỳ vọng = 1.5
cat(' Phương sai = ', VarX)
##  Phương sai =  3

Có rất nhiều phân phối xác suất, nhưng phổ biến nhất đối với biến liên tục đó là các phân phối Normal, Chi bình phương, Student, và F. Do đó, ta sẽ bàn đến các hàm trong R đối với những phân phối này liên quan đến các bài toán phổ biến: mật độ, xác suất và quantile.

Mọi phân phối xác suất mà R xử lý có 4 hàm chính theo format: chữ cái đầu + tên. Chẳng hạn, với phân phối Normal, tên trong Rnorm và 4 chữ cái đầu ứng với công dụng như sau:

  • d cho density - hàm PDF
  • p cho probability - hàm CDF
  • q cho quantile - hàm quantile (tức CDF nghịch đảo)
  • r cho random - phương pháp tạo lập số liệu ngẫu nhiên (random number generator)

Như vậy R có 4 hàm chính cho phân phối Normal: dnorm(), pnorm(), qnorm(), rnorm().

2.5.1 Phân phối Normal

Nhờ vào phân phối Standard Normal và Định lý giới hạn tập trung Central Limit Theorem, phân phối Normal trở thành phân phối quan trọng bậc nhất. Phân phối này được quy định bởi hai tham số, trung bình \(\mu\) và độ lệch chuẩn \(\sigma\), ký hiệu \(N(\mu,\sigma)\). PDF của phân phối Normal như sau

\[ f(x) = \frac{1}{\sqrt{2\pi} \sigma } \exp \left( \frac{-(x-\mu)^2}{2\sigma^2 } \right) \]

Với phân phối Standard Normal thì $ = 0 $, và $ = 1$. Biến Standard Normal được ký hiệu là \(Z\). Thông thường, PDF của Standard Normal ký hiệu là \(\phi\) và CDF ký hiệu là \(\Phi\). Do đó:

\[ \phi(x) = \Phi '(x) \] trong đó: \(\Phi(x) = P(Z \le x )\), và $ Z N(0,1) $

Trong R, ta có thể vẽ PDF của phân phối Normal như sau:

curve(dnorm(x), xlim = c(-3.5, 3.5), ylab = "Density", xlab = "Standard Normal")

Ta có thể tính toán PDF của Normal cho nhiều giá trị cùng một lúc:

dnorm(x = c(-1.96, 0, 1.96))
## [1] 0.05844094 0.39894228 0.05844094

Bài tập: Vẽ CDF của phân phối Normal

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjQ0RGIiwic29sdXRpb24iOiJwYmlub20ifQ==

Ta cũng có thể dùng R để tính toán xác suất events liên quan đến biến Standard Normal.

Giả sử ta cần tính $ P(Z 1.337) $. \[ P(Z \le 1.337) = \int_{-\infty}^{1.337} \phi (x)dx \]

Không hề có một đáp án đúng cho tích phân trên. Tuy nhiên R cung cấp một ước lượng xấp xỉ bằng hàm integrate() như ta đã biết.

f <- function(x) {
  1/(sqrt(2 * pi)) * exp(-0.5 * x^2)
}

quants <- c(-1.96, 0, 1.96)
f(quants)
## [1] 0.05844094 0.39894228 0.05844094
f(quants) == dnorm(quants)
## [1] TRUE TRUE TRUE
integrate(f, lower = -Inf, upper = 1.337)
## 0.9093887 with absolute error < 1.7e-07

R cho thấy kết quả hàm f() mà ta đã code bằng với hàm dnorm() có sẵn. Và xác suất quan sát biến \(Z \le 1.337\) vào khoảng \(90.93 \%\). Một cách khác ngoài integrate() đó là dùng hàm pnorm().

Một kết quả thường dùng trong Kinh tế lượng chính là \(95\%\) xác suất biến ngẫu nhiên Standard Normal nằm trong vùng giá trị \([-1.96,1.96]\), tức khoảng cách vào khoảng 2 lần độ lệch chuẩn so với giá trị trung bình.

Phân phối Standard Normal hoặc Normal nói chung có hình chuông đối xứng. Sự đối xứng được thể hiện qua ví dụ sau:

\[P(-1.96 \le Z \le 1.96) = 1 - 2 \times P(Z \le -1.96)\]

1 - 2*(pnorm(-1.96))
## [1] 0.9500042

Một phân phối mở rộng của Normal trong nền tảng đơn biến chính là phân phối Normal đa biến multivariate Normal distribution. Hàm mật độ hợp joint PDF có dạng như sau:

\[\begin{aligned}g_{X,Y}(x,y) =\frac{1}{2\pi \sigma_X \sigma_Y \sqrt{1-\rho^2}} \exp\left[ \frac{1}{-2(1-\rho^2)}\left( \left(\frac{x-\mu_X}{\sigma_X}\right)^2 \\ -2\rho\left(\frac{x-\mu_X}{\sigma_X}\right)\left(\frac{y-\mu_Y}{\sigma_Y}\right)+\left(\frac{y-\mu_Y}{\sigma_Y}\right)^2 \right) \right] \end{aligned}\]

Phuơng trình trên chính là PDF nhị biến Normal \(X\)\(Y\) với độ tương quan \(\rho\). Nếu \(X\)\(Y\) đều tuân theo Standard Normal và độ tương quan giữa chúng bằng 0 thì joint PDF trở thành \[ g_{X,Y}(x,y) = \frac{1}{2\pi}\exp \left[ -\frac{1}{2} (x^2+y^2) \right] \]

Hình vẽ dưới đây cho thấy hình hoạ hoá hàm PDF trên. Trỏ chuột vào hình để thấy sự kết hợp giữa hai biến Standard Normal \(X\) \(Y\).

2.5.2 Phân phối Chi bình phương

Đây là phân phối hay được sử dụng trong các bài toán kiểm định mô hình hồi quy.

Tổng của \(M\) bình phương biến ngẫu nhiên Standard Normal sẽ tuân theo phân phối này với bậc tự do là \(M\): \[ Z_1^2 + ... + Z_M^2 = \sum_{m=1}^M Z_m^2 \sim \chi^2_M \]

trong đó \(Z_m \sim i.i.d. \phi(0,1)\).

Phân phối Chi bình phương bậc tự do \(M\) sẽ có kỳ vọng tại \(M\), mode tại \(M-2\) với \(M \ge 2\) và phương sai \(2\times M\).

Sử dụng code dưới đây để vẽ PDF và CDF của phân phối \(\chi_3^2\) để làm ví dụ.

curve(dchisq(x, df = 3),xlim=c(0,10), ylim = c(0,1), col = "blue", ylab = "", main = "PDF và CDF của Chi bình phương, M=3")
curve(pchisq(x, df = 3), xlim =c(0,10), add = TRUE, col = 'red')
legend('topleft', c('PDF','CDF'), col = c('blue', 'red'), lty = c(1,1))

Bởi vì kết quả phân phối Chi bình phương luôn dương, nên tập giá trị của PDF và CDF phải là \(\mathbb{R}_{\ge 0}\).

Chú ý rằng kỳ vọng và phương sai đều phụ thuộc vào số bậc tự do, do đó khi \(M\) thay đổi thì hình dạng phân phối cũng sẽ thay đổi.

curve(dchisq(x, df=1), xlim = c(0,10), xlab = 'x', ylab='Mật độ', main = 'Phân phối Chi bình phương')
for (M in 2:7){
  curve(dchisq(x,df = M), xlim = c(0,10), add = TRUE, col = M)
}
legend('topright', as.character(1:7), col = 1:7, lty = 1, title = 'Bậc tự do')

Việc tăng bậc tự do sẽ dịch chuyển phân phối sang phải (vì mode tăng dần), và tăng độ phân tán (vì độ lệch chuẩn tăng dần).

2.5.3 Phân phối Student t

Với \(Z\) tuân theo Standard Normal, \(W\) tuân theo \(\chi_M^2\)\(Z,W\) độc lập với nhau, thì biến \(X\) tuân theo phân phối Student t được định nghĩa: \[ X = \frac{Z}{\sqrt{W/M}}\]

Khi đó \(X\sim t_M\) với \(M\) là bậc tự do.

Giống với phân phối Chi bình phương, hình dạng \(t_M\) sẽ thay đổi phụ thuộc vào \(M\), tuy nhiên nhìn chung đều có dạng chuông đối xứng. Khi \(M\) càng lớn, thì hình dạng \(t_M\) càng giống với phân phối Normal. Tuy nhiên, nó không đồng nghĩa với việc: với \(M\) đủ lớn thì phân phối \(t_M\) được xấp xỉ bởi phân phối Standard Normal. Ngược lại, khi \(M \ge 30\) thì phân phối \(t_M\) xấp xỉ phân phối Standard Normal. Cụ thể, \(t_\infty\) chính là phân phối Standard Normal.

Với điều kiện \(M>1\) thì phân phối này có kỳ vọng và với \(M>2\) thì tồn tại phương sai: \[E(X) = 0 , M>1\] \[\text{Var}(X) = \frac{M}{M-2}, M>2\]

Bây giờ ta dùng R để vẽ vài hình hoạ PDF của các phân phối \(t\) và so sánh chúng với Standard Normal.

curve(dnorm(x), xlim = c(-4,4), xlab = 'x', lty = 2, ylab = 'Mật độ', main = 'Mật độ phân phối t')
for (M in c(2,4,25)){
  curve(dt(x, df = M), xlim = c(-4,4), add = TRUE, col = M)
}
legend('topright', c('N(0,1)','M=2', 'M=4', 'M=25'), col = c(1,2,4,25), lty = c(2,1,1,1))

Đồ thị cho thấy khi \(M\) càng tăng, hình dạng \(t_M\) càng gần giống Standard Normal. Khi \(M\) nhỏ, phân phối \(t_M\) có đuôi nặng hơn so với Standard Normal, ta gọi là dạng chuông mập.

2.5.4 Phân phối F

Một tỷ số giữa các biến ngẫu nhiên cũng quan trọng trong kinh tế lượng đó là tỷ số giữa hai biến độc lập \(\chi^2\): \[ \frac{W/M}{V/n} \sim F_{M,n},\text{ trong đó } W \sim \chi^2_M, V \sim \chi^2_n \]

Giả sử \(Y\) tuân theo \(F_{3,14}\) và ta đang quan tâm đến \(P(Y \ge 2)\). R sẽ tính như sau:

pf(2,3,14, lower.tail = FALSE)
## [1] 0.1603538

Chú ý argument lower.tail đặt FALSE nếu ta tính đuôi phải, mặc định là TRUE để tính đuôi trái. Hình sau minh hoạ xác suất này.

2.6 Lấy mẫu ngẫu nhiên

Để làm rõ khái niệm lấy mẫu ngẫu nhiên, ta xem xét ví dụ tung xúc xắc.

Giả sử ta tung \(n\) lần. Ở đây ta đang quan tâm giá trị của \(Y_1, Y_2,..., Y_n\). Vì xúc xắc được tung ngẫu nhiên, do vậy bản thân biến \(Y\) là một biến ngẫu nhiên, và tập giá trị của nó được thu thập khi ta xây dựng một mẫu, nghĩa là ta tung \(n\) lần. Mỗi quan sát được chọn ngẫu nhiên từ tập tổng thể từ \(1\) đến \(6\), kéo theo phân phối mỗi giá trị \(Y\) là giống nhau.

Ta biết rằng giá trị bất kỳ biến \(Y_i\) không ảnh hưởng đến giá trị các biến \(Y_j\) khác của mẫu, nghĩa là việc thay đổi giá trị biến \(Y_i\) không làm thay đổi bản chất xác suất của biến \(Y\). Do đó mọi giá trị \(Y_i\) tuân theo một phân phối độc lập.

Kết quả từ hai luận cứ trên, cho thấy \(Y \sim i.i.d\). Như vậy trong lấy mẫu ngẫu nhiên đơn thuần, \(n\) đối tượng được chọn ngẫu nhiên từ một tổng thể. Mỗi đối tượng đều có khả năng xuất hiện như nhau. Khi đó, \(Y_1,...,Y_n\) tuân theo phân phối giống nhau và độc lập.

Một câu hỏi được đặt ra là: điều gì xảy ra nếu ta gán một hàm cho một mẫu ngẫu nhiên thay vì một tổng thể? Một mẫu bao gồm 2 lần lấy mẫu từ tập \(\{1,2,3,4,5,6\}\). Bất kỳ một hàm số nào đối với hai biến này cũng sẽ là ngẫu nhiên. Ví dụ, ta thực hiện đoạn code sau nhiều lần.

sum(sample(1:6,2,replace = TRUE))
## [1] 8

Kết quả của phép tổng, ta gọi là \(S\), là một biến ngẫu nhiên vì nó phụ thuộc vào hai biến ngẫu nhiên khác, giả sử \(D_1\)\(D_2\). Tập giá trị có thể có của \(S\) từ \(6^2 = 36\) cặp kết hợp có thể có của hai biến \(D_i\) là: \(\{2,3,4,5,6,7,8,9,10,11,12\}\). Bảng phân phối được tổng hợp như sau.

S 2.0000000 3.0000000 4.0000000 5.0000000 6.0000000 7.0000000 8.0000000 9.0000000 10.0000000 11.0000000 12.0000000
PS 0.0277778 0.0555556 0.0833333 0.1111111 0.1388889 0.1666667 0.1388889 0.1111111 0.0833333 0.0555556 0.0277778

Ta quan tâm đến kỳ vọng và phương sai của \(S\).

S <- 2:12
PS <- c(1:6, 5:1)/ 36
ES <- S %*% PS
VarS <- (S-c(ES))^2 %*% PS
d <- cbind(ES,VarS)
colnames(d) <- c('ES','VarS')
d
##      ES     VarS
## [1,]  7 5.833333

Toán tử %*% dùng để nhân một số với các phần tử của một vector.

Ta có thể thấy rằng phân phối \(S\) khác hẳn phân phối của \(D\). Đoạn code R dưới đây dùng để minh hoạ sự khác biệt này.

par(mfrow = c(1,2))
probability <- rep(1/6, 6)
names(probability) <- 1:6
barplot(PS, 
        ylim = c(0, 0.2), 
        xlab = "S", 
        ylab = " Xác suất", 
        col = "steelblue", 
        space = 0, 
        main = "Tổng hai lần tung")
barplot(probability, 
        ylim = c(0, 0.2), 
        xlab = "D", 
        col = "steelblue", 
        space = 0, 
        main = "Một lần tung")

2.7 Trung bình và phương sai của trung bình mẫu

Giả sử \(Y_1, ..., Y_n\)\(i.i.d\) với \(\mu_Y\)\(\sigma_Y^2\) là trung bình và phương sai của \(Y_i\). Trung bình mẫu của biến \(Y\) được định nghĩa như sau: \[ \bar{Y} = \frac{1}{n} \sum_{i=1}^nY_i \]

\(\bar{Y}\) là một biến ngẫu nhiên, do đó kỳ vọng và phương sai của nó chính là:

\[E(\bar{Y}) = E \left( \frac{1}{n} \sum^n Y_i \right) = \frac{1}{n} E \left(\sum^n Y_i \right) = \frac{1}{n} \sum^n E(Y_i) = \frac{1}{n}\times n \times \mu_Y = \mu_Y\]

$$ \[\begin{aligned} \text{Var}(\bar{Y}) &= \text{Var} \left( \frac{1}{n} \sum^n Y_i \right) \\ & = \frac{1}{n^2}\sum^n \text{Var}(Y_i) + \frac{1}{n^2} \sum^n_{i=1} \sum^n_{j=1, j \ne i} \text{cov}(Y_i,Y_j) \\ & = \frac{\sigma^2_Y}{n} \end{aligned}\]

$$

Kết quả trên đúng với bất kể phân phối của \(Y\) là gì. Ta có thể xác nhận với trường hợp ví dụ là phân phối Normal thông qua giả lập Monte Carlo trong R. Tiến trình giả lập như sau:

  • Chọn quy mô mẫu n và số lượng mẫu cần lấy reps.
  • Lấy n quan sát từ một mẫu, lặp lại reps lần.
  • Tính giá trị trung bình của từng mẫu.
n <- 10
reps <- 10000
set.seed(123)
samples <- replicate(reps, rnorm(n))
sample.avgs <- colMeans(samples)
hist(sample.avgs, ylim = c(0,1.4), col = "steelblue", freq = F, breaks = 20)
curve(dnorm(x, sd = 1/sqrt(n)), col = "red", lwd = 2, add = T)
legend('topright','Sample mean', col = 'red', lwd = 2)

Dễ thấy rằng phân phối thực tiễn rất gần với phân phối lý thuyết \(N(0,0.1)\) nếu các biến \(Y\) tuân theo \(\phi(0,1)\)\(n=10\).

Ta có thể lấy ví dụ khác là phân phối Chi bình phương \(\chi^2_3\). Ta biết rằng, phân phối Chi bình phương bậc \(M\) được phát sinh từ phân phối của tổng của \(M\) biến ngẫu nhiên Standard Normal độc lập nhau. Do đó thuật toán Monte Carlo trong trường hợp này là:

  • Chọn bậc tự do, DF, và số lượng mẫu reps.
  • Lặp lại reps lần việc lấy mẫu quy mô DF từ phân phối Standard Normal qua hàm replicate().
  • Với mỗi mẫu, bình phương kết quả và cộng tất cả lại bằng hàm colSums().
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjU3RlcCAxXG5cblxuI1N0ZXAgMlxuc2V0LnNlZWQoMTIzKVxuXG4jU3RlcCAzXG5cbiNIaXN0b2dyYW0iLCJzb2x1dGlvbiI6InJlcHMgPC0gMTAwMDBcbkRGIDwtIDNcbnNldC5zZWVkKDEyMylcblogPC0gcmVwbGljYXRlKHJlcHMsIHJub3JtKERGKSlcblggPC0gY29sU3VtcyhaXjIpXG5oaXN0KFgsIGNvbCA9IFwic3RlZWxibHVlXCIsIGZyZXEgPSBGLCBicmVha3MgPSA0MClcbmN1cnZlKGRjaGlzcSh4LCBkZiA9IERGKSwgdHlwZSA9ICdsJywgY29sID0gXCJyZWRcIiwgbHdkID0gMiwgYWRkID0gVClcbmxlZ2VuZCgndG9wcmlnaHQnLCdTYW1wbGUgbWVhbicsIGNvbCA9ICdyZWQnLCBsd2QgPSAyKSJ9

2.8 Định lý Giới hạn tập trung

Khi nói về các phân phối lấy mẫu sampling distributions thì sẽ có hai hướng tiếp cận: chính xácxấp xỉ.

Ví dụ về cách tiếp cận chính xác nằm ở ví dụ phần trên, với mục tiêu là tìm ra phân phối đúng nhất cho mọi quy mô mẫu \(n\). Với cách tiếp cận này, ta có thể viết ra giấy công thức tính toán phân phối của một hàm số bất kỳ. Tuy nhiên điều kiện ràng buộc hay gọi là giả định rằng, ta biết rõ phân phối thành phần. Trong thực tế ta thường không may mắn như vậy. Do vậy cách tiếp cận thứ hai được sử dụng nhiều hơn.

Ở cách tiếp cận thứ hai, mục tiêu là tìm ra phân phối xấp xỉ đúng trong trường hợp \(n\) đủ lớn. Ý tưởng ở đây chính là, khi \(n \rightarrow \infty\) thì phân phối xấp xỉ này chính là phân phối đúng. Hai cách tiếp cận cho cùng một kết quả. Khi làm thí nghiệm nghiên cứu, cần lấy mẫu sao cho sự khác biệt giữa hai cách tiếp cận này không nhiều.

Quy luật số lớn nói rằng, khi mẫu lớn, \(\bar{Y}\) sẽ tiệm cận \(\mu_Y\) (giá trị đúng) với xác suất rất cao. Định lý giới hạn tập trung nói rằng, phân phối mẫu của biến được chuẩn hoá, \(\frac{\bar{Y}-\mu_Y}{\sigma_{\bar{Y}}}\) tiệm cận phân phối Normal.

Để dễ hình dung quy luật số lớn, giả sử ta liên tục lặp đi lặp lại thí nghiệm tung đồng xu. Ta biết rằng \(Y_i\) tuân theo phân phối Bernoulli với xác suất \(p\) nếu là mặt head, ứng với điều kiện \(Y_i=1\). Khi đó giá trị kỳ vọng đúng của biến \(Y\) là: \(\mu_Y = 0.5\).

Gọi \(R_n\) là phần trăm head quan sát được trong \(n\) lần tung đầu tiên, nghĩa là \[R_n = \frac{1}{n}\sum_{i=1}^n Y_i\]

Theo quy luật số lớn: $ R_n _Y $ khi $ n $. Đoạn code R sau sẽ thể hiện điều này với giá trị \(n\) rất lớn, có thể xem như \(\infty\).

set.seed(1)
N <- 30000
Y <- sample(0:1, N, replace = T)
S <- cumsum(Y)
R <- S/(1:N)
plot(R, ylim = c(0.3,0.7), type = 'l', col = 'steelblue', lwd = 2, xlab = 'n', ylab = 'R_n', main = 'Sự hội tụ của R khi lặp lại thí nghiệm')
lines(c(0,N), c(0.5,0.5), col = 'darkred', lty = 2, lwd = 1)

Theo CLT, phân phối của trung bình mẫu \(\bar{Y}\) của phân phối Bernoulli sẽ xấp xỉ phân phối Normal với trung bình \(0.5\) và phương sai \(0.5\times 0.5 = 0.25\) với \(n\) rất lớn. Hệ quả, trung bình mẫu chuẩn hoá \(\frac{\bar{Y}- 0.5}{0.5/\sqrt{n}}\) sẽ xấp xỉ tuân theo \(\phi(0,1)\). Ta sẽ tiến hành giả lập để xác nhận trường hợp này như sau.

  • Chọn một lượng mẫu lớn, chẳng hạn \(10000\), có quy mô \(n\) từ phân phối Bernoulli sau đó tính trung bình từng mẫu.
  • Chuẩn hoá giá trị trung bình này
  • So sánh histogram của các giá trị trung bình này với phân phối \(\phi(0,1)\)
  • Lặp lại quá trình với \(n\) lớn hơn
reps <- 10000
sample.sizes <- c(2,10,50,1000)
par(mfrow= c(2,2))
for (n in sample.sizes){
  sample.means <- rep(0, reps)
  std.sample.means <- rep(0, reps)
  set.seed(49)
  for (i in 1:reps){
    x <- rbinom(n,1,0.5)
    sample.means[i] <- mean(x)
    std.sample.means[i] <- sqrt(n)*(mean(x) - 0.5) / 0.5
  }
  hist(std.sample.means, col = 'steelblue', freq = FALSE, xlim = c(-3,3), ylim = c(0,0.4), xlab = paste('n =',n), main = '')
  curve(dnorm(x), lwd = 2, col = 'darkred', add = T)
}

Ta có thể thấy phân phối lấy mẫu giả định của trung bình chuẩn hoá tiến rất gần với \(\phi(0,1)\) khi \(n\) tăng từ \(50\) đến \(1000\). Tuy nhiên với quy mô mẫu nhỏ thì lại khác hoàn toàn so với \(\phi(0,1)\)