Statistical Inference: Tiệm cận (Asymptotics)

Coin plot 10000

Trong bài viết này ta sẽ thảo luận về tiệm cận (asymptotics), làm thế nào để miêu tả dáng điệu của thống kê khi kích thước mẫu ngày càng tiến đến vô cùng. Giả định kích thước mẫu và kích thước quần thể là vô cùng hữu ích cho việc suy diễn thống kê và xấp xĩ.

Luật số lớn

Luật số lớn (Law of Large Numbers (LLN)) phát biểu rằng giá trị trung bình (mean) có xu hướng tiệm cận về giá trị mà nó đang ước lượng. Ta quan sát những mô phỏng của chúng ta rằng kích thước mẫu dữ liệu càng lớn thì ước lượng quần thể càng chính xác. Như chúng ta tung một con xúc sắc nhiều lần liên tiếp, nó có xu hướng hội tụ về giá trị xác suất 0.5. LLN tạo cơ sở cho kiểu suy luận dựa trên tần xuất. Lấy ví dụ minh hoạ với hàm coinPlot(). Hàm này nhận vào đối số n là số lượt tung xúc sắc. Sau mỗi lần tung, hàm này sẽ cộng dồn các giá trị 1 (head) và 0 (tail) sau đó tính trung bình cộng lại. Sau đó, hàm này sẽ xuất ra biểu đồ biến thiên của giá trị trung bình vừa mới tính được này.

coinPlot <- function(n){
  means <- cumsum(sample(0 : 1, n , replace = TRUE)) / (1  : n)
  g <- ggplot(data.frame(x = 1 : n, y = means), aes(x = x, y = y))
  g <- g + geom_hline(size=1.5 ,yintercept = 0.5,alpha=0.6,
                      linetype="longdash") + geom_line(size = 1)
  if(n<100){
    g <- g + geom_point(colour="red",size=3,alpha=0.8)
  }
  g <- g + labs(x = "Number of obs", y = "Cumulative mean")
  g <- g + scale_x_continuous(breaks=seq(0,n+1,ceiling(n/10)))
  print(g)
  invisible()
}

Ví dụ tung xúc sắc 10 lần

Coin plot 10

Coin plot 10

Biểu đồ xuất ra phụ thuộc vào giá trị ngẫu nhiên mà R tạo được. Vì vậy, biểu đồ trên có vẻ khác một chút khi chạy trên máy tính của bạn. Nếu bạn thử chạy hàm trên coinPlot(10) vài lần bạn sẽ thấy mỗi lần như vậy sẽ cho ra kết quả khác một chút. Bây giờ, ta thử tung xúc sắc 10,000 lần

Coin plot 10000

Coin plot 10000

Bạn đã thấy được sự khác biệt phải không nào. Đây là minh chứng cho lý thuyết tiệm cận (asymptotics). Đường biến thiên hội tụ tại giá trị 0.5. Ta bảo rằng một ước lượng thống kê như vậy thì đồng nhất nếu nó hội tụ tại giá trị mà nó đang ước lượng. Định lý số lớn bảo rằng trung bình mẫu của biến ngẫu nhiên iid thì đồng nhất với trung bình của quần thể. Tương tự như vậy, phương sai mẫu (sample variance) cũng có xu hướng hội tụ về phương sai của quần thể (population variance).

Định lý giới hạn trung tâm

