Chương 3 Thống kê cơ bản

Trong section này ta sẽ tìm hiểu các nội dung cơ bản xung quanh Thống kê:

  • Ước lượng tham số chưa biết tổng thể
  • Kiểm định
  • Khoảng tin cậy

Những phương pháp này được sử dụng cực kỳ nhiều trong Kinh tế lượng. Các vấn đề kinh tế lượng đa phần đều xảy ra trong trường hợp chưa biết tổng thể là gì.

3.1 Ước lượng trung bình tổng thể

Estimators là những hàm số suy luận dữ liệu tổng thể từ dữ liệu mẫu. Estimates là những giá trị số học được tính bởi estimator dựa trên dữ liệu mẫu. Estimator là ngẫu nhiên vì chúng là hàm số của dữ liệu ngẫu nhiên. Trong khi đó estimate thì không ngẫu nhiên.

Giả sử đứng dưới góc độ kinh tế, gọi \(Y\) là thu nhập theo giờ của sinh viên. Giả sử ta đang quan tâm \(\mu_Y\) trung bình của \(Y\). Để xác định chính xác \(\mu_Y\) ta cần phải phỏng vấn mọi sinh viên trong nền kinh tế. Điều này không bao giờ đạt được vì nguồn lực về tài chính và thời gian. Tuy nhiên, ta có thể lựa mẫu ngẫu nhiên của \( n\) quan sát \(i.i.d\) và ước lượng \(\mu_Y\) bằng cách sử dụng công thức sau:

\[ \bar{Y} = \frac{1}{n} \sum_{i=1}^n Y_i \]

Liệu đây có là một estimator tốt của \(\mu_Y\). Ta so sánh với một estimator đơn giản hơn, quan sát đầu tiên \(Y_1\). Giả sử \(Y_1 \sim \chi_{12}^2\) vì thu nhập thì không âm và phân phối thường nghiêng phải, những đặc điểm quá phù hợp với phân phối Chi bình phương. Cùng xem qua PDF của phân phối này.

Ta có thể tạo một mẫu quy mô \(100\) quan sát và lấy quan sát đầu tiên cho \(Y_1\) như sau.

set.seed(1)
y_1 = rchisq(n=100, df = 12)[1]
y_1
## [1] 8.257893

Ước lượng \(8.26\) không quá xa so với \(\mu_Y=12\). Vậy giữa \( {Y}\) và \(Y_1\) estimator nào tốt hơn trong việc phản ánh \(\mu_Y\)?

Một estimator tốt cần phải đảm bảo những yêu cầu sau:

  • Không thiên lệch unbiasedness: trung bình của phân phối lấy mẫu của một estimator phải bằng trung bình đúng: \(E(\hat{\mu_Y}) = \mu_Y\).
  • Tính nhất quán consistency: sự sai lệch của estimator phải giảm dần khi lượng quan sát mẫu tăng lên. Nghĩa là, xác suất mà ước lượng \(\hat{\mu_Y}\) phải nằm trong một khoảng hẹp xung quanh \(\mu_Y\) tiến gần về \(100\%\) khi \(n\) tăng.
  • Hiệu quả efficiency: giả sử ta có hai estimator, \(\hat{\mu_Y}\)\(\tilde{\mu_Y}\). Với cùng quy mô mẫu \(n\)\(E(\hat{\mu_y}) = E(\tilde{\mu_Y}) = \mu_Y\), nếu \(\text{Var}(\hat{\mu_Y}) < \text{Var}(\tilde{Y})\) thì \(\hat{\mu_Y}\) tốt hơn.

Để xem xét trung bình mẫu như một estimator cho trung bình tổng thể, ta thực hiện thuật toán R như sau:

  • tạo một tổng thể gồm \(10000\) quan sát từ phân phối \(N(10,1)\).
  • lấy mẫu ngẫu nhiên từ tổng thể này và tính \({Y}\) cho từng mẫu.
  • lặp lại \(N=25000\) lần.
pop <- rnorm(10000,10,1)

est1 <- replicate(expr = mean(sample(x = pop, size = 5)), n=25000) 
est2 <- replicate(expr = mean(sample(x = pop, size = 25)), n=25000)
fo <- replicate(expr = sample(x=pop, size = 5)[1], n=25000)
plot(density(fo),
  col = 'green', lwd = 2, ylim = c(0,2), xlab = 'Ước lượng', main = 'Phân phối lấy mẫu', ylab = 'Mật độ')
