4.8 산불 문제
이것은 Conway의 life game이 단순화된 형태이다. (https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life)
행x열이 41 x 41인 나무가 심어져 있는 산(숲)이 있다고 하자.
중심의 한 나무에 불이 붙었고, 어느 나무에 불이 붙었을 때 상하 또는 좌우의 어느 한 나무에 불이 옮겨 붙을 확률이 p라면 평균적으로 얼마 동안 불이 나고, 몇 그루의 나무가 타겠는가?
불이 붙은 나무는 다음 단위시간이 지난 후에 나무가 다 타서 불이 꺼지며, 이미 탄 나무에는 다시 불이 붙지 않는다. 불은 대각선으로도 옮겨 붙지 않는다고 가정하자.
한 번 simulation하기 위한 함수는 다음과 같다.
= function(p=0.5, nRow=41, nCol=41, PLOT = T)
Fire
{= matrix(0, nrow=nRow, ncol=nCol)
mTree1 median(1:nRow), median(1:nCol)] = 1 # 0=Alive, 1=onFire, 2=Burnt
mTree1[= (sum(mTree1 == 1) > 0)
fAlive = 0
fireTime = mTree1
mTree2 while (fAlive) {
for (j in 1:nRow) {
for (k in 1:nCol) {
if (mTree1[j, k] == 1) {
= 2
mTree2[j, k] if (j < nRow) if (mTree1[j + 1, k] == 0) mTree2[j + 1, k] = rbinom(1, 1, p)
if (j > 1) if (mTree1[j - 1, k] == 0) mTree2[j - 1, k] = rbinom(1, 1, p)
if (k < nCol) if (mTree1[j, k + 1] == 0) mTree2[j, k + 1] = rbinom(1, 1, p)
if (k > 1) if (mTree1[j, k - 1] == 0) mTree2[j, k - 1] = rbinom(1, 1, p)
}
}
}= mTree2
mTree1 = fireTime + 1
fireTime = (sum(mTree1 == 1) > 0)
fAlive if (PLOT) {
filled.contour(mTree1, nlevels=3, col=c("green", "red", "black", "black"), axes=F)
Sys.sleep(0.1)
}
}= c(burntTree = sum(mTree1 == 2), fireTime = fireTime)
Res return(Res)
}
다음과 같이 하면 1회 simulation 결과를 볼 수 있다.
Fire()
만약 1000회 simulation한다면 다음과 같이 할 수 있다.
= 1000
N = data.frame(nTree = rep(NA, N), Time = rep(NA, N))
Res = Sys.time()
t1 for (i in 1:N) Res[i, ] = Fire(PLOT=F)
Sys.time() - t1
Time difference of 2.218893 secs
hist(Res[, "Time"], breaks=1:max(Res[,"Time"]), xlab="Duratino of fire", main="")
Figure 4.10: Histogram of fire duration
#dev.new()
hist(Res[, "nTree"], breaks=1:max(Res[,"nTree"]), xlab="Count of burnt tree", main="")
Figure 4.11: Histogram of burnt trees
#dev.new()
plot(nTree ~ Time, Res, xlab="Duration of fire", ylab="Count of burnt tree")
Figure 4.12: Burnt trees vs. duration of fire