정적분의 계산 (2007년 9월 평가원)

(문제) $\int_0^2 |x^2(x-1)| dx $ 의 값은?

(1) $\frac{3}{2}$

(2) $2$

(3) $\frac{5}{2}$

(4) $3$

(5) $\frac{7}{2}$

(풀이)

x=seq(from=0, to=2, by=0.01)
y=abs(x**2 * (x-1))

plot(x,y,type='l')

- 아이디어: (1) 가로가 2 세로가 4인 작사각형에 임의의 점을 뿌린다. (2) 전체점의 수와 함수 아래에 위치한 점의 갯수를 센다. (3) 점의 갯수를 바탕으로 넓이를 추론한다.

xx=runif(10000) # 0~1 사이의 점을 랜덤으로 10000개 뽑음

xx = xx*2 # 0~2사이의 값을 랜덤으로 10000개 뽑은것과 같은 효과

yy=runif(10000)*4 # 0~4 사의의 점을 랜덤으로 10000개 뽑음

plot(xx,yy)

가로가 2, 세로가 4인 직사각형에 무작위로 점이 흩뿌려져 있다.

plot(xx,yy)
lines(x,y,col='red',lwd=3)

붉은선 아래의 점들이 몇개일까?

test = function(xx,yy){
    yy < abs(xx**2 * (xx-1))
}

print(c(xx[1],yy[1])) ## xx[1], yy[1] 의 값 = 무작위로 찍힌 값 
print(abs(xx[1]**2 *(xx[1]-1))) ## xx[1]에 해당하는 붉은점의 값 
test(xx[1],yy[1])

[1] 0.830618 2.372124
[1] 0.1168611

[1] FALSE

plot(xx,yy,col='gray')
lines(x,y,col='red',lwd=3)
points(xx[1],yy[1],col='blue')
points(xx[1],abs(xx[1]**2 * (xx[1]-1)),col='red')

tst = c()
for (i in 1:10000) tst[i] = test(xx[i],yy[i])

head(tst)

[1] FALSE FALSE FALSE FALSE FALSE FALSE

plot(xx,yy,col='gray')
lines(x,y,col='red',lwd=3)
points(xx[tst],yy[tst],col='red')

붉은점들의 갯수는? tst==TRUE인 점들의 수

sum(tst)

[1] 1829

전체 10000개 중에서 1829개정도가 붉은점임. 따라서

$$\frac{s}{8} \approx \frac{1829}{10000}$$

이므로 구하는 넓이는 대략적으로

1829/10000 * 8

[1] 1.4632

따라서 답은 1번 3/2 이다.

- 이와 같은 방법을 몬테카를로 방법이라고 한다.

- 꽤 복잡한 함수에 대하여서도 정적분의 값을 근사적으로 계산해낼 수 있다는 장점이 있음.

(풀이2)

- 사실 그냥 아래와 같이 풀어도 된다.

sum(yy  <  abs(xx^2 * (xx-1)) )

[1] 1829

랜덤변수

랜덤변수를 생성하는 방법

예제1: 균등분포

- 목표: 구간 [0,1]에서 하나의 임의의 값을 뽑는것 (이때 구간 [0,1]에서 각 점이 뽑힐 확률은 동일)

runif(1)

[1] 0.7638336

- 목표: 구간 [0,1]에서 10개의 값을 뽑고 싶다면?

runif(10)

 [1] 0.43916197 0.58808151 0.01061738 0.57311258 0.62008149 0.79359762
 [7] 0.43567517 0.06111325 0.66865770 0.98463477

- 10000개의 값을 뽑아보고 평균이 0.5근처인지 체크해보자.

runif(10000) %>% mean

[1] 0.5001962

예제2: 균등분포

- 목표: 구간 [0,2]에서 하나의 값을 임의로 뽑고 싶다. (이때 구간 [0,2]에서 각 점이 뽑힐 확률은 동일)

- 방법1

runif(n=1,min=0,max=2)

[1] 0.7713044

- 방법2

runif(1)*2

[1] 1.197119