lines(density(est1),
  col = 'steelblue', lwd = 2, bty = 'l')
lines(density(est2),
  col = 'red2', lwd=2)
abline(v=10, lty=2)
curve(dnorm(x, mean = 10),
  lwd = 2, lty = 2, add = T)
legend('topleft', legend = c('N(10,1)',
  expression(Y[1]),
  expression(bar(Y) ~ n==5),
  expression(bar(Y) ~ n==25)
  ), lty = c(2,1,1,1), col = c('black', 'green', 'steelblue','red2'), lwd = 2
  )

Đầu tiên, mọi phân phối lấy mẫu (đường nét liền) đều tập trung xung quanh \(\mu=10\). Đây là bằng chứng cho thấy tính unbiasedness của \(Y_1, \bar{Y}_5,\bar{Y}_{25}\).

Tiếp tục nhìn độ phân tán của phân phối.

  1. Phân phối lấy mẫu của \(Y_1\) (màu xanh lá cây) rất khớp với PDF \(N(10,1)\) lý thuyết (đường gạch). Điều này không có gì ngạc nhiên vì \(Y_1\) được chọn ngẫu nhiên \(25000\) lần từ phân phối \(N(10,1)\).
  2. Cả hai phân phối lấy mẫu \(\bar{Y}\) đều cho thấy độ phân tán thấp hơn rất nhiều so với phân phối của \(Y_1\). Điều này cho thấy phương sai của \(\bar{Y}\) thấp hơn so với của \(Y_1\). Như vậy estimator \(\bar{Y}\) tốt hơn.
  3. Khi \(n\) tăng thì độ phân tán giảm dần, cho thấy estimator tốt hơn.

3.1.1 Bình phương nhỏ nhất

Giả sử ta có quan sát \(Y_1,...,Y_n\) của biến \(Y \sim N(10,1)\) (đương nhiên là ta không biết thông tin về tổng thể). Ta muốn tìm một estimator \(m\) nào đó sao cho dự đoán quan sát tốt nhất. Nói nôm na, ta chọn \(m\) sao cho tổng chênh lệch bình phương giữa giá trị dự đoán và giá trị quan sát là nhỏ nhất: \(\sum_{i=1}^n(Y_i - m)^2\).

Ta hãy xem \(Y_i-m\) như là một lỗi sai khi dự đoán \(Y_i\) bằng estimator \(m\) và nhiệm vụ của ta là tìm \(m\) để lỗi sai này nhỏ nhất. Ta cũng có thể tối thiểu hoá tổng của giá trị tuyệt đối chênh lệch này nhưng tối thiểu hoá tổng của bình phương chênh lệch về mặt toán học sẽ tiện hơn. Đó là lý do tại sao estimator có tên là bình phương nhỏ nhất. Và \(m=\bar{Y}\) là một estimator làm được điều này. Ta chứng minh thông qua R như sau.

sqm <- function(m){
  sum((y-m)^2)
}

sqm <- Vectorize(sqm)
y <- rnorm(100,10,1)
mean(y)
## [1] 9.918857
curve(sqm(x),from = -50, to = 70, xlab = 'm', ylab = 'Tổng bình phương')
abline(v = mean(y), lty = 2, col = 'darkred')
text(x = mean(y), y = 0, labels = paste(round(mean(y),2)))

Vì hàm số bình phương nhỏ nhất là bậc 2 nên chỉ có duy nhất một giá trị nhỏ nhất. Hình vẽ cho thấy giá trị này nằm chính xác tại trung bình mẫu của dữ liệu mẫu.

Cho đến hiện tại, ta giả định ngầm rằng dữ liệu quan sát \(Y_i\) đều được lấy từ quy trình lấy mẫu thoả mãn điều kiện lấy mẫu ngẫu nhiên. Giả định này thường được dùng khi ước lượng trung bình tổng thể bằng estimator \(\bar{Y}\). Nếu không thoả mãn, ước lượng sẽ bị thiên lệch.

