drunkard’s walk

---

【模擬 1：模擬醉漢走路過程(單次，紀錄位置變化與步數)】

POSITION <- 0 ; TIMES <- 0 
repeat{POSITION <- POSITION + sample(c(-1,1),1);print(POSITION) ; TIMES <- TIMES + 1; if(POSITION==-1 || POSITION==2) {break}}

## [1] 1
## [1] 2

print(TIMES)

## [1] 2

【模擬 2：模擬醉漢走路100,000次，共十回，估算醉漢撞到左邊牆壁的機率】

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:100000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(-1,1),1) ; TIMES <- TIMES + 1; if(POSITION==-1 || POSITION==2) {break}}
    
    ifelse(POSITION==-1, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/100000)
}

## [1] 0.66687
## [1] 0.66725
## [1] 0.66929
## [1] 0.66841
## [1] 0.6671
## [1] 0.66648
## [1] 0.66581
## [1] 0.66857
## [1] 0.66603
## [1] 0.66356

【計算理論值(利用古典破產理論)】

u0 <- (3-1)/3
u0

## [1] 0.6666667

【模擬 3：模擬醉漢走路100,000次(亦可修改參數i，讓模擬次數增加到1,000,000次或更多)，共十回，估算在左右方三步之處有水溝，且醉漢離家100步的情況下，醉漢最終能安然返家的機率】

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:100000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(-1,1),1) ; TIMES <- TIMES + 1; if(POSITION==-3 || POSITION==3 || TIMES==100) {break}}
    
    ifelse(TIMES==100, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/100000)
}

## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0

【計算理論值(利用一步轉移機率矩陣的100次方，即100步轉移機率矩陣來計算)】

P <- matrix(c(1,0,0,0,0,0,0,.5,0,.5,0,0,0,0,0,.5,0,.5,0,0,0,0,0,.5,0,.5,0,0,0,0,0,.5,0,.5,0,0,0,0,0,.5,0,.5,0,0,0,0,0,0,1),nrow=7)
P10 <- P%*%P%*%P%*%P%*%P%*%P%*%P%*%P%*%P%*%P
P100<- P10%*%P10%*%P10%*%P10%*%P10%*%P10%*%P10%*%P10%*%P10%*%P10
sum(P100[2:6,4])

## [1] 7.550955e-07

【計算理論值(利用小岢猜想)】

p100 <- 0
for (i in 0:floor((100-3)/2)) { p100 <- p100 + 3^(i)*2/2^(2*i+3)}
1-p100

## [1] 7.550955e-07

【模擬 4：模擬醉漢走路100,000次，共十回，估算在左右方三步之處有水溝，且醉漢離家僅僅四步的情況下，醉漢最終能安然返家的機率】

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:100000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(-1,1),1) ; TIMES <- TIMES + 1; if(POSITION==-3 || POSITION==3 || TIMES==4) {break}}
    
    ifelse(TIMES==4, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/100000)
}

## [1] 0.74888
## [1] 0.75131
## [1] 0.74995
## [1] 0.75045
## [1] 0.75152
## [1] 0.75064
## [1] 0.75019
## [1] 0.75036
## [1] 0.74975
## [1] 0.75303

【模擬 5：模擬醉漢走路100,000次，共十回，分別估算在左右方三步之處有水溝，且醉漢離家為10步與50步的情況下，醉漢最終能安然返家的機率】

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:100000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(-1,1),1) ; TIMES <- TIMES + 1; if(POSITION==-3 || POSITION==3 || TIMES==10) {break}}
    
    ifelse(TIMES==10, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/100000)
}

## [1] 0.3179
## [1] 0.31627
## [1] 0.31586
## [1] 0.31819
## [1] 0.31672
## [1] 0.31419
## [1] 0.31591
## [1] 0.31757
## [1] 0.31707
## [1] 0.3176

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:100000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(-1,1),1) ; TIMES <- TIMES + 1; if(POSITION==-3 || POSITION==3 || TIMES==50) {break}}
    
    ifelse(TIMES==50, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/100000)
}

## [1] 0.00088
## [1] 0.00094
## [1] 0.00113
## [1] 0.00092
## [1] 0.00082
## [1] 0.00101
## [1] 0.00097
## [1] 0.00096
## [1] 0.00085
## [1] 0.00094

【計算理論值(利用一步轉移機率矩陣的10次方、50次方，即10步、50步轉移機率矩陣來計算)】

P <- matrix(c(1,0,0,0,0,0,0,.5,0,.5,0,0,0,0,0,.5,0,.5,0,0,0,0,0,.5,0,.5,0,0,0,0,0,.5,0,.5,0,0,0,0,0,.5,0,.5,0,0,0,0,0,0,1),nrow=7)
P10 <- P%*%P%*%P%*%P%*%P%*%P%*%P%*%P%*%P%*%P
sum(P10[2:6,4])

## [1] 0.3164062

P50<- P10%*%P10%*%P10%*%P10%*%P10
sum(P50[2:6,4])

## [1] 0.001003391

【計算理論值(利用小岢猜想)】

p10 <- 0
for (i in 0:floor((10-3)/2)) { p10 <- p10 + 3^(i)*2/2^(2*i+3)}
1-p10

## [1] 0.3164062

p50 <- 0
for (i in 0:floor((50-3)/2)) { p50 <- p50 + 3^(i)*2/2^(2*i+3)}
1-p50