Định lý giới hạn trung tâm (Central Limit Theorem (CLT)) là một trong những định lý quan trọng trong thống kê. Định lý phát biểu rằng phân phối của các giá trị trung bình của các biến ngẫu nhiên iid (đã được chuẩn hoá – normalized) hội tụ về các giá trị chuẩn khi kích thước mẫu dữ liệu tăng lên. Ta hãy cùng nhau “mổ xẻ” phát biểu trên. Đầu tiên “được chuẩn hoá” nghĩa là ta cần biến đổi trung bình mẫu X’ bằng cách trừ đi trung bình quần thể mu, sau đó chia cho sigma/sqrt(n). Trong đó, sigma là độ lệch chuẩn (standard deviation) của quần thể và n là kích thước mẫu dữ liệu. Thứ hai, CLT phát biểu rằng khi n lớn thì biễn ngẫu nhiên đã được chuẩn hoá (normalized variable) (X’-mu)/(sigma/sqrt(n)) này có xu hướng trở thành phân phối chuẩn với mean bằng 0 và variance bằng 1. Ta có thể phát biểu lại CLT như sau. Giả sử X_1, X_2, … X_n độc lập với nhau và là các biến ngẫu nhiên được phân phối như nhau (iid) từ quần thể có kích thước vô hạn với mean mu và variance sigma^2. Khi đó, nếu n lớn, trung bình (mean) của X đặt là X’, sẽ xấp xĩ với mean mu và variance sigma^2/n. Ta kí hiệu X’~N(mu,sigma^2/n). Để minh hoạ CLT trực quan hơn ta quan sát hình sau.

nosim <- 1000
cfunc <- function(x, n) sqrt(n) * (mean(x) - 3.5) / 1.71
dat <- data.frame(
  x = c(apply(matrix(sample(1 : 6, nosim * 10, replace = TRUE),
                     nosim), 1, cfunc, 10),
        apply(matrix(sample(1 : 6, nosim * 20, replace = TRUE),
                     nosim), 1, cfunc, 20),
        apply(matrix(sample(1 : 6, nosim * 30, replace = TRUE),
                     nosim), 1, cfunc, 30)
  ),
  size = factor(rep(c(10, 20, 30), rep(nosim, 3))))
g <- ggplot(dat, aes(x = x, fill = size)) + geom_histogram(alpha = .20, binwidth=.3, colour = "black", aes(y = ..density..))
g <- g + stat_function(fun = dnorm, size = 2)
g <- g + facet_grid(. ~ size)
print(g)
CLT dice

CLT dice

Hình trên biểu diễn 3 histogram của 1000 giá trị trung bình khi thực hiện 1000 lần lấy mẫu trên kích thước mẫu n (10, 20, 30). Mỗi giá trị trung bình của n mẫu (n=10,20,30) được chuẩn hoá bằng cách trừ cho mean (3.5) sau đó chia cho độ lỗi chuẩn (standard error), sqrt(2.92/n). Việc chuẩn hoá này khiến cho các histogram có dạng như phân phối chuẩn (với mean bằng 0 và độ lệch chuẩn bằng 1). Xem thêm video: https://youtu.be/JNm3M9cqWyc Chú ý rằng CLT không nói gì về phân bố ban đầu của quần thể là phân phối chuẩn. Ta có thể giả định tính phân phối chuẩn của trung bình quần thể bất chấp phân bố của quần thể ta đang có. Chỉ cần kích thước mẫu dữ liệu đủ lớn và các mẫu dữ liệu độc lập với nhau. Ta hãy quan sát định lý này hoạt động như thế nào với thí nghiệm tung đồng xu. Nhắc lại rằng nếu xác suất tung được mặt phải (head) là p thì xác suất tung được mặt trái (tail) là 1-p. Giá trị trung bình E(X) là p và phương sai Var(X) là p-p^2 hoặc p(1-p). Giả sử ta tung đồng xu n lần và gọi p’ là trung bình cộng của n lần tung đồng xu đó. Ta chuẩn hoá p’ bằng cách trừ cho mean p và chia cho độ lệch chuẩn sqrt(p(1-p)/n). Ta thực hiện thí nghiệm này 1000 lần và xuất ra biểu đồ histogram.

nosim <- 1000
cfunc <- function(x, n) 2 * sqrt(n) * (mean(x) - 0.5)
dat <- data.frame(
  x = c(apply(matrix(sample(0:1, nosim * 10, replace = TRUE),
                     nosim), 1, cfunc, 10),
        apply(matrix(sample(0:1, nosim * 20, replace = TRUE),
                     nosim), 1, cfunc, 20),
        apply(matrix(sample(0:1, nosim * 30, replace = TRUE),
                     nosim), 1, cfunc, 30)
  ),
  size = factor(rep(c(10, 20, 30), rep(nosim, 3))))
