9.5 타원 그리기

R에서 bivariate normal distribution을 시각적으로 표현하고 싶을 때는 먼저 타원을 그릴 수 있어야 한다.

다음과 같이 개발중인 math package내의 ellipse라는 함수를 이용하여 보자.

math package는 아직 CRAN에 올라가 있지 않으므로 다음과 같이 설치해야 한다.

install.packages("math", repos="http://r.acr.kr")

require(math)

ellipse = function(center=c(0, 0), radius=c(2, 1), alpha=0, npoints=100, add=F, ...)
{
  theta = seq(0, 2*pi, length=npoints + 1)
  x0 = radius[1]*cos(theta)
  y0 = radius[2]*sin(theta)
  x = x0*cos(alpha) - y0*sin(alpha) + center[1]
  y = x0*sin(alpha) + y0*cos(alpha) + center[2]
  LongAxis = max(radius)
  xlm = center[1] + c(-1, 1)*LongAxis 
  ylm = center[2] + c(-1, 1)*LongAxis
  if (add) {
    lines(x, y, xlim=xlm, ylim=ylm, ...)    
  } else {
    plot(x, y, xlim=xlm, ylim=ylm, type="l", ...)
  }
  invisible(cbind(x, y))
}

위 함수는 타원의 중심(두 중심의 중점)의 좌표, 반지름 두 개(긴 것과 짧은 것, 장반경과 단반경), 중심축의 각도(두 중심을 이은 직선이 x축과 이루는 각도, radian 단위)와 타원을 그리는데 사용할 점의 갯수를 입력 인자(argument)로 받는다.

Default는 중심(0,0)에 반지름(2, 1)이고, 중심축이 수평(0 radian)인 타원을 그려준다.

아래에서 asp=1 옵션을 준 이유는 가로세로가 정확한 scale로 그려지도록 하기 위해서이다.

ellipse(asp=1)
points(0, 0, pch="*")
abline(h = 0, lty=3)
abline(v = 0, lty=3)

Figure 9.2: An example of ellipse

장반경과 단반경을 3과 2로 늘리고, 45도 (\(pi/4\) radian)로 회전한 타원은 다음과 같이 그릴 수 있다.

ellipse(c(0, 0), c(3, 2), pi/4, asp=1)
points(0, 0, pch="*")
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.3: Rotation of ellipse

만약 중심의 위치도 (1,1) 좌표로 옮긴다면 다음과 같이 될 것이다.

ellipse(c(1, 1), c(3, 2), pi/4, asp=1)
points(1, 1, pch="*")
abline(a = 0, b = 1, lty=2)
abline(a = 2, b = -1, lty=2)
abline(h = 0, lty=3)
abline(v = 0, lty=3)

Figure 9.4: Rotation and shifting of ellipse

중심축이 회전한 타원을 그리게 되면, x축에서 바라본 range와 y축에서 바라본 range가 장반경과 단반경의 정확한 2배가 되지 않는다. 따라서 이를 구하는 함수는 다음과 같다.

ellipRange = function(center=c(0, 0), radius=c(2, 1), alpha=0)
{
  x0max = sqrt(radius[1]^2*cos(alpha)^2 + radius[2]^2*sin(alpha)^2) 
  y0max = sqrt(radius[1]^2*sin(alpha)^2 + radius[2]^2*cos(alpha)^2)
  xrange = center[1] + c(-1, 1)*x0max
  yrange = center[2] + c(-1, 1)*y0max
  return(cbind(xrange, yrange))
}

첫번째로 그린 타원은 축이 수평이므로 장반경과 단반경 만으로 바로 구할 수 있다.

ellipRange()

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

하지만, 회전한 경우에는 다음과 같이 나오게 된다.

ellipRange(c(0, 0), c(3, 2), pi/4)

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

ellipRange(c(1, 1), c(3, 2), pi/4)

       xrange   yrange
[1,] -1.54951 -1.54951
[2,]  3.54951  3.54951

ellipRange(c(1, 1), c(3, 2), pi/6)

        xrange    yrange
[1,] -1.783882 -1.291288
[2,]  3.783882  3.291288