몬테카를로 시뮬레이션(MC Simulation) 시뮬레이션을 이용하는 이유는 수리적으로 계산이 불가능하거나 너무 어려운 경우에 시뮬레이션을 이용하면 비교적 쉽게 계산이 가능하다.

표준정규분포에서 100개의 랜덤 변수를 발생시켜라

x <- rnorm(100, mean=0, sd=1)
summary(x) #평균이 0에 가깝지만 0은 아니다.
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.34800 -0.67040  0.08683  0.00861  0.62000  2.14700
x <- rnorm(10000, 0, 1)
summary(x) #n에 커지면서 평균0에 훨씬 밀접한 결과가 나왔다.
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.92000 -0.64860  0.02608  0.02106  0.68610  3.22300
hist(x)

랜덤변수 x가 평균이 10이고 표준편차가 5인 정규분포를 따를때 (x^2 + 2x + log(x^2 + 3)) / sqrt(5abs(x) + exp(x)) 의 평균과 분산을 계산하기

x <- rnorm(10000, 10,5)
hx <- (x^2 + 2*x + log(x^2 + 3)) / sqrt(5*abs(x) + exp(x))
mean(hx)
## [1] 1.124611
var(hx)
## [1] 1.021346

표본수 n을 증가하면 MC 추정치의 분산이 줄어들고 좀 더 정확한 추정이 가능하며 많은 경우에 이런 MC추정치 계산은 한번만 하지 않고 여러번 반복

sim.n <- 1000
myhx <- rep(0,sim.n)
for(i in 1:sim.n){
  x <- rnorm(10000, 10,5)
  hx <- (x^2 + 2*x + log(x^2 + 3)) / sqrt(5*abs(x) + exp(x))
  myhx[i] <- mean(hx)
}

mean(myhx)
## [1] 1.130038
var(myhx)
## [1] 8.812085e-05

중심극한 정리(Central Limit Theorem) 통계학에서 가장 중요한 정리중에 하나이다. 표본의 수가 커지면 표본평균의 분포가 모집단의 분포에 상관없이 정규분포로 근사한다. - 시뮬레이션을 통해서 모집단의 분포가 정규분포가 아니더라도 표본평균의 분포가 정말 정규분포로 근사되는지 확인

카이제곱분포를 따르는 pdf및 표본수 증가에 따른 표본평균의 분포의 변화

x <- seq(0,15, by=0.1)
plot(x,dchisq(x,3),type="l", main="pdf of chisq(3)")

표본크기가 5,10,20,50인 4가지 경우에 표본평균의 분포가 어떻게 되는지 비교하기

sim.n <- 1000
sam.n <- c(5,10,20,50)
x.mean <- matrix(0,sim.n,4)
for(i in 1:4){
  for(j in 1:sim.n){
    x <- rchisq(sam.n[i],3)
    x.mean[j,i] <- mean(x)
  }
}
par(mfrow=c(2,2))
hist(x.mean[,1], main="n=5")
hist(x.mean[,2], main="n=10")
hist(x.mean[,3], main="n=20")
hist(x.mean[,4], main="n=50")

표본수가 커짐에 따라 정규분포에 가까와지는 모습을 확인할 수 있다.

수치적분(Numerical Integration) 기본적인 아이디어는 적분 구간을 나눈 다음에 구간 내에서 적분해야하는 함수를 다항식 함수로 근사한 다음에 적분하는것. R에서는 integrate 함수 사용

integrate(dnorm, -1.96,1.96)
## 0.9500042 with absolute error < 1e-11
integrate(dnorm, -Inf,Inf)
## 1 with absolute error < 9.4e-05

(x^2 + 2x + log(x^2 + 3)) / sqrt(5abs(x) + exp(x)) 적분

mytestfp1 <- function(x){
  hx <- (x^2 + 2*x + log(x^2 + 3)) / sqrt(5*abs(x) + exp(x))  
  res <- hx * dnorm(x,10,5)
  return(res)
}

integrate(mytestfp1, -Inf, Inf)
## 1.129548 with absolute error < 7.1e-06

MC 시뮬레이션 결과와 소수점 4째자리까지 일치. 일반적으로 수치적분이 더 정확한 값을 주고 (오류율까지 제공) 계산도 더 빠른 경우가 많음

최적화 기법 통계학이나 수학에서 최적화란 어떤 함수의 최대값이나 최소값을 찾는것. optimize 함수 이용

f(x) = (x-2)^2 + 3을 최소화 하는 값과 그때의 함수의 값을 찾아라.

myfn1 <- function(x){
  res <- (x-2)^2 + 3
  return(res)
}
optimize(myfn1,c(-10,10))
## $minimum
## [1] 2
## 
## $objective
## [1] 3

결과값에서 minimum은 x의 최소값. objective는 x가 최소값일때 함수의 최소값을 의미.

f(x) = log(x) / (1+x)의 함수츼 최대값을 구하라.

myfn2 <- function(x){
  res <- log(x) / (1+x)
  return (res)
}
optimize(myfn2,c(0,20),maximum = T)
## $maximum
## [1] 3.591121
## 
## $objective
## [1] 0.2784645
y <- myfn2(1:20)
par(mfrow=c(1,1))
plot(y,type="l")