g <- ggplot(dat, aes(x = x, fill = size)) + geom_histogram(binwidth=.3, colour = quot;blackquot;, aes(y = ..density..))
g <- g + stat_function(fun = dnorm, size = 2)
g <- g + facet_grid(. ~ size)
print(g)
CLT fair coin

CLT fair coin

Ta thấy khi n bằng 30, histogram có xu hướng hội tụ tại mean 0. Tương tự với ví dụ trên, lần này ta sẽ mô phỏng đồng xu không công bằng với xác suất xảy ra mặt phải (head) là 0.9. Lúc này E(x) = 0.9 và độ lỗi chuẩn (standard error) bằng sqrt(0.09/n).

nosim <- 1000
cfunc <- function(x, n) sqrt(n) * (mean(x) - 0.9) / sqrt(.1 * .9)
dat <- data.frame(
  x = c(apply(matrix(sample(0:1, prob = c(.1,.9), nosim * 10, replace = TRUE),
                     nosim), 1, cfunc, 10),
        apply(matrix(sample(0:1, prob = c(.1,.9), nosim * 20, replace = TRUE),
                     nosim), 1, cfunc, 20),
        apply(matrix(sample(0:1, prob = c(.1,.9), nosim * 30, replace = TRUE),
                     nosim), 1, cfunc, 30)
  ),
  size = factor(rep(c(10, 20, 30), rep(nosim, 3))))
g <- ggplot(dat, aes(x = x, fill = size)) + geom_histogram(binwidth=.3, colour = "black";, aes(y = ..density..))
g <- g + stat_function(fun = dnorm, size = 2)
g <- g + facet_grid(. ~ size)
print(g)
CLT unfair coin

CLT unfair coin

Lần nữa, ta thấy khi n lấy mẫu càng lớn phân bố của các giá trị trung bình càng xấp xĩ với phân phối chuẩn, mặc dù ta đã chủ quan hoá đồng xu.

Khoảng tin cậy

Bây giờ chúng ta sẽ nói về khoảng tin cậy (confidence intervals). Từ CLT với n đủ lớn, trung bình mẫu là mean mu và độ lệch chuẩn sigma/sqrt(n). Ta biết rằng 95% diện tích nằm phía dưới đường cong nằm trong khoảng độ lệch chuẩn này.

Standard normal

Standard normal

Biểu đồ trên là phân phối chuẩn với mu=0 và sigma=1, để minh họa điều trên, toàn bộ bóng đỏ mô tả diện tích nằm trong độ lệch chuẩn 2, bóng đậm hơn mô tả 68% diện tích nằm trong độ lệch chuẩn 1. Theo đó, 5% diện tích còn lại không được tô bóng đỏ. Do đường cong đối xứng nên chỉ có 2.5% diện tích lớn hơn mean + độ lệch chuẩn 2 (mu+2sigma/sqrt(n)) và 2.5% diện tích nhỏ hơn mean – độ lệch chuẩn 2 (mu-2sigma/sqrt(n)). Vậy xác suất trung bình mẫu X’ lớn hơn mu+2sigma/sqrt(n) hay nhỏ hơn mu-2sigma/sqrt(n) là 5%. Tương tự, xác suất nằm trong giới hạn này là 95%. Tất nhiên ta có thể có nhiều khoảng tin cậy khác nhau. Nếu ta muốn giá trị khác 95, ta có thể sử dụng phân vị (quantile) thay vì 2. Phân vị X’ +/- 2 sigma/sqrt(n) được gọi là 95% khoảng tính từ mu. 95% nghĩa là nếu ta liên tục lấy mẫu với kích thước n thì khoảng 95% giá trị này nằm trong khoảng có chứa mu (số lượng chúng ta cần ước lượng). Lưu ý để tính khoảng tin cậy 95% ta chia (100-95) cho 2 (vì đường tiệm cận hai bên biểu đồ đối xứng) và công kết quả cho 95, ta sẽ có phân vị cần tính. Phân vị 97.5 thật ra là 1.96 nhưng để đơn giản ta làm tròn thành 2. Nếu ta muốn tìm khoảng tin cậy 90% thì phân vị chúng ta được tính như sau (100-90)/2 + 90 = 95. Ta dùng hàm qnorm của R để tìm phân vị 95 cho phân phối chuẩn. Hàm này lấy đối số là xác suất. Các đối số còn lại ta có thể lấy mặc định