왜 방법2처럼 해도되는가? $X \sim U(0,1)$ 이면 $2X \sim U(0,2)$ 이므로
이론적인 내용은 수리통계학에서 배울것임

- 10000개의 값을 뽑아보고 평균이 1인지 살펴보자.

runif(n=10000,min=0,max=2) %>% mean

[1] 0.9979827

(runif(10000)*2) %>% mean

[1] 0.9942539

예제3: 베르누이분포 ($n=1$인 이항분포)

- 목표: 동전을 한번 던져서 앞면이 나오는 경우를 생성하고 싶음 (이때 동전은 공평한 동전이라 가정하자. 즉 앞면이 나올 횟수와 뒷면이 나올 횟수는 같다고 하자.)

rbinom(1,size=1,prob=0.5)

[1] 0

앞면이 1, 뒷면이 0이라고 생각

- 잘 되는것이 맞는지 한번 체크해보자. (10000개 정도의 동전을 던져서 결과를 저장하고 평균이 0.5인지 체크하자)

rbinom(10000,size=1,prob=0.5) %>% mean

[1] 0.5022

- 동전이 공평하지 않다고 하자. 예를들어 앞면이 나올 확률이 0.8이라고 생각하자.

rbinom(10000,size=1,prob=0.8) %>% mean

[1] 0.8083

예제4: 이항분포

- 목표: 동전을 10회던져 앞면이 나올 총 횟수를 생성하자.

동전을 10회던져 앞면이 나온 횟수를 $X$라고 하자. $X$를 생성하라.

rbinom(1,size=10,prob=0.5)

[1] 3

- 평균적으로 5가 나와야 한다. (10000번 정도 뽑아보고 관찰하자)

rbinom(10000,size=10,prob=0.5) %>% mean

[1] 5.0029

- 베르누이 분포는 이항분포의 특수한 경우이다.

- 직관적으로 생각하면 아래와 같이 이항분포를 생성해도 될 것 같다.

rbinom(10,size=1,prob=0.5) %>% sum

[1] 4

X=c() 
for (i in 1:10000) X[i] = rbinom(10,size=1,prob=0.5) %>% sum

