[R 통계분석] 추정

Section 01. 점추정

 

> ## 01) 유효성
> 
> # 두 추정량의 분산 비교
> # Y1 = (X1 + X2 + X3) / 3
> # Y2 = (X1 + 2X2 + 3X3) / 6

> 
> x <- seq(-3, 3, by = 0.01)
> y <- dnorm(x)
> # Y1 : N(1, sqrt(1/3))
> y.1 <- dnorm(x, sd = sqrt(1/3))
> # Y2 : N(1, sqrt(7/18))
> y.2 <- dnorm(x, sd = sqrt(7/18))
> 
> # P(-0.1<X<=0.1)
> pnorm(0.1, sd = sqrt(1/3)) - pnorm(-0.1, sd = sqrt(1/3))
[1] 0.1375098
> pnorm(0.1, sd = sqrt(7/18)) - pnorm(-0.1, sd = sqrt(7/18))
[1] 0.1273999
> # Y1의 분포의 확률값이 더 높아 모평균 주변으로 더 많이 모여 있다고 할 수 있다.
> 
> # 분포 확인
> plot(x, y, type = "l", ylim = c(0, 0.8), axes = F, ylab = "", lwd = 3, col = "yellow")
> lines(x, y.1, col = "red", lwd = 3)
> lines(x, y.2, col = "green", lty = 2, lwd = 3)
> axis(1)
> legend("topright", c("Y1","Y2"), lty=1:2, cex=0.7, col=c("red", "green"))

> # 두 추정량의 평균이 0으로 동일하고, Y1의 분포가 Y2의 분포보다 상대적으로 덜 퍼져있음을 알 수 있다.
> 

> ## 02) 유효성 모의실험
> 
> options(digits = 3)
> set.seed(1)
> # Y1 : mean
> # Y2 : mean.seq
> mean.seq <- function(x) {
+   n <- length(x)
+   sum <- 0
+   n2 <- 0
+   for (i in 1:n) {
+     newx <- i * x[i]
+     sum <- sum + newx
+     n2 <- n2 + i
+   }
+   return(sum / n2)
+ }
> 
> y1 <- rep(NA, 1000)
> y2 <- rep(NA, 1000)
> for (i in 1:1000) {
+   smp <- rnorm(3)
+   y1[i] <- mean(smp)
+   y2[i] <- mean.seq(smp)
+ }
> 
> n1 <- length(y1[(y1 > -0.1) & (y1 < 0.1)])
> n2 <- length(y2[(y2 > -0.1) & (y2 < 0.1)])
> data.frame(mean = mean(y1), var = var(y1), n = n1)
     mean  var   n
1 -0.0042 0.36 134
> data.frame(mean = mean(y2), var = var(y2), n = n2)
     mean   var   n
1 -0.0113 0.427 120
> # -0.1과 0.1 사이의 횟수가 Y1을 사용한 경우가 더 크다.
> 
> par(mfrow = c(1,2))
> hist(y1, probability = T, xlim = c(-2, 2), ylim = c(0, 0.65),
+      main = "(x1+x2+x3)/3", xlab = "", col = "orange", border = "red")
> hist(y2, probability = T, xlim = c(-2, 2), ylim = c(0, 0.65),
+      main = "(1*x1+2*x2+3*x3)/6", xlab = "", col = "orange", border = "red")
> # 미세한 차이지만 Y1의 히스토그램이 Y2보다 평균 주변의 높이가 높음을 확인할 수 있다.

> 

> ## 03) 모비율에 대한 점추정
> 

> # 표본의 크기는 3이고 복원으로 추출하여(주사위를 세번 굴려) 표본비율들의 분포를 구한다.

> 
> library(prob)
> n <- 3
> smps.all <- rolldie(n)
> str(smps.all)
'data.frame': 216 obs. of  3 variables:
 $ X1: int  1 2 3 4 5 6 1 2 3 4 ...
 $ X2: int  1 1 1 1 1 1 2 2 2 2 ...
 $ X3: int  1 1 1 1 1 1 1 1 1 1 ...
> head(smps.all, 3)
  X1 X2 X3
1  1  1  1
2  2  1  1
3  3  1  1
> 
> # 짝수 여부를 확인하는 함수
> is.even <- function(x) return(!x%%2)
> # 모분산을 구하는 함수
> var.p <- function(x) {
+   return( sum((x-mean(x))^2) / length(x) )
+ }
> # 표본 중에서 짝수의 비율을 구하는 함수
> p.even <- function(x, s.size=3) {
+   return(sum(is.even(x)) / s.size)
+ }
> 
> # 행별 짝수의 비율
> phat <- apply(smps.all, 1, p.even)
> 
> # 표본 비율의 기대값
> mean(phat)
[1] 0.5
> # 모집단에서 짝수의 비율
> ( p.p <- 0.5 )
[1] 0.5
> # 표본 비율의 분산
> var.p(phat)
[1] 0.0833
> # 모비율을 이용하여 구한 표본 비율의 분산
> ( p.p*(1-p.p)/3 )
[1] 0.0833
> # 표준오차 = sqrt(표본분산)
> sqrt(var.p(phat))
[1] 0.289
> 

 