Quay lại tổng thể pop đã được định nghĩa ở trên, gồm \(10000\) quan sát và ta có thể tính trung bình tổng thể thông qua hàm mean(). Bây giờ ta sẽ tính estimator \(\bar{Y}\) cho mẫu quy mô \(10\) được lặp lại \(25000\) lần. Tuy nhiên lần này thay vì lấy mẫu ngẫu nhiên (nghĩa là xác suất mỗi giá trị trong tổng thể xuất hiện ngang nhau) thì ta sẽ phân bổ xác suất cao hơn đối với \(2500\) quan sát nhỏ nhất.

est3 <- replicate(n=25000, expr = mean(sample(x=sort(pop), size = 10, prob = c(rep(4,2500),rep(1,7500)))))
plot(density(est2), col = 'steelblue', lwd = 2, xlim = c(8,11), xlab = '', ylab = '', main = '')
lines(density(est3), col = 'red2', lwd = 2)
abline(v = mean(pop), lty = 2, col = 'darkred', lwd = 2)
text(x = mean(pop), y = 0, labels = paste(mean(pop)))
legend('topleft', legend = c(expression(bar(Y)[n==25000] ~ ', không i.i.d'),
  expression(bar(Y)[n==25000] ~ 'i.i.d'), 'trung bình tổng thể'), lty = c(1,1,2), col = c('red2', 'steelblue', 'darkred'), lwd = 2)

Hình vẽ cho thấy, lấy mẫu không ngẫu nhiên (hay giả định \(i.i.d\) không được thoả mãn) thì ta sẽ ước lượng thấp giá trị trung bình tổng thể.

3.2 Kiểm định

3.2.1 Giả thuyết và kiểm định thống kê

Trong kiểm định, ta muốn khai thác thông tin có trong một mẫu để làm bằng chứng ủng hộ hay bác bỏ giả thuyết. Những giả thuyết là những câu hỏi nhị phân đơn giản mà được trả lời yes hay no. Trong kiểm định ta thường giải quyết hai giả thuyết:

  • Giả thuyết null \(H_0\): là giả thuyết mà ta cần kiểm định.
  • Giả thuyết alternative \(H_a/H_A/H_1\): giả thuyết được cho là đúng nếu \(H_0\) bị bác bỏ.

Đối với bài toán liên quan đến kiểm định trung bình tổng thể, cặp giả thuyết được quy định như sau

\[ \begin{cases} H_0: E(Y) = \mu_{Y,0} \\ H_1: E(Y) \ne \mu_{Y,0} \end{cases} \]

Đây gọi là kiểm định hai đuôi.

3.2.2 p-value

Giả sử rằng \(H_0\) đúng, thì \(p\)-value là giá trị xác suất của việc ta giết nhầm, nghĩa là ta bác bỏ \(H_0\), một giả thuyết đúng. Cụ thể: \[ p-\text{value} = P_{H_0}\left[ |\bar{Y}-\mu_{Y,0}| > |\bar{Y}^{act} - \mu_{Y,0}|\right] \]

trong đó, \(\bar{Y}^{act}\) là trung bình mẫu thực tính. Như vậy, để tính được \(p\)-value, ta cần biết về phân phối lâý mẫu của \(\bar{Y}\) khi \(H_0\) đúng. Tuy nhiên, hầu hết các trường hợp phân phối này ta không biết. Nhưng theo CLT thì \(\bar{Y} \approx N(\mu_{Y,0},\sigma^2_{\bar{Y}}), \text{ }\sigma^2_{\bar{Y}} = \frac{\sigma^2_Y}{n}\). Do đó:

\[ \frac{\bar{Y} - \mu_{Y,0}}{ \sigma_Y / \sqrt{n}} \sim \phi{N}(0,1)\]

Như vậy đối với mẫu quy mô lớn, \(p\)-value có thể được tính mà không cần biết chính xác phân phối lấy mẫu của \(\bar{Y}\).

Giả sử \(\sigma_{\bar{Y}}\) đã được biết. Khi đó ta có thể viết lại công thức \(p\)-value như sau: \[ \begin{aligned} p-\text{value} &= P_{H_0} \left[ \Big| \frac{\bar{Y}- \mu}{\sigma_{\bar{Y}}} \Big| > \Big| \frac{\bar{Y}^{act}- \mu}{\sigma_{\bar{Y}}} \Big| \right] \\ &= 2\times \Phi\left[-\Big|\frac{\bar{Y}^{act}- \mu}{\sigma_{\bar{Y}}} \Big| \right] \end{aligned}\]

