【무한대의 기대치? 】 상트 페테르부르크의 역설을 R로 시뮬레이션

요비노리 타쿠미씨의 동영상 「 기대치가 무한대인 베팅(상트페테르부르크의 패러독스) 」를 보고 재미있을 것 같았으므로 시뮬레이션해 보았습니다.

문제 설정

문제의 게임은
확률 1/2의 코인을 표가 나올 때까지 던지고, k회째에 표가 나오면 2^(k-1)엔상금을 받을 수 있다고 하는 게임입니다.

예를 들어, 뒷뒤 표, 즉 네 번째로 표가 나오면
2^(4-1) = 8(엔)
됩니다.

이 게임, 절반의 확률로 1엔, 1000엔 이상 받을 확률은 0.1%에도 못 미친다.
한 번 10만엔으로 할 수 있다고 해도 절대로 하고 싶지 않은 베팅입니다만,
기대치를 계산하면 무한대이며 분명 직감에 반하고 있습니다.

이 문제를 해결하기 위해 다양한 방법이 생각되었지만,
윌리엄 펠러는 "표본 추출"이라는 통계학의 개념을 사용했습니다.
기대치는 「모집단의 평균치」이며, 실제로 획득하는 상금액은 「표본 평균」이라고 생각하는 것으로,
참여 횟수 n에 따라 표본 평균이 변경되었음을 나타냅니다.

시뮬레이션

이번 시뮬레이션하는 것은, 처음 동영상 의 20:30 정도부터 그리고 있는 그래프입니다.

n회 참가할 때의 1회당의 획득 상금은 n을 크게 하면 log2(n)에 가까워진다(정확하게는 확률 수렴)인 것 같기 때문에 해 봅니다.

# サンクトペテルブルクのパラドックス

# 賞金を出力する関数を作る
reward <- function(){

      k <- 0 # 試行回数のカウント
    result <- 0 # コイン投げの結果：0を裏、1を表とする

        ## 成功するまで処理を繰り返す 
        while (result == 0) {
            k <- k + 1

            ##　コイン投げ
            result <- sample(0:1, size = 1, replace = TRUE, prob = c(0.5, 0.5))
    }
    return(2^(k-1)) # 賞金
}

# 賭けをやってみよう
cat("賞金は", reward(), "円")
# >>賞金は 2 円

# 賭けにn回参加した時の一回当たりの利得(S(n)/n)を計算する関数
rew_per_game <- function(n){
    S_n <- 0 # 合計利得
    for (i in 1:n){
        S_n <- S_n + reward()
    }
    return(S_n/n)
}

# グラフ描画

plotGraph <- function(m, n, ylim, by){

    par(new =F) # 新しいグラフを描画

    plot(log2(0:n+1), type = "l", col = "red", xlim = c(0,n+1), ylim = c(0,ylim),
         xlab = '', ylab = '', pch = 20) # 理論値log2(n)のグラフ
    par(new = TRUE) # 次に描画するグラフを重ねる
    part = seq(1,n+1,by=by) # シミュレーションする参加回数のベクトル
    df <- data.frame(part = part, rew_per_game = 0, mean = 0) # シミュレーション結果を入れるデータフレーム
    for(i in 1:m){

        df <- data.frame(part = part, rew_per_game = 0, mean = df$mean) # データフレームの初期化

        df$rew_per_game <- lapply(df$part ,rew_per_game) # 各参加回数に対してシミュレーション

        plot(df$part, df$rew_per_game, xlim = c(0,n+1), ylim = c(0,ylim), pch = 20, cex = 0.01,
            xlab = '', ylab = '', axes = FALSE) # 結果のプロット
        par(new = TRUE)

        df$mean <- unlist(df$mean) + unlist(df$rew_per_game)/m # 平均値の計算
    }
    plot(df$part, df$mean, type = "b", col = "blue", xlim = c(0,n+1), ylim = c(0,ylim),
         xlab = '', ylab = '') # 平均値のグラフ
}

# パラメータの設定
m = 10 # サンプル数
n = 1000 # 参加回数をn回まで描画
ylim = 30 # y軸最大値の設定
by = 100 # 参加回数何回ごとにシミュレーションするか

plotGraph(m, n, ylim, by)

m, n, ylim, by 의 파라미터를 놀면서 여러가지 놀 수 있습니다.


↑이것은 m=500, n=10000, ylim=50, by=1000


↑시간이 걸리지만 이런 것도 할 수 있습니다(m=10, n=100000, ylim=30, by = 100)

실제로 시뮬레이션하면 변동이 너무 크고 전혀 결과가 안정되지 않는 것을 실감할 수 있습니다.
그런데 이것, m을 늘리면 1회당의 획득 상금도 늘어 버리는 것처럼 생각됩니다만 어떻습니까.