X

   [1]  6  5  7  2  2  5  4  5  5  2  3  5  5  4  5  5  4  6  5  5  3  3  5  3
  [25]  5  3  9  4  8  3  3  5  5  3  6  4  7  4  7  6  6  4  5  5  7  3  5  2
  [49]  6  5  5  5  4  6  8  5  2  4  4  6  6  4  7  5  4  1  4  5  7  7  4  5
  [73]  4  5  4  7  4  4  3  4  6  6  5  4  6  5  5  5  6  4  5  6  4  4  8  6
  [97]  5  5  3  8  4  4  5  2  4  3  2  3  5  4  5  7  4  4  9  5  4  5  5  5
 [121]  7  4  6  6  6  5  3  5  7  6  7  4  3  6  4  3  4  6  5  6  6  5  6  9
 [145]  4  6 10  2  6  3  7  3  6  2  2  5  6  6  5  3  3  5  5  7  5  5  5  5
 [169]  6  3  5  6  3  6  3  6  4  6  5  5  4  5  6  3  5  6  3  5  6  6  4  7
 [193]  5  4  3  8  5  5  5  4  7  5  5  5  5  4  6  6  5  4  3 10  5  4  5  4
 [217]  4  5  4  9  4  5  3  4  4  3  5  3  4  4  5  7  3  4  6  5  8  5  4  3
 [241]  5  4  4  6  7  5  7  5  4  8  3  6  5  5  3  5  5  3  2  5  5  4  6  6
 [265]  6  4  6  5  3  3  6  6  5  5  4  5  5  4  7  6  5  4  3  4  4  4  5  5
 [289]  3  6  6  3  6  4  6  7  4  6  3  4  5  3  5  3  5  7  4  2  6  4  5  4
 [313]  5  7  8  1  7  4  5  8  5  3  6  6  6  7  5  8  4  6  6  6  4  6  3  6
 [337]  5  4  5  6  8  7  3  5  6  6  5  5  4  3  8  4  6  4  3  2  5  8  5  6
 [361]  2  4  5  5  5  5  6  5  4  5  4  4  2  6  5  7  3  3  6  6  9  6  4  6
 [385]  4  5  8  4  7  5  6  5  5  8  6  4  5  3  6  6  5  7  6  8  3  3  6  7
 [409]  8  5  5  5  6  6  5  4  6  3  4  4  5  5  6  6  5  5  8  4  5  6  2  4
 [433]  3  1  4  6  7  5  3  5  9  8  5  6  2  4  5  5  6  6  5  5  7  6  6  4
 [457]  7  4  4  7  7  2  3  7  3  7  4  5  6  4  4  4  5  7  2  4  8  7  6  6
 [481]  5  4  4  6  5  5  3  7  2  3  8  4  4  3  5  3  6  4  6  6  5  4  4  7
 [505]  3  6  3  5  2  3  4  5  6  4  3  4  7  5  7  6  4  6  7  2  6  8  5  2
 [529]  7  5  3  8  7  6  5  3  6  5  7  2  6  3  5  5  3  4  6  2  2  4  6  3
 [553]  5  6  8  4  4  7  5  6  2  8  7  5  6  5  1  7  7  2  5  4  5  7  4  6
 [577]  7  3  6  7  6  4  4  7  5  5  5  5  4  1  5  5  5  5  6  1  5  2  3  2
 [601]  3  5  6  0  5  4  3  5  6  2  5  3  5  6  6  6  2  4  4  4  3  5  7  5
 [625]  4  6  4  9  4  7  4  5  7  8  6  5  5  6  7  3  4  4  4  5  6  5  6  3
 [649]  6  7  5  5  3  4  5  3  4  6  3  4  3  7  3  5  6  6  5  6  6  3  3  5
 [673]  6  4  2  6  6  6  6  4  5  5  5  3  7  4  6  3  4  6  6  7  3  1  8  1
 [697]  7  6  8  7  4  7  3  3  3  4  3  8  4  4  6  3  6  2  2  5  3  5  8  6
 [721]  2  6  5  7  6  5  7  4  7  6  3  5  4  4  6  7  5  6  4  5  4  4  6  5
 [745]  4  7  6  5  6  7  5  7  2  5  2  5  3  6  5  7  7  6  4  3  2  1  9  6
 [769]  5  6  4  8  6  5  4  5  5  3  5  3  5  6  8  5  5  6  2  5  3  5  3  4
 [793]  4  5  5  4  4  6  4  4  4  5  5  5  6  7  5  3  7  6  5  7  7  5  6  7
 [817]  3  6  4  3  6  5  6  3  7  5  6  6  6  6  4  7  6  5  3  3  6  5  3  5
 [841]  5  6  5  4  6  6  5  6  5  4  5  6  3  4  5  5  4  8  3  6  4  5  6  8
 [865]  5  6  9  5  7  7  4  4  6  7  4  5  6  4  4  5  6  6  4  5  2  5  6  3
 [889]  5  5  5  5  3  5  5  7  3  6  6  7  5  2  5  5  8  4  3  4  6  6  3  5
 [913]  7  7  5  6  5  6  4  3  7  4  9  9  7  4  3  6  4  3  4  3  5  6  5  4
 [937]  5  5  4  4  2  3  4  7  6  6  4  5  5  7  5  8  6  8  5  5  6  6  4  6
 [961]  6  2  5  8  8  5  7  6  6  7  7  8  3  5  3  4  6  6  6  5  7  3  6  4
 [985]  1  5  5  2  4  7  4  8  5  6  3  5  4  6  6  5
 [ reached getOption("max.print") -- omitted 9000 entries ]

직관적으로 생각하면 X는 rbinom(10000,size=10,prob=0.5) 와 같은 효과 (랜덤성 때문에 같은 값은 아님) 일 것 같은데 이는 실제로 그러하다. (좀 더 엄밀한 증명은 수리통계 시간에 배운다)

X %>% mean()