Khi đó \(p\)-value được nhìn theo hình sau:

curve(dnorm(x),
  xlim = c(-4,4),
  main = 'p-value',
  yaxs = 'i', ylab = '', xlab = 'z', lwd = 2, axes = F)
axis(1,
  at = c(-1.5,0,1.5),
  padj = 0.75,
  labels = c(expression(-frac(bar(Y)^'act' ~-~ bar(mu)[Y,0],sigma[bar(Y)])),
    0,
    expression(frac(bar(Y)^'act' ~-~ bar(mu)[Y,0],sigma[bar(Y)]))))

polygon(x = c(-6, seq(-6,-1.5,0.01), -1.5),
  y = c(0,dnorm(seq(-6,-1.5,0.01)),0), col = 'steelblue')
polygon(x = c(1.5, seq(1.5,6,0.01), 6),
  y = c(0,dnorm(seq(1.5,6,0.01)),0), col = 'steelblue')

3.2.3 Phương sai mẫu, độ lệch chuẩn mẫu và sai số chuẩn

Nếu \(\sigma_Y^2\) chưa biết, nó phải được ước lượng thông qua phương sai mẫu: \[s^2_Y = \frac{1}{n-1} \sum_{i=1}^n (Y_i -\bar{Y})^2\]

Tương tự, độ lệch chuẩn mẫu \(s_Y = \sqrt{s^2_Y}\)estimator của độ lệch chuẩn tổng thể. Trong R độ lệch chuẩn mẫu được tính bằng hàm sd().

Sử dụng R ta có thể chứng minh \(s_Y\) là một estimator tốt cho độ lệch chuẩn tổng thể.

n <- c(10000, 5000, 2000, 1000, 500)
sq_y <- replicate(n = 10000, expr = sd(rnorm(n[1],10,10)))
plot(density(sq_y),
  main = expression('Phân phối lấy mẫu' ~ s[Y]),
  xlab = expression(s[Y]),
  lwd = 2)
for (i in 2:length(n)){
  sq_y <- replicate(n=10000, expr = sd(rnorm(n[i],10,10)))
  lines(density(sq_y), col = i, lwd = 2)
}
abline(v = 10, col = 6, lwd = 2)
text(x = 10, y=0,labels = '10')
legend('topleft', legend = c(expression(n==10000),
  expression(n==5000),
  expression(n==2000),
  expression(n==1000),
  expression(n==500),
  'Phương sai tổng thể'), col = 1:6, lwd = 2)

Đồ thị cho thấy phân phối của \(s_Y\) sẽ hẹp hơn và xoay quanh giá trị đúng \(\sigma_Y = 10\) khi \(n\) tăng.

Hàm số ước lượng độ lệch chuẩn của một estimator được gọi là sai số chuẩn. Với một biến \(Y \sim i.i.d\) thì phân phối lấy mẫu của nó là \(\frac{\sigma_Y^2}{n}\). Sai số chuẩn của \(\bar{Y}\), ký hiệu là \(SE(\bar{Y})\) là một estimator của độ lệch chuẩn của \(\bar{Y}\): \[SE(\bar{Y} = \hat{\sigma_Y} = \frac{s_Y}{n})\]

Ta tiếp tục dùng R để khảo sát \(SE(\bar{Y})\) với \(Y\) tuân theo phân phối Bernoulli \(i.i.d\)\(p=0.1\). Ta biết rằng \(E(Y) = 0.1\)\(\text{Var}(Y) = p(1-p)=0.09\). \(E(Y)\) được ước lượng bởi \(\bar{Y}\) và nếu ta xem xét mẫu quy mô \(n=100\) thì nó có phương sai: \[\sigma_{\bar{Y}}^2 = p(1-p)/n = 0.0009)\]

và độ lệch chuẩn: \[\sigma_{\bar{Y}} = \sqrt{p(1-p)/n} = 0.03 \]

mean_est <- numeric(10000)
se_est <- numeric(10000)
for (i in 1:10000){
  s <- sample(0:1, size = 100, prob = c(0.9,0.1),
    replace = T)

  mean_est[i] <- mean(s)
  se_est[i] <- sqrt(mean(s)*(1-mean(s))/100)
}
mean(mean_est)
## [1] 0.100065
mean(se_est)
## [1] 0.02956867

