9.6 Joint Confidence Region

9.6.1 Drawing bivariate normal using distribution function and sampled data

다음은 bivariate normal distribution에서 95% 신뢰영역(confidence region)을 그리는 R script이다. 이것은 joint confidence region이라고도 한다.

ci = 0.95
npoints = 100
dimR = 2 # plotting dimension
npara = 5 # mu: 2, mCov: 3

mu0 = c(0, 0)
mCov0 = matrix(c(1, 0.5, 0.5, 1), nrow=2)
eg0 = eigen(mCov0)
alpha0 = atan(eg0$vectors[2, 1]/eg0$vectors[1, 1])
radius0 = sqrt(qchisq(ci, dimR)*eg0$values) # from normal distribution

r1 = ellipse(mu0, radius0, alpha0, col="red", asp=1)
points(mu0[1], mu0[2], pch="*", col="red")
abline(a=0, b=1, lty=2)
abline(a=0, b=-1, lty=2)
abline(h=0, lty=3)
abline(v=0, lty=3)

Figure 9.5: Joint confidence region of bivariate normal

이제 위의 이론적인 그림에 난수를 이용해서 구한 점들과 이로 추정된 평균 벡터, 분산-공분산 행렬을 이용하여 그려보자.

require(MASS)
Data = mvrnorm(npoints, mu0, mCov0)
head(Data, 10)

            [,1]        [,2]
 [1,] -0.1964644 -1.54205679
 [2,] -0.9954189 -0.37614152
 [3,]  0.7990320 -0.04081758
 [4,]  1.4542735 -0.92468583
 [5,]  0.1334372  0.60586253
 [6,]  0.6128279  0.42606646
 [7,] -0.6776387 -1.16541165
 [8,] -0.8820948 -0.64158363
 [9,]  0.4946249 -1.23449997
[10,] -0.1812580  1.20569189

mu = colMeans(Data) ; mu

[1] -0.03290263 -0.13541921

mCov = cov(Data) ; mCov

          [,1]      [,2]
[1,] 0.8349144 0.3025257
[2,] 0.3025257 0.8042766

eg = eigen(mCov) # sorted as descening order
alpha = atan(eg$vectors[2,1]/eg$vectors[1,1]) ; alpha*180/pi

[1] 43.5506

radius = sqrt(dimR*qf(ci, dimR, npoints - npara)*eg$values) # from sample data

ellipse(mu0, radius0, alpha0, col="red", asp=1)
ellipse(mu, radius, alpha, asp=1, add=T)
points(Data)
points(mu0[1], mu0[2], pch="*", col="red")
points(mu[1], mu[2], pch="+")
abline(a=0, b=1, lty=2, col="red")
abline(a=0, b=-1, lty=2, col="red")
abline(h=0, lty=3)
abline(v=0, lty=3)

Figure 9.6: Simulated data and confidence region

상당히 유사하게 겹치는 것을 볼 수 있다.

추정된 평균벡터와 분산공분산 행렬을 이용해서 각 축에서 바라본 95% 신뢰구간을 구하면 다음과 같다.

ellipRange(mu0, radius0, alpha0) # True value

        xrange    yrange
[1,] -2.447747 -2.447747
[2,]  2.447747  2.447747

ellipRange(mu, radius, alpha) # Estimated value

        xrange    yrange
[1,] -2.305230 -2.365665
[2,]  2.239425  2.094827