[R] 결합적률생성함수(JMGF)와 수치계산,기각 샘플링 시뮬레이션

GPT의 도움을 빌려서 각 확률변수의 적률생성함수의 결과를 구해두었다. 

### 1. Closed form MGF 정의
Mx1_cf <- function(t){
    if(any(t >= 2)) stop("M_X1(t)는 t < 2 에서만 수렴.")
    2 / (2 - t)
}

Mx2_cf <- function(t){
    if(any(t >= 1)) stop("M_X2(t)는 t < 1 에서만 수렴.")
    2 / ((1 - t) * (2 - t))
}

### 2. 수치적분을 이용한 계산
Mx1_num <- function(t, upper2 = 50) {
    if (t >= 2)
    
    inner_scalar <- function(x2) { #한번 x1에 대해서 적분한 결과는 x2만의 함수로 남는다.
        f_x1 <- function(x1) 2 * exp(-(x1 + x2)) * exp(t * x1) #우선 x1에 대해서 적분하기 위해 함수를 이중으로 이용.
        integrate(f_x1, lower = 0, upper = x2)$value
    }
    
    inner_vec <- Vectorize(inner_scalar, vectorize.args = "x2")
    #벡터로 만드는데 있어 벡터인자를 x2만 쓰고 t같은 것은 스칼라로 받음.
    #벡터화 안시키면 스칼라로 함수가 받아져서 오류발생
    
    # 유한 상한 upper2 사용
    integrate(inner_vec, lower = 0, upper = upper2)$value
}

Mx2_num <- function(t, upper2 = 50) {
    if (t >= 1) stop
    
    inner_scalar <- function(x1) { 
        f_x2 <- function(x2) 2 * exp(-(x1 + x2)) * exp(t * x1)
        integrate(f_x2, lower = 0, upper = x1)$value
    }
    
    inner_vec <- Vectorize(inner_scalar, vectorize.args = "x1")
    #벡터로 만드는데 있어 벡터인자를 x1만 쓰고 t같은 것은 스칼라로 받음.
    
    # 유한 상한 upper2 사용
    integrate(inner_vec, lower = 0, upper = upper2)$value
}
Mx2_num(0.3, upper2 = 50)

### 4. 비교 실행
t1 <- 0.5  # for X1
t2 <- 0.3  # for X2

cat("=== M_X1(t=0.5) ===\n"
    ,"닫힌식      :", Mx1_cf(t1), "\n",
    "수치적분    :", Mx1_num(t1), "\n")
    
cat("=== M_X2(t=0.3) ===\n",
"닫힌식      :", Mx2_cf(t2), "\n",
"수치적분    :", Mx2_num(t2), "\n")

 

> cat("=== M_X1(t=0.5) ===\n"
+     ,"닫힌식      :", Mx1_cf(t1), "\n",
+     "수치적분    :", Mx1_num(t1), "\n")
=== M_X1(t=0.5) ===
 닫힌식      : 1.333333 
 수치적분    : 1.333333
 
 > cat("=== M_X2(t=0.3) ===\n",
+     "닫힌식      :", Mx2_cf(t2), "\n",
+     "수치적분    :", Mx2_num(t2), "\n")
=== M_X2(t=0.3) ===
 닫힌식      : 1.680672 
 수치적분    : 1.680672

 

시뮬레이션은 기각샘플링을 이용하겠다. 아래내용부터는 상당히 딥한 내용이다. 

지식적한계에 부딪혀서 이론에 대한 내용은 내려놓고 코드구현으로 마치겠다.

아래 알고리즘표와 구글링을 이용해서 구현하였다.

출처: https://en.wikipedia.org/wiki/Rejection_sampling

u : 균일분포, 부등식을 체크해서 부등호가 성립하면 유지, 아니면 거절하는 방식의 알고리즘이다.

여기서 f(y)를 우리의 목표인 타겟분포(Target dist.), g(y)를 우리가 생각해서 제안한 분포(=Proposal dist.)라고 부르겠다.

또한, M값은

M값이 뭐든간에 답은 나오지만 시간효율성이 안좋다는 얘기가 있다. 제안분포에 따라 M값이 달라지고 효율성만 달라지지 이에 대한 일치성,불편성에 대한 문제는 없다. 이에 대한 의구심이 든다면 아래를 확인하자.

https://en.wikipedia.org/wiki/Rejection_sampling

Rejection sampling - Wikipedia

알고리즘 구현방식은 위의 내용으로 진행하겠다.

rejection_sampler_general <- function(N = 1e5) {
    # 출력용 빈 행렬(샘플 저장)
    X <- matrix(NA, nrow = N, ncol = 2)
    n <- 0        # 현재까지 받은 샘플 수
    while (n < N) {
        # 3. 제안분포 q에 대한 값
        x0 <- rexp(2, rate = 1)        # 독립 Exp(1) 두 개 → (x1, x2)
        
        # 4. u ~ U(0,1)
        u <- runif(1)
        
        # 5. 수락조건: u < f(x0)/(M * g(x0))
        #    f(x0) = 2 * exp(-(x1 + x2)) * I(x1 < x2)
        #    g(x0) =        exp(-x1 - x2)
        #    따라서 f/(M g) = [2 exp(-x1-x2) I(x1<x2)] / [2 exp(-x1-x2)] = I(x1 < x2)
        if (u < ( (x0[1] < x0[2]) * 1 )) {
            # 6. 수락 -> 저장
            n <- n + 1
            X[n, ] <- x0
        }
        
      
    }
    return(X)
}


set.seed(4121)
samples <- rejection_sampler_general(N = 1e5)

# X1에 대한 MGF 추정 (t=0.5)
t <- 0.5
mgf_hat <- mean(exp(t * samples[,1]))
cat(sprintf("Estimated MGF of X1 t=%.2f: %.5f\n", t, mgf_hat))

 

> cat(sprintf("Estimated MGF of X1 t=%.2f: %.5f\n", t, mgf_hat))
Estimated MGF of X1 t=0.50: 1.33443