Section 02. 구간추정

 

> # 표본평균의 분포에서의 신뢰수준(1-α)
> 
> par(mar=c(0,1,0,1))
> x <- seq(-3, 3, by=0.01)
> y <- dnorm(x)
> plot(x, y, axes=F, type="l", ylim=c(-0.1, 0.5), xlab="", ylab="", col = "blue")
> abline(h=0, col = "blue")
> ll <- qnorm(0.025)
> ul <- qnorm(0.975)
> polygon(c(-3, x[x<ll], ll), c(0, y[x<ll], 0), density=20, col = "red")
> polygon(c(ul, x[x>ul], 3), c(0, y[x>ul], 0), density=20, col = "red", angle=135)
> 
> text(0, 0.2, expression(1-alpha))
> text(-2.0, 0.1, expression(plain(P)(Z<z) == over(alpha, 2)), cex=0.7)
> text(2.0, 0.1, expression(plain(P)(Z>z) == over(alpha, 2)), cex=0.7)
> text(-1.96, -0.02, expression(-z[over(alpha, 2)]), cex=0.8)
> text(1.96, -0.02, expression(z[over(alpha, 2)]), cex=0.8)

> 

> ## 04) 모평균에 대한 95% 신뢰구간
> 
> # 표본정규분포로부터 10개의 표본을 뽑아 95% 신뢰구간을 구하는 것을 100번 반복했을 때, 몇 개의 신뢰구간이 모평균 0을 포함할지 확인해본다.
> 
> set.seed(9)
> n <- 10
> x <- 1:100
> y <- seq(-3, 3, by = 0.01)
> 
> # 각 시행별로 10개씩 표본을 추출하기 위해 행렬 생성
> smps <- matrix(rnorm(n * length(x)), ncol = n)
> 
> # 각 행별로 평균을 구한다.
> xbar <- apply(smps, 1, mean)
> # 표준오차를 구한다. (모집단이 표준정규분포)
> se <- 1 / sqrt(10)
> # 오류확률(α)
> alpha <- 0.05
> # 하한의 z와 상한의 z 사이의 면적(확률)이 (1-α)가 되는 z 값을 구한다.
> z <- qnorm(1 - alpha/2)
> # 신뢰구간(하한값, 상한값)을 구한다.
> ll <- xbar - z * se
> ul <- xbar + z * se
> 
> # 100번의 표본추출별로 신뢰구간을 그린다.
> # 또한 각 표본별 신뢰구간이 모평균을 포함하는지에 따라 선의 색상을 다르게 한다.
> plot(x, type = "n", xlab = "표본추출", ylab = "z",
+      main = "95% Confidence Interval for Population mean",
+      xlim = c(1, 100), ylim = c(-1.5, 1.5), cex.lab = 1.8)
> abline(h=0, col = "red", lty = 2)
> # 각각의 신뢰구간이 평균(0)을 포함하면 하한은 음수, 상한은 양수이므로 곱한 값은 음수가 된다.
> l.c <- rep(NA, length(x))
> l.c <- ifelse(ll * ul > 0, "red", "black")
> # 화살표의 머리 각을 90도(angle=90)로 하고 화살표의 시작점과 끝점에 머리가 생기도록(code=3) 하여
> # 신뢰구간을 표현한다.
> arrows(1:length(x), ll, 1:length(x), ul, code = 3,
+        angle = 90, length = 0.02, col = l.c, lwd = 1.5)

> 

> ## 05) 모평균에 대한 95% 신뢰구간 (모분산을 모를 때)
> 
> # 모회사에서 초등학교 신입생들을 대상으로 모자를 제작하기 위해
> # 부모의 동의를 얻은 10명을 임의 추출하여 머리 둘레를 측정하였다.
> # 모평균에 대한 신뢰구간을 구해본다.
> ci.t <- function(x, alpha=0.05) {
+   n <- length(x)
+   m <- mean(x)
+   s <- sd(x)
+   t <- qt(1 - alpha/2, df = n-1)
+   ll <- m - t * (s/sqrt(n))
+   ul <- m + t * (s/sqrt(n))
+   ci <- c(1-alpha, ll, m, ul)
+   names(ci) <- c("Confidence Level", "Lower limit", "Mean", "Upper limit")
+   return(ci)
+ }
> 
> smp <- c(520, 498, 481, 512, 515, 542, 520, 518, 527, 526)
> # 95% 신뢰구간
> ci.t(smp)
Confidence Level      Lower limit             Mean      Upper limit 
            0.95           503.98           515.90           527.82 
> # 90% 신뢰구간
> ci.t(smp, 0.1)
Confidence Level      Lower limit             Mean      Upper limit 
             0.9            506.2            515.9            525.6 
> # 앞의 95% 신뢰구간보다 그 폭이 줄어든 것을 확인할 수 있다.

 

 

 

 

 

 

출처 : 이윤환, ⌜제대로 알고 쓰는 R 통계분석⌟, 한빛아카데미, 2020