## [1] 0.001003391

【模擬 6：模擬爸爸和小晴對決100,000次，共十回，估算在小晴口袋裡有3元，且爸爸口袋裡有10元的情況下，小晴最終能讓爸爸破產的機率】

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:100000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(-1,1),1) ; TIMES <- TIMES + 1; if(POSITION==-3 || POSITION==10) {break}}
    
    ifelse(POSITION==10, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/100000)
}

## [1] 0.22979
## [1] 0.22969
## [1] 0.23088
## [1] 0.23079
## [1] 0.23133
## [1] 0.23131
## [1] 0.22875
## [1] 0.23157
## [1] 0.23053
## [1] 0.23112

【計算理論值(利用古典破產理論)】

u3 <- 3/(10+3)
u3

## [1] 0.2307692

【模擬 7：模擬爸爸和賭場對決100次(或想成100位好賭的爸爸和賭場對決)，共十回，估算在爸爸的口袋裡有10元，且賭場的資產是爸爸口袋內錢數的10,000倍(100,000)的情況下，爸爸最終能贏走賭場所有資產的機率】

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:100) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(-1,1),1) ; TIMES <- TIMES + 1; if(POSITION==-10 || POSITION==100000) {break}}
    
    ifelse(POSITION==100000, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/100)
}

## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0

【計算理論值(利用古典破產理論)】

u10 <- 10/(100000+10)
u10

## [1] 9.999e-05

【模擬 8：模擬爸爸和前述的賭場對決10次(或想成10位好賭的爸爸和賭場對決)，觀察在經歷過多少次賭局後，好賭的爸爸才輸光了口袋裡的10元】

for (i in 1:10) {
  POSITION <- 0 ; TIMES <- 0 
  repeat{POSITION <- POSITION + sample(c(-1,1),1) ; TIMES <- TIMES + 1; if(POSITION==-10 || POSITION==100000) {print(TIMES) ; break}}
  ifelse(POSITION==100000, INDEX <- INDEX + 1, INDEX <- INDEX)
}

## [1] 506
## [1] 15128
## [1] 924
## [1] 190
## [1] 292
## [1] 26
## [1] 3408
## [1] 192
## [1] 3479820
## [1] 214

print(INDEX)

## [1] 0

【計算理論值(利用古典破產理論)】

D10 <- 10*(100010-10)
D10

## [1] 1e+06

【模擬 9：模擬爸爸和前述的賭場對決4次(或想成4位好賭的爸爸和賭場對決)，透過繪製折線圖，觀察好賭的爸爸們如何輸光了口袋裡10元的詳細過程】

par(mfrow=c(2,2))
for (i in 1:4){
  POSITION <- 0 ; TIMES <- 0 ; x <- c() 
  repeat{POSITION <- POSITION + sample(c(-1,1),1) ; x <- rbind(x, POSITION); TIMES <- TIMES + 1; if(POSITION==-10 || POSITION==100000) {print(TIMES) ; break}}
  
  plot(1:length(x),x,xlab="次數",ylab="金額",type="l")
}    

## [1] 2904

## [1] 414

## [1] 4718

## [1] 140

【模擬 10：模擬爸爸和莊家對決10,000次(或想成10,000位好賭的爸爸和莊家對決)，共十回，估算在爸爸的口袋裡有10元，莊家的資產有100元，但爸爸的勝率分別只有499/1000、49/100、2/5的情況下，爸爸最終能讓莊家破產的機率】

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:10000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(rep(1,499),rep(-1,501)),1) ; TIMES <- TIMES + 1; if(POSITION==-10 || POSITION==100) {break}}
    
    ifelse(POSITION==100, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/10000)
}

## [1] 0.0774
## [1] 0.0731
## [1] 0.0741
## [1] 0.0709
## [1] 0.0724
## [1] 0.0731
## [1] 0.0679
## [1] 0.0772
## [1] 0.0735
## [1] 0.0718

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:10000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(rep(1,49),rep(-1,51)),1) ; TIMES <- TIMES + 1; if(POSITION==-10 || POSITION==100) {break}}
    
    ifelse(POSITION==100, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/10000)
}

## [1] 0.006
## [1] 0.0059
## [1] 0.0063
## [1] 0.0057
## [1] 0.0046
## [1] 0.0062
## [1] 0.0063
## [1] 0.0073
## [1] 0.0054
## [1] 0.0065

for (j in 1:10) {
  INDEX <- 0
  for (i in 1:10000) {
    POSITION <- 0 ; TIMES <- 0 
    repeat{POSITION <- POSITION + sample(c(-1,1,-1,1,-1),1) ; TIMES <- TIMES + 1; if(POSITION==-10 || POSITION==100) {break}}
    
    ifelse(POSITION==100, INDEX <- INDEX + 1, INDEX <- INDEX)
  }
  print(INDEX/10000)
}

## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0
## [1] 0

【計算理論值(利用古典破產理論)】

ua <- 1-((501/499)^110-(501/499)^10)/((501/499)^110-1)
ua

## [1] 0.07383794

ub <- 1-((51/49)^110-(51/49)^10)/((51/49)^110-1)
ub

## [1] 0.006110713

uc <- 1-((3/2)^110-(3/2)^10)/((3/2)^110-1)
uc

## [1] 0

---