[1] 5.0216

어떠한 집합에서 샘플링

예제5: 비복원추출

- 주머니에서 빨간공 2개와 파란공 3개가 있다고 하자.

set = c(rep("red",2),rep("blue",3))
set

[1] "red"  "red"  "blue" "blue" "blue"

- 이러한 집합에서 랜덤으로 3개의 원소를 고르고 싶다면? (단, 뽑은 공은 다시 넣지 않는다)

sample(set,size=3)

[1] "red"  "blue" "red"

- 50번정도 해보자.

for (i in 1:50) sample(set,size=3) %>% print

[1] "red"  "blue" "blue"
[1] "blue" "red"  "red" 
[1] "blue" "blue" "red" 
[1] "blue" "blue" "blue"
[1] "blue" "blue" "red" 
[1] "red"  "blue" "blue"
[1] "red"  "red"  "blue"
[1] "red"  "blue" "blue"
[1] "red"  "red"  "blue"
[1] "blue" "red"  "blue"
[1] "blue" "blue" "red" 
[1] "red"  "red"  "blue"
[1] "red"  "blue" "blue"
[1] "blue" "red"  "blue"
[1] "blue" "blue" "red" 
[1] "blue" "blue" "red" 
[1] "red"  "red"  "blue"
[1] "blue" "red"  "blue"
[1] "red"  "blue" "blue"
[1] "blue" "blue" "red" 
[1] "red"  "blue" "blue"
[1] "red"  "blue" "blue"
[1] "blue" "blue" "blue"
[1] "blue" "red"  "blue"
[1] "blue" "red"  "red" 
[1] "blue" "blue" "red" 
[1] "red"  "blue" "blue"
[1] "blue" "red"  "red" 
[1] "blue" "red"  "blue"
[1] "blue" "blue" "blue"
[1] "blue" "blue" "blue"
[1] "blue" "blue" "red" 
[1] "red"  "blue" "blue"
[1] "blue" "red"  "red" 
[1] "red"  "blue" "blue"
[1] "blue" "red"  "blue"
[1] "blue" "blue" "red" 
[1] "red"  "blue" "red" 
[1] "blue" "red"  "blue"
[1] "red"  "blue" "blue"
[1] "blue" "red"  "blue"
[1] "blue" "red"  "blue"
[1] "blue" "blue" "blue"
[1] "red"  "red"  "blue"
[1] "blue" "blue" "blue"
[1] "blue" "blue" "red" 
[1] "blue" "blue" "blue"
[1] "blue" "red"  "red" 
[1] "blue" "blue" "blue"
[1] "blue" "red"  "red"

빨간공을 3번 뽑지는 못함

예제6: 복원추출

- 뽑은공을 다시 넣는다고 가정한다면?

set = c(rep("red",2),rep("blue",3))
set

[1] "red"  "red"  "blue" "blue" "blue"

for (i in 1:50) sample(set, size=3,replace = TRUE) %>% print

[1] "blue" "red"  "blue"
[1] "red"  "blue" "red" 
[1] "blue" "blue" "red" 
[1] "red"  "blue" "blue"
[1] "blue" "blue" "red" 
[1] "red"  "red"  "blue"
[1] "red"  "blue" "blue"
[1] "red"  "blue" "blue"
[1] "blue" "blue" "red" 
[1] "red" "red" "red"
[1] "blue" "red"  "red" 
[1] "blue" "red"  "blue"
[1] "blue" "blue" "red" 
[1] "blue" "blue" "blue"
[1] "blue" "red"  "blue"
[1] "blue" "red"  "blue"
[1] "blue" "blue" "red" 
[1] "red"  "red"  "blue"
[1] "red"  "blue" "red" 
[1] "red"  "blue" "red" 
[1] "red"  "blue" "red" 
[1] "blue" "red"  "red" 
[1] "blue" "red"  "blue"
[1] "blue" "blue" "blue"
[1] "red" "red" "red"
[1] "red"  "red"  "blue"
[1] "blue" "red"  "blue"
[1] "red" "red" "red"
[1] "blue" "blue" "blue"
[1] "blue" "red"  "red" 
[1] "blue" "blue" "red" 
[1] "blue" "red"  "red" 
[1] "red"  "blue" "red" 
[1] "blue" "red"  "blue"
[1] "red"  "blue" "blue"
[1] "blue" "red"  "blue"
[1] "blue" "red"  "blue"
[1] "blue" "blue" "blue"
[1] "red"  "blue" "blue"
[1] "red"  "red"  "blue"
[1] "blue" "red"  "red" 
[1] "blue" "blue" "blue"
[1] "blue" "red"  "red" 
[1] "red"  "blue" "blue"
[1] "blue" "blue" "blue"
[1] "blue" "red"  "blue"
[1] "red"  "red"  "blue"
[1] "blue" "blue" "red" 
[1] "red"  "red"  "blue"
[1] "blue" "blue" "blue"