Ta thấy \(SE(\bar{Y}) = 0.0295687\) rất gần với \(\sigma_{\bar{Y}} = 0.03\) do đó nó là một estimator có tính nhất quán.

Trong trường hợp phương sai tổng thể chưa biết, \(p\)-value được tính như sau: \[ p-\text{value} = 2\times \Phi\left[-\Big|\frac{\bar{Y}^{act}- \mu}{SE(\bar{Y})} \Big| \right]\]

Trong R ta thực hiện như sau để tính:

samplemean_act <- mean(sample(0:1,
  prob = c(0.9,0.1), replace = T, size = 100))

mean_h0 <- 0.1

SE_samplemean <- sqrt(samplemean_act*(1-samplemean_act)/100)

pvalue <- 2*pnorm(-abs(samplemean_act - mean_h0)/ SE_samplemean)
pvalue
## [1] 0.2396777

3.2.4 Thống kê \(t\)

Trong kiểm định, trung bình mẫu chuẩn hoá được gọi là thống kê \(t\): \[t = \frac{\bar{Y}-\mu_{Y,0}}{SE(\bar{Y})}\]

Hiện tại, ta đã tính xong thống kê \(t\) trong R trong ví dụ trên: (samplemean_act - mean_h0)/ SE_samplemean. Nếu \(H_0\) đúng, thì thống kê \(T\) sẽ tuân theo phân phối xấp xỉ \(N(0,1)\) khi \(n\) lớn.

t_stat <- numeric(10000)
n <- 300
for (i in 1:10000){
  s <- sample(0:1, size = n, prob = c(0.9,0.1), replace = T)
  t_stat[i] = (mean(s)-0.1)/sqrt(var(s)/n)
}
plot(density(t_stat), xlab = 'Thống kê t', main = '', ylab = '',lwd = 2, col = 'steelblue', xlim = c(-4,4))
curve(dnorm(x), add = T, lty = 2, lwd = 2)

Ta thấy rõ ràng việc chọn giá trị \(n=300\) khá hợp lý vì phân phối của t_stat rất xấp xỉ phân phối \(N(0,1)\).

3.2.5 Kiểm định thống kê với mức ý nghĩa cho trước

Trong kiểm định, có hai loại lỗi sai có thể xảy ra:

  • Lỗi sai loại I: ta bác bỏ \(H_0\) mặc dù nó đúng.
  • Lỗi sai loại II: ta không bác bỏ \(H_0\) mặc dù nó sai.

Mức ý nghĩa của kiểm định là xác suất rơi vào lỗi sai loại I mà chúng ta chấp nhận trước. Mức ý nghĩa được chọn trước khi kiểm định được thực hiện.

Tính đến hiện tại, ta ra quyết định thống kê dựa trên việc so sánh hai tham số \(p\)-value và mức ý nghĩa. Chú ý cả hai đều liên quan đến lỗi sai loại I. Nếu \(p-\text{value} <\) mức ý nghĩa thì ta bác bỏ \(H_0\).

Một cách khác để ra quyết định chính là so sánh thống kê kiểm định và giá trị tới hạn của thống kê kiểm định. Giá trị tới hạn được xác định dựa trên mức ý nghĩa đã đề ra. Giá trị tới hạn này sẽ xác định ranh giới giữa hai vùng: vùng chấp nhậnvùng bác bỏ. Vùng chấp nhận chưa mọi giá trị của thống kê kiểm định mà không bác bỏ \(H_0\) trong khi đó vùng bác bỏ chứa mọi giá trị bác bỏ \(H_0\).

Với thống kê \(t\), giá trị tới hạn của nó tại ý nghĩa \(5\%\)\(1.96\).

Xác suất mà kiểm định bác bỏ giả thuyết đúng được gọi là quy mô kiểm định, chính là mức ý nghĩa. Trong khi đó, xác suất mà kiểm định bác bỏ đúng giả thuyết sai được gọi là sức mạnh kiểm định.

3.2.6 Kiểm định một đuôi

