메뉴 건너뛰기

웹에서 하는 R 통계

함수곡선 아래의 면적 구하기

2016.12.22 13:49

cardiomoon 조회 수:304

R에서 함수곡선 아래의 면적을 구하려면 어떻게 할까? 예를 들어 y=sin(x) 함수 곡선에서 π4부터 34×π 까지의 면적은 얼마일까?

이 면적을 수동으로 계산하지 않더라도 R의 내장함수인 integrate()를 사용하면 쉽게 구할 수 있다.

integrate(sin,pi/4,pi*3/4)
## 1.414214 with absolute error < 1.6e-14

이번 장에서는 integrate()함수를 쓰지 않고 Newton-Cotes Formula를 이용하여 면적을 구하는 방법을 다룬다. 이 글의 전반부는 Hadley Wickam의 Advanced-R에 나온 내용에서 인용하였다.

10.5 Case study: numerical integration

Calculation with midpoint and trapezoid

값을 구하고자 하는 범위의 시작을 a, 끝을 b라고 하자.

  1. 먼저 midpoint는 a와 b의 중간값인 a+b2에 대한 함수 값을 높이로 하고 b-a를 밑변으로 하는 사각형의 면적으로 구하는 방법이다.
    midpoint <- function(f, a, b) {
      (b - a) * f((a + b) / 2)
    }
    
    (result=midpoint(sin,pi/4,pi*3/4))
    [1] 1.570796

    1. trapezoid는 a의 함수값 f(a)와 b의 함수값 f(b)의 평균을 높이로 하는 사각형의 높이로 구하는 방법이다.
    trapezoid1<-function(f,a,b){
        (b-a)*(f(a)+f(b))/2
    }
    
    (result=trapezoid1(sin,pi/4,pi*3/4))
    ## [1] 1.110721

    1. 보다 정밀한 방법으로 a부터 b까지를 몇 개의 구간으로 나누고 각 구간에 해당하는 사각형의 면적을 합해서 근사치를 만드는 방법을 생각해 볼 수 있다. 여러 가지 방법들이 알려져 있는데 이러한 방법들을 Issac Newton과 Roger Cotes의 이름을 따서 Newton-Cotes rules라고 이야기 하며 널리 알려진 formula는 다음과 같은 것들이 있다.

    Newton-Cotes Rules

    (https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas)

    Simpson의 방법은 전체구간을 3개의 사각형으로 나누어 계산하는 방법이다.

    Boole의 방법은 전체를 90개의 구간으로 나누고 5개의 사각형으로 면적을 측정한다.

    traepzoid방법도 전체를 10개의 구간으로 나누어 적용한다면 보다 정밀해진다.

    trapezoid_composite<-function(f,a,b,n=10){
        
        points <- seq(a,b,length=n+1)
        width <- (b-a)/n
        
        area<-0
        for(i in seq_len(n)){
            area<-area+width*f((points[i]+points[i+1])/2)
        }
        area
    }
    
    (result=trapezoid_composite(sin,0,pi/2))
    ## [1] 1.001029

    Function for general Newton-Cotes Rule:

    Hadley가 쓴 advanced-R 책에는 이와 같은 Newton-Cotes Rule을 일반화하는 함수를 소개하고 있다.

    composite=function(f,a,b,n=10,rule){
        points<-seq(a,b,length=n+1)
        area<-0
        result=list()
        for(i in seq_len(n)){
            result[[i]]<-rule(f,points[i],points[i+1])
            area<-area+rule(f,points[i],points[i+1])
        }
        area
    }
    
    newton_cotes<-function(coef,open=FALSE){
        n<-length(coef)+open
        
        result=function(f,a,b){
            pos<- function(i) a+i*(b-a)/n
            points<-pos(seq.int(0,length(coef)-1))
            
            (b-a)/sum(coef)*sum(f(points)*coef)
            
        }
    }
    
    trapezoid2<-newton_cotes(c(1,1))

    Recalcuation with generalized function

    이 함수를 이용하여 0부터 pi까지의 값을 계산해보면 다음과 같고 이 값은 trapezoid1(sin,0,pi)의 값은 같아야 한다.

    # It should be equal to trapezoid(sin,pi/4,pi*3/4)
    composite(sin,0,pi,n=1,rule=trapezoid2)
    [1] 1.570796
    trapezoid1(sin,0,pi)
    [1] 1.923671e-16

    왜 값이 다를까?

    값이 왜 다른지, 어디에서 에러가 나는지 보기위해 위키의 내용을 다시 리뷰해보았다.

    (https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas)

    Newton-Cote rule에서 자유도 n 은 계수의 갯수-1로 계산해야 한다. 하지만 Hadley도 신이 아니라 사람인지라 실수한 것으로 생각된다. 함수에 다음과 같이 고치니 바라는 대로 동작하였다.

    newton_cotes1<-function(coef,open=FALSE){
        n<-length(coef)-1+open
        
        result=function(f,a,b){
            pos<- function(i) a+i*(b-a)/n
            points<-pos(seq.int(0,length(coef)-1))
            
            (b-a)/sum(coef)*sum(f(points)*coef)
            
        }
    }
    trapezoid3<-newton_cotes1(c(1,1))
    composite(sin,0,pi,n=1,rule=trapezoid3)
    [1] 1.923671e-16

    Plot for the numerical integration

    이 함수을 이용하여 plot을 그리기 위해 함수의 내용을 약간 변경하여 다음과 같은 함수를 제작하였다.

    newton_cotes<-function(coef,open=FALSE){
        n<-length(coef)-1+open
        
        result=function(f,a,b){
            pos<- function(i) a+i*(b-a)/n
            points<-pos(seq.int(0,length(coef)-1))
            
            area=(b-a)/sum(coef)*sum(f(points)*coef)
            list(a=a,b=b,y=f(points),coef=coef,area=area)
        }
        result
    }
    
    composite=function(f,a,b,n=10,rule){
        points<-seq(a,b,length=n+1)
        area<-0
        result=list()
        for(i in seq_len(n)){
            result[[i]]<-rule(f,points[i],points[i+1])
            area<-area+rule(f,points[i],points[i+1])$area
        }
        list(result,f=f,area=area)
    }
    require(ggplot2)
    plot_quadrature=function(result,color="red",fill="blue",alpha=0.2,size=0.2){
        test=result$f
        cadd <- function(x) Reduce("+", x, accumulate = TRUE)
        p<-ggplot(data.frame(x=c(0,pi)),aes(x))+stat_function(fun=test)
     
        for(i in 1:length(result[[1]])){
            subdata=result[[1]][[i]]
            coef=subdata$coef
            x=seq(subdata$a,subdata$b,length=(sum(coef)+1))
            h=subdata$y
            start=c(1,cadd(coef)+1)
            for(j in 1:length(coef)){
                p<-p+geom_rect(xmin=x[start[j]],xmax=x[start[j+1]],
                               ymin=0,ymax=h[j],colour=color,fill=fill,alpha=alpha,size=size)
            }
        }
        p
    }

    Plots for general Newton-Cotes Rule

    이 함수를 이용하여 0부터 pi까지 사인함수 곡선의 면적을 측정하는 plot을 그리면 다음과 같다.

    boole<-newton_cotes(c(7,32,12,32,7))
    milne<-newton_cotes(c(2,-1,2),open=TRUE)
    trapezoid<-newton_cotes(c(1,1))
    simpson<-newton_cotes(c(1,4,1))
    simpson8<-newton_cotes(c(1,3,3,1))
    
    rules=list(trapezoid,simpson,simpson8,boole)
    
    result=lapply(rules,function(x) composite(sin,0,pi,n=10,rule=x))
    plot_quadrature(result[[1]])+ggtitle("Rule : Trapezoid, n=10")+
        geom_text(x=pi/2,y=sin(pi/2)/2,label=result[[1]]$area)

    plot_quadrature(result[[2]])+ggtitle("Rule : Simpson's , n=10")+
        geom_text(x=pi/2,y=sin(pi/2)/2,label=result[[2]]$area)

    plot_quadrature(result[[3]])+ggtitle("Rule : Simpson's 3/8, n=10")+
        geom_text(x=pi/2,y=sin(pi/2)/2,label=result[[3]]$area)