qnorm(.95)
[1] 1.644854

Như ta đã thấy trước đó, trong phân phối nhị phân (binomial distribution), p đại diện cho xác suất biến cố xảy ra và phương sai sigma^2 được tính bằng p(1-p). Do đó, độ lỗi chuẩn của trung bình mẫu p’ là sqrt(p(1-p)/n), trong đó n là kích thước mẫu. Khoảng tin cậy 95% của p là p’ +/- 2sqrt(p(1-p)/n). Ở đây 2 thể hiện cho phân vị 97.5%. Điểm mấu chốt ở đây là ta không biết được giá trị thật của p, giá trị mà ta cần ước lượng. Làm thế nào ta tính được khoảng tin cậy khi không biết p(1-p) là bao nhiêu? Ta sẽ thận trọng và cố gắng cực đại hóa giá trị này để có được khoảng tin cậy cao nhất. Trong tích phân p(1-p) đạt cực đại khi p=1/2 nên ta có khoảng tin cậy lớn nhất 95% khi ta lấy giá trị p=1/2 trong công thức p’+/- 2sqrt(p(1-p)/n). Nếu p=1/2 ta có công thức cho khoảng tin cậy 95% là p’+/- 1/sqrt(n). Một kĩ thuật khác để tính khoảng tin cậy cho phân phối nhị phân đó là thế p bằng p’. Đây gọi là khoảng tin cậy Wald. Ta có thể dùng hàm qnorm để có giá trị phân vị chính xác (xấp xĩ 1.96) thay vì làm tròn 2. Với công thức p’+/- qnorm(.975)*sqrt(p'(1-p’)/100), với p’=0.6 ta có thể tính được khoảng tin cậy.

.6 + c(-1,1)*qnorm(.975)*sqrt(.6*.4/100)
[1] 0.5039818 0.6960182

Ta có thể tin khoảng tin cậy Wald bằng cách sử dụng hàm binom.test với tham số là 60 và 100 (100 người bỏ phiếu, trong đó có 60 phiếu ủng hộ) và các tham số khác mặc định. Hàm này thực hiện kiểm định giải thuyết rỗng (null hypothesis) cho kiểm định Bernoulli (đảm bảo độ phủ, sử dụng rất nhiều tính toán và không dựa vào CLT). Hàm này trả về khá nhiều thông tin, nếu chúng ta cần lấy thông tin khoảng tin cậy thì chỉ cần lấy x$conf.int.

binom.test(60,100)

	Exact binomial test

data:  60 and 100
number of successes = 60, number of trials = 100, p-value = 0.05689
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.4972092 0.6967052
sample estimates:
probability of success
                   0.6 

# chỉ lấy giá trị khoảng tin cậy
binom.test(60,100)$conf.int
[1] 0.4972092 0.6967052
attr(,quot;conf.levelquot;)
[1] 0.95

Ta hãy chuyển qua phân phối Poisson và khoảng tin cậy của nó. Nhắc lại phân phối Poisson áp dụng cho tổng số (counts) và tần xuất (rates). Về sau, ta viết X~Poisson(lambdat) với lambda là tổng số trung bình trên một đơn vị thời gian và t là tổng số thời gian theo dõi. Giả sử một máy bơm hạt nhân thất bại 5 lần trong 94.32 ngày và chúng ta muốn biết khoảng tin cậy tỉ lệ thất bại mỗi ngày là bao nhiêu. Số thất bại X là phân phối Poisson với thông số (lambdat). Chúng ta không quan sát tần xuất thất bại mà ước lượng x/t. Gọi ước lượng của chúng ta là lambda_hat, lambda_hat=x/t. Theo lý thuyết, phương sai của ước lượng tần xuất thất bại là lambda/t. Một lần nữa, ta không quan sát lambda nên ta sẽ sử dụng ước lượng. Ta xấp xĩ phương sai của lambda_hat bằng lambda_hat/t. Trong ví dụ trên, ta ước lượng trung bình x/t là 5/94.32 tỉ lệ thất bại trên một đơn vị thời gian.

lamb <- 5/94.32

Như vậy biến lamb là ước lượng trung bình và lamb/t là ước lượng cho phương sai. Công thức chúng ta sử dụng để tính khoảng tin cậy 95% là mean + c(-1,1)qnorm(.975)sqrt(est var).

lamb +c(-1,1)*qnorm(.975)*sqrt(lamb/94.32)
[1] 0.006545667 0.099476386

Ta có thể kiểm tra kết quả trên bằng hàm poisson.test trong R với đối số là 5 và 94.32. Đây là một phép kiểm chính xác nên nó đảm bảo độ phủ (coverage). Như phép kiểm trên phân phối nhị phân, ta chỉ cần quan tâm giá trị khoảng tin cậy bằng cách sử dụng x$conf.

poisson.test(5,94.32)$conf
[1] 0.01721254 0.12371005
attr(,quot;conf.levelquot;)
[1] 0.95

Kết quả gần đúng phải không nào. Bây giờ, để kiểm tra độ phủ ước lượng của chúng ta, ta sẽ chạy mô phỏng phân phối nhị phân. Ta sẽ biến thiên giá trị lambda của chúng ta từ 0.005 đến 0.1 với giá trị tăng từng bước (steps) là 0.01 (do đó, ta có 10 giá trị trong khoảng này). Với từng giá trị thu được, ta sẽ phát sinh 1000 mẫu cho phân phối Poisson với mean lambda*t. Ta sẽ tính các trung bình mẫu và sử dụng chúng để tính khoảng tin cậy 95%. Sau đó, ta đếm tần xuất của 1000 trung bình mẫu mô phỏng này (lambda) chứa trong khoảng tin cậy vừa tính toán trên.

Poisson Demo

Poisson Demo

Biểu đồ trên biểu diễn cho kết quả tính toán được. Ta thấy rằng độ phủ (coverage) được cải thiện khi lambda càng lớn. Bây giờ, ta hãy xem độ phủ được cải thiện nhiều hơn khi ta tăng đơn vị thời gian t=1000. Dưới đây là biểu đồ như ví dụ trên với t=1000. Ta thấy rằng độ phủ tốt hơn nhiều với hầu hết các giá trị của lambda trừ những giá trị quá nhỏ.

Poisson improvement

Poisson improvement

Tổng kết

Ta hãy điểm qua một số vấn đề đã nêu ở trên. Định lý số lớn phát biểu cho chúng ta biết rằng các giá trị trung bình của mẫu dữ liệu iid hội tụ về các giá trị trung bình của quần thể. Định lý giới hạn trung tâm phát biểu cho chúng ta biết rằng các giá trị trung bình xấp xĩ phân phối chuẩn khi kích thước lấy mẫu càng lớn. Và phân phối này có trọng tâm đặt tại trung bình của quần thể, độ lệch chuẩn xấp xỉ với độ lỗi chuẩn. Để tính khoảng tin cậy (confidence interval) cho mean ta lấy trung bình mẫu +/- cho tích của phân vị chuẩn (normal quantile) và độ lỗi chuẩn (standard error). Với khoảng tin cậy là 95% ta cần +/- cho độ lỗi chuẩn 2 lần hay chính xác là 1.96 lần. Khi dữ liệu của bạn có phương sai nhỏ thì khoảng tin cậy của chúng ta cũng nhỏ. Nếu kích thước mẫu dữ liệu càng lớn thì khoảng tin cậy của chúng ta càng nhỏ (khi kích thước mẫu càng lớn, phương sai mẫu càng nhỏ). Như vậy chúng ta đã hoàn tất bài học về tiệm cận (asymptotics). Hy vọng chúng ta đã tự tin hơn (confidence) và tiệm cận hơn (asymtomatic) với kiến thức về thống kê của mình.

Nguồn tham khảo:

Advertisements

Trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com Đăng xuất / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Đăng xuất / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Đăng xuất / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Đăng xuất / Thay đổi )

Connecting to %s