Đôi lúc chúng ta quan tâm đến kiểm định liệu trung bình nhỏ hơn hay lớn hơn một giá trị nào đó. Giả sử ta đang khảo sát sự khác biệt giữa mức lương ra trường của hai nhóm sinh viên, một được đào tạo bài bản và một ít được đào tạo tốt. Cặp giả thuyết dưới đây là một ví dụ điển hình cho kiểm định đuôi phải: \[ H_0: \mu_Y = \mu_{Y,0} \text{ vs } H_1: \mu_Y > \mu_{Y,0}\]

trong đó \(\mu_Y\) là trung bình lương theo giờ của nhóm được đào tạo tốt và \(\mu_{Y,0}\) là lương theo giờ trung bình của nhóm còn lại. Ta sẽ bác bỏ \(H_0\) nếu giá trị thống kê lớn hơn giá trị tới hạn \(1.64\) (thay vì \(1.96\) của kiểm định hai đuôi), đây là quantile \(95\%\) của \(N(0,1)\). Vùng bác bỏ như sau:

3.3 Khoảng tin cậy

Ta sẽ không bao giờ đo lường chính xác giá trị trung bình tổng thể của \(Y\) bằng cách sử dụng mẫu ngẫu nhiên. Tuy nhiên, ta có thể tính khoảng tin cậy cho trung bình tổng thể. Đó là khoảng có chứa trung bình tổng thể với một xác suất cho trước. Khoảng tin cậy \(95\%\) của trung bình tổng thể \(\mu_Y\) là: \[ \mu_Y = [ \bar{Y} \pm 1.96 \times SE(\bar{Y})\]

Trong R để tính khoảng tin cậy ta có thể dùng hàm t.test()ls() được xây sẵn trong gói stats.

set.seed(1)
sampledata <- rnorm(100,10,10)
ls(t.test(sampledata))
## [1] "alternative" "conf.int"    "data.name"   "estimate"    "method"     
## [6] "null.value"  "p.value"     "parameter"   "statistic"
t.test(sampledata)$'conf.int'
## [1]  9.306651 12.871096
## attr(,"conf.level")
## [1] 0.95

Kết quả cho thấy khoảng tin cậy \(95\%\)\([9.31, 12.87]\), khoảng này chứa giá trị đúng của trung bình là \(10\). Sau đây ta sẽ xem xét những thông tin từ hàm t.test().

## 
##  One Sample t-test
## 
## data:  sampledata
## t = 12.346, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   9.306651 12.871096
## sample estimates:
## mean of x 
##  11.08887

Ta thấy t.test() không chỉ tính khoảng tin cậy mà còn tự động thực hiện kiểm định hai đuôi cho giả thuyết \(H_0: \mu_Y = 0\) tại mức ý nghĩa \(5\%\). Trong ví dụ này, ta kết luận trung bình tổng thể có ý nghĩa thống kê (nghĩa là khác \(0\) về mặt thống kê) với mức \(5\%\).

3.4 Kiểm định so sánh

Giả sử ta quan tâm đến trung bình của hai tổng thể khác nhau, ký hiệu là \(\mu_1\)\(mu_2\), câu hỏi đặt ra là liệu chúng có khác nhau không. Cặp giả thuyết ở đây chính là: \[H_0: \mu_1 - \mu_2 = d_0 \text{ vs } H_1: \mu_1 - \mu_2 \ne d_0\]

trong đó \(d_0\) là khoảng chênh lệch giả định. Giả thuyết \(H_0\) được kiểm định nhờ vào thống kê \(t\): \[t = \frac{(\bar{Y_1} - \bar{Y_2}-d_0)}{SE(\bar{Y_1} - \bar{Y_2})}\]

trong đó: \(SE(\bar{Y_1} - \bar{Y_2}) = \sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}\).

Đối với mẫu quy mô lớn thì thống kê \(t\) tuân theo Standard Normal.

Trong R kiểm định này (mặc định \(d_0=0\)) được thực hiện như sau.

set.seed(1)
sample_pop1 <- rnorm(100,10,10)
sample_pop2 <- rnorm(100,10,20)
t.test(sample_pop1, sample_pop2)
## 
##  Welch Two Sample t-test
## 
## data:  sample_pop1 and sample_pop2
## t = 0.872, df = 140.52, p-value = 0.3847
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.338012  6.028083
## sample estimates:
## mean of x mean of y 
## 11.088874  9.243838

Ta thấy rằng không có bằng chứng bác bỏ \(H_0\).

3.5 Hệ số tương quan