빨간공을 3번 뽑을 수 있음

(인덱스,확률) 정보가 주어졌을 경우 샘플링

예제7: 예제6의 다른구현

- 예제6의 코드는 아래와 같다.

set = c('red','blue') ## 인덱스.. 
p = c(0.4,0.6) ## 인덱스에 해당하는 확률

sample(set,size=3, prob=p, replace = T)

[1] "red"  "blue" "blue"

- 아래와 같은 코드로도 응용가능하다.

sample(set,size=20,prob=c(0.95,0.05),replace=T)

 [1] "blue" "blue" "red"  "red"  "red"  "red"  "red"  "red"  "red"  "red" 
[11] "red"  "red"  "red"  "red"  "red"  "red"  "red"  "red"  "red"  "red"

전체집합에서 서로 다른 확률로 원소가 뽑히는 경우

예제8

- 예제7의 코드에서 replace를 제외한다면?

set = c('red','blue') ## 인덱스..

sample(set,size=20,prob=c(0.95,0.05))

Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'
Traceback:

1. sample(set, size = 20, prob = c(0.95, 0.05))
2. sample.int(length(x), size, replace, prob)

당연히 에러가남 (원소가 2개 있는 set에서 20개를 뽑으라고 하니..)

- 아래는 에러가 나지 않을 것임

sample(set,size=2,prob=c(0.95,0.05))

[1] "red"  "blue"

- 의미는?

rslt = matrix(rep("-",1000*2),ncol=2)
head(rslt)

     [,1] [,2]
[1,] -    -   
[2,] -    -   
[3,] -    -   
[4,] -    -   
[5,] -    -   
[6,] -    -

for (i in 1:1000) rslt[i,] = sample(set,size=2,prob=c(0.95,0.05))

head(rslt)

     [,1] [,2]
[1,] red  blue
[2,] red  blue
[3,] red  blue
[4,] red  blue
[5,] red  blue
[6,] red  blue

sum(rslt[,1]=='red')/1000

[1] 0.947

대략 95%

- 첫번째가 붉은 공이 뽑힐 확률은 95%임. 그리고 두번째 공은 첫번째 공과 다른색이 뽑힘

숙제: 징검다리 건너기

(유리,강화유리)의 쌍으로 이루어진 징검다리가 총 5개 있다고 하자. (따라서 징검다리는 모두 10개이다)

강화유리로 된 징검다리를 밟으면 살아남지만 유리로 된 징검다리를 밟으면 죽는다.

따라서 강화유리로 된 징검다리를 계속 골라야 살아남을 수 있다.

A씨는 유리공장에서 20년 근무한 장인으로 유리와 강화유리를 구분할수 있는 능력을 가졌다고 하자.

그래서 강화유리로 된 징검다리를 고를 확률이 80%라고 하자.

A씨가 살아남을 확률을 시뮬레이션을 통하여 구하여라.

Hint: 총 1000번의 시뮬레이션을 수행하고 그중에서 A씨가 살아남는 케이스가 몇회정도 되는지 찾아라.

sample = rbinom(1000,size=5,prob=0.8)
sample

X <- sample == 5
X

sum(X)/1000*100