Sheldon
November 29, 2015
이번에는 R을 가지고 기본적인 시뮬레이션을 하는 방법을 소개하겠습니다. 앞으로 다룰 확률과 통계에 대한 내용엔 결국 random number(난수)를 가지고 시뮬레이션을 하는 것이 이해에 도움이 됩니다. 시뮬레이션을 하려면 먼저 난수(random number)를 생성해야 됩니다. 난수를 생성하는 방법은 많이 있지만 여기서는 runif()를 소개할게요. runif는 uniform distribution을 따르는 random numbers를 생성하는 명령어입니다. Uniform distribution은 다음에 확률에 대해 다룰 때 자세히 언급할게요. runif(10)이라고 하면 0에서 1사이의 숫자 중 랜덤하게 10개를 골라서 출력하게 됩니다.
runif(10)## [1] 0.13894121 0.09478444 0.56983378 0.74138869 0.45956052 0.65444516
## [7] 0.46554138 0.96851162 0.45030918 0.75814276
runif(5)## [1] 0.4717690 0.9361607 0.6793190 0.7472224 0.8143992
runif(20, min=5, max=10)## [1] 6.934429 7.007927 7.635405 8.733614 8.013959 8.562446 6.772075
## [8] 8.961107 5.509713 6.067835 6.598740 5.122696 5.763497 5.433184
## [15] 7.047651 8.696878 7.405097 8.665665 6.887924 6.822792
runif(30, min=0, max=100)## [1] 96.11735 90.49486 42.64492 85.77272 23.34318 18.95856 39.43801
## [8] 16.64120 79.34150 64.55650 61.01691 80.27949 28.07886 8.89806
## [15] 39.38371 24.47642 36.64586 62.24883 85.11644 23.18318 89.06437
## [22] 99.75858 37.46490 81.85986 41.29159 41.39004 71.12814 55.93925
## [29] 74.07582 73.09222
그럼 한 번 이것을 이용해서 동전을 n번 던지는 사건을 시뮬레이션 해 봅시다. 이 때, 만약 0에서 1까지의 runif결과가 0.5보다 크면 앞면, 0.5보다 작거나 크면 뒷면이 되는 것으로 사건을 정의합시다. 그리고 앞면의 경우는 1, 뒷면의 경우는 0으로 나타나기로 했다 합시다. 그럼 동전 던지기 n번을 실행하는 함수 simcoin을 만들어 봅시다.
simcoin<-function(n){
result<-numeric()
for (i in 1:n){
toss<-runif(1)
if (toss<0.5){
result<-c(result,1)
}
else{
result<-c(result,0)
}
}
return (result)
}
simcoin(5) #가상으로 동전을 5번 던졌을 때의 결과입니다.## [1] 1 1 0 1 1
simcoin(10) #가상으로 동전을 10번 던졌을 때의 결과입니다.## [1] 1 1 0 1 0 1 0 1 0 0
mean(simcoin(10))#이렇게 하면 수학적으로 동전을 10번 던졌을 때 앞면이 나올 확률을 알 수 있습니다.## [1] 0.3
mean(simcoin(100))## [1] 0.45
mean(simcoin(1000))## [1] 0.502
mean(simcoin(10000))## [1] 0.4905
동전을 여러 번 던질수록 앞면이 나올 확률은 가상이나마 0.5에 가까워집니다. 이번에는 좀 더 복잡하게 주사위 2개를 던져 나온 숫자의 합을 구하는 시뮬레이션을 해 봅시다. 실행하기 전에 시뮬레이션에 유용한 또 다른 함수 sample()을 소개하겠습니다. 앞에서 말했지만 이렇게 새로운 함수를 알아갈 때는 help()를 눌러서 꼭 찾아보세요. sample은 기본적으로는 벡터를 입력변수로 삼습니다. 그리고 그 벡터안의 구성 요소를 샘플을 추출하듯이 랜덤하게 뽑는 것이지요. sample()의 인수 중 replace는 기본적으로는 FALSE를 기본값으로 가집니다. 그러면 한 번 뽑힌 것은 다시 뽑히지 않죠. 하지만, replace=TRUE로 하면 한 번 뽑힌 원소도 그 다음 추출에 다시 뽑힐 수 있습니다. size는 몇 개나 sample을 추출할 지를 결정합니다. 그런 한 번 실제로 함수를 만져서 익숙해져 봅시다.
x<-1:10
sample(x, size=10)## [1] 1 4 5 8 3 10 6 2 7 9
sample(x, size=10, replace=TRUE)## [1] 9 3 6 10 6 3 9 4 3 1
sample(x, size=30, replace=TRUE)## [1] 4 6 9 5 10 7 4 1 8 4 1 2 3 7 3 6 3 3 1 10 7 10 5
## [24] 2 6 3 6 6 3 1
우리는 이 함수를 이용해서 주사위를 굴릴 것입니다. 만약 우리가 sample(1:6,size=1,replace=TRUE)를 시행한다면 이것은 눈금이 1에서 6까지 있고 각각의 확률이 동일한 주사위를 1번 굴린 것과 동일한 결과를 낼 것입니다. 만약 size를 10으로 한 다면 10번 굴린 것이지요.
sample(1:6, size=1, replace=TRUE)## [1] 4
sample(1:6, size=10, replace=TRUE)## [1] 2 4 6 4 3 3 5 5 5 4
그러면 주사위 2개를 던져 나온 숫자의 합을 구하는 시뮬레이션을 해 봅시다.
sumdice<-function(n){
dice1<-sample(1:6,size=n,replace=TRUE)
dice2<-sample(1:6,size=n,replace=TRUE)
result<-dice1+dice2
return(result)
}여기서 주사위 합이 7이 될 확률은 어떻게 구할까요? 논리 연산자의 산술을 이용하면 다음과 같이하여 주사위 합이 7이 될 확률을 구할 수 있습니다.
y<-sumdice(100)
mean(y==7)## [1] 0.18
이 수치가 수학적으로 계산했을 때 예측할 수 있는 기대값과 맞아 떨어지나요? 잘 모르겠네요. 하지만 시뮬레이션의 횟수를 늘리면 거의 그 값에 근접해지겠죠.
y<-sumdice(10000)
mean(y==7)## [1] 0.1615
굳이 합이 7이 될 확률뿐만 아니라 이 상황에서 여러가지 다른 확률들도 시뮬레이션을 이용해서 구할 수 있습니다. 해당 코드가 어떤 확률을 의미하는 지 한 번 생각해 보세요.
y<-sumdice(10000)
mean(y==3)## [1] 0.0589
mean(y>4)## [1] 0.8295
mean(y<5 & y==4)## [1] 0.0834