Đồ thị điểm rải sử dụng dữ liệu hai chiều, \(n\) quan sát \(X_i\)\(Y_i\). Để hình dung rõ hơn, xem đoạn code sau.

set.seed(123)
X <- runif(n=100, min=18,max = 70)
Y <- X + rnorm(n=100, 50, 15)
plot(X,Y, type = 'p', main = 'Đồ thị điểm rải', col = 'steelblue', pch = 19)

Hiệp phương sai và hệ số tương quan của hai biến sẽ liên quan đến phân phối xác suất hợp của các biến này. Estimator ước lượng hiệp phương sai chính là hiệp phương sai mẫu. \[s_{XY} = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X})(Y_i - \bar{Y})\]

Tương tự vậy ta có hệ số tương quan mẫu \[r_{XY} = \frac{s_{XY}}{s_X s_Y}\]

Trong R để tính toán các đại lượng này ta làm như sau.

cov(X,Y)
## [1] 213.934
cor(X,Y)
## [1] 0.706372
cov(X,Y)/(sd(X)*sd(Y))
## [1] 0.706372

Kết quả cho thấy \(X\)\(Y\) tương quan khá tốt với nhau. Đồ thị dưới đây cho phép ta khảo sát dấu của hệ số tương quan.

3.6 Câu hỏi lập trình

3.6.1 Câu 1: Thiên lệch

Xem xét một estimator khác cho \(\mu_Y\): \[\tilde{Y} = \frac{1}{n-1} \sum_{i=1}^n Y_i\]

Chứng minh rằng estimator này bị thiên lệch. Sử dụng set.seed(123). Thực hiện các bước sau:

  1. Định nghĩa Y_tilde
  2. Lấy ngẫu nhiên \(5\) quan sát từ \(N(10,25)\) và tính ước lượng Y_tilde. Lặp lại \(10000\) lần
  3. Vẽ histogram
  4. Vẽ đường trung bình tổng thể
eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiIjIFN0ZXAgMVxuWV90aWxkZSA8LSBmdW5jdGlvbih4KXtcblxufVxuXG5zZXQuc2VlZCgxMjMpXG4jIFN0ZXAgMlxuXG4jIFN0ZXAgM1xuXG4jIFN0ZXAgNCIsInNvbHV0aW9uIjoiWV90aWxkZSA8LSBmdW5jdGlvbih4KXtcbiAgKDEvKGxlbmd0aCh4KS0xKSkqc3VtKHgpXG59XG5zZXQuc2VlZCgxMjMpXG5lc3RfYmlhc2VkIDwtIHJlcGxpY2F0ZShuPTEwMDAwLGV4cHIgPSBZX3RpbGRlKHJub3JtKG49NSwgbWVhbiA9IDEwLCBzZCA9IDUpKSlcbmhpc3QoZXN0X2JpYXNlZClcbmFibGluZSh2PTEwLCBjb2wgPSAncmVkJykifQ==

3.6.2 Câu 2: Tính nhất quán

Estimator ở trên có nhất quán không? Sử dụng set.seed(123)

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiJzZXQuc2VlZCgxMjMpXG4jIFRcdTFlZWIgcGhcdTAwZTJuIHBoXHUxZWQxaSBOKDEwLDI1KSBsXHUxZWFieSBuZ1x1MWVhYnUgbmhpXHUwMGVhbiAxMDAwIHF1YW4gc1x1MDBlMXQsIHRcdTAwZWRuaCBnaVx1MDBlMSB0clx1MWVjYiBZX3RpbGRlLCBsXHUxZWI3cCBsXHUxZWExaSAxMDAwMCBsXHUxZWE3bi5cblxuXG5cbiMgVlx1MWViZCBoaXN0b2dyYW1cblxuIyBUaFx1MDBlYW0gXHUwMTExXHUwMWIwXHUxZWRkbmcgdGhcdTFlYjNuZyBcdTAxMTFcdTFlZTluZyB0aFx1MWVjMyBoaVx1MWVjN24gdHJ1bmcgYlx1MDBlY25oIHRcdTFlZDVuZyB0aFx1MWVjMyBcdTAxMTFcdTAwZmFuZy4iLCJzb2x1dGlvbiI6IllfdGlsZGUgPC0gZnVuY3Rpb24oeCl7XG4gICgxLyhsZW5ndGgoeCktMSkpKnN1bSh4KVxufVxuc2V0LnNlZWQoMTIzKVxuZXN0X2NvbnNpc3RlbnQgPC0gcmVwbGljYXRlKG49MTAwMDAsZXhwciA9IFlfdGlsZGUocm5vcm0obj0xMDAwLCBtZWFuID0gMTAsIHNkID0gNSkpKVxuaGlzdChlc3RfY29uc2lzdGVudClcbmFibGluZSh2ID0gMTAsIGNvbCA9ICdyZWQnKSJ9

3.6.3 Câu 3: Tính hiệu quả

Cho hai estimators như sau: \[\hat{\mu_Y} = \frac{1}{n}\sum_{i=1}^nY_i\]

và: \[\tilde{\mu_Y} = \sum_{i=1}^n a_iY_i\]

trong đó \(b_i\) cho \(n/2\) quan sát đầu tiên có trọng số lớn hơn với \(n\) chẵn.

Giữa hai estimators này, cái nào hiệu quả hơn? Cho biết vector trọng số của \(b_i\) như sau c(rep((1+0.5)/n,n/2),rep((1-0.5)/n,n/2))

eyJsYW5ndWFnZSI6InIiLCJzYW1wbGUiOiJuIDwtIDEwMFxudyA8LSBjKHJlcCgoMSswLjUpL24sbi8yKSxyZXAoKDEtMC41KS9uLG4vMikpXG5cbm11X3RpbGRlIDwtIGZ1bmN0aW9uKHgpe1xuXHQjQlx1MWVhZnQgXHUwMTExXHUxZWE3dSBjb2RlIFx1MWVkZiBcdTAxMTFcdTAwZTJ5OiBcdTAxMTFcdTFlY2JuaCBuZ2hcdTAxMjkgaFx1MDBlMG0gbXVfdGlsZGVcblxuXG5cblx0I0tcdTFlYmZ0IHRoXHUwMGZhYyBjb2RlIFx1MWVkZiBcdTAxMTFcdTAwZTJ5XG59XG5cbnNldC5zZWVkKDEyMylcbiMgQ29kZTogdFx1MWVlYiBwaFx1MDBlMm4gcGhcdTFlZDFpIE4oNSwxMCkgY2hcdTFlY2RuIG1cdTFlYWJ1IHF1eSBtXHUwMGY0IDEwMCwgbFx1MWViN3AgbFx1MWVhMWkgMTAwMDAgbFx1MWVhN25cblxuXG4jIENvZGU6IHRcdTAwZWRuaCBwaFx1MDFiMFx1MDFhMW5nIHNhaSBtXHUxZWFidSBjaG8gdFx1MWVlYm5nIGVzdGltYXRvciIsInNvbHV0aW9uIjoibiA8LSAxMDBcbncgPC0gYyhyZXAoKDErMC41KS9uLG4vMikscmVwKCgxLTAuNSkvbixuLzIpKVxuXG5tdV90aWxkZSA8LSBmdW5jdGlvbih4KXtcbiAgc3VtKHcqeClcbn1cblxuc2V0LnNlZWQoMTIzKVxuZXN0X2JhciA8LSByZXBsaWNhdGUobj0xMDAwMCwgZXhwciA9IG1lYW4ocm5vcm0oMTAwLDUsMTApKSlcbmVzdF90aWxkZSA8LSByZXBsaWNhdGUobj0xMDAwMCwgZXhwciA9IG11X3RpbGRlKHJub3JtKDEwMCw1LDEwKSkpXG5cbnZhcl9iYXIgPC0gdmFyKGVzdF9iYXIpXG52YXJfdGlsZGUgPC0gdmFyKGVzdF90aWxkZSlcbmQgPC0gY2JpbmQodmFyX2Jhcix2YXJfdGlsZGUpXG5jb2xuYW1lcyhkKSA8LSBjKGV4cHJlc3Npb24oaGF0KG11KSksIGV4cHJlc3Npb24odGlsZGUobXUpKSlcbmxpYnJhcnkoa25pdHIpXG5rYWJsZShkKSAlPiVcbiAga2FibGVfc3R5bGluZyhib290c3RyYXBfb3B0aW9ucyA9IFwic3RyaXBlZFwiLCBmdWxsX3dpZHRoID0gRikifQ==