單元摘要： 程式、機率與決策


pacman::p_load(dplyr, vcd)

【A】實務案例一：老忠實賭局

🗿 問題討論：
假如賭場老闆從零開始把每10秒鐘設為一個區間，每個區間的賭金是五塊錢 …
 ■ 你要怎麼押注，才能獲得最高的期望值呢？
 ■ 你的賭金和期望值各是多少？
將獲利的期望值和賭金的比值稱為「期望投資報酬率」 …
 ■ 「最大期望投資報酬率」和「最大期望報酬」的投資策略是一樣的嗎？
 ■ 你認為哪一個策略目標比較好呢？ 為什麼？

D = faithful$eruptions        # copy data to a short name
Adjust = 0.5                       # bandwidth adjustment
DEN = density(D, adjust = Adjust)      # densitu function
PDF = approxfun(DEN$x, DEN$y, yleft=0, yright=0)    # PDF
x = seq(1, 6, 1/6)
px = sapply(2:length(x), function(i) integrate(PDF, x[i-1], x[i])$value)
payoff = 100 * px - 5
df = data.frame(x2=x[-1], px, payoff)  # x is longer than px by 1 
df = df[order(-df$payoff),]            # order by decreasing payoff
df$cumsum =  cumsum(df$payoff)         # cumulative sum
round(df[df$payoff > 0,], 3)           # round to the 3rd place
      x2    px payoff cumsum
21 4.500 0.096  4.577  4.577
22 4.667 0.089  3.938  8.515
20 4.333 0.089  3.932 12.447
6  2.000 0.082  3.226 15.673
19 4.167 0.077  2.654 18.327
7  2.167 0.070  1.968 20.294
23 4.833 0.068  1.843 22.137
5  1.833 0.064  1.367 23.504
18 4.000 0.058  0.834 24.338
🗿： 至於哪一個策略目標比較好，標準答案是“It Depends!” Therefore, the real question is “What does it depends?”

【B】實務案例二：天后宮的擲筊

🗿 問題討論：
在cup.csv檔案裡面是天后宮一千個信徒的擲筊次數，假定每一個人都是擲到3次成功才停止，請問：
  ■ 這個筊的成功機率大約是？
  ■ 請畫出用這個筊擲10次成功之前，失敗次數的機率分布
  ■ 用這個筊擲15次還不能有5次成功的機率是？

失敗次數是Poisson嗎？

X = read.csv("cup.csv")$x
fit =  goodfit(table(X - 3), "poisson")
summary(fit)

     Goodness-of-fit test for poisson distribution

                    X^2 df    P(> X^2)
Likelihood Ratio 1605.4 33 5.0416e-317

已知成功次數為3，用Negative Binomial試試看

fit =  goodfit(table(X - 3), "nbinom", par=list(size=3))
summary(fit)

     Goodness-of-fit test for nbinomial distribution

                    X^2 df P(> X^2)
Likelihood Ratio 37.564 33  0.26795
fit$par
$size
[1] 3

$prob
[1] 0.248

這個筊的成功機率大約是？ 0.248
畫出用這個筊擲10次成功之前，失敗次數的機率分布

par(cex=0.7)
dnbinom(0:80, 10, fit$par$prob) %>% 
  barplot(names=0:80, xlab="No. Failures before 10 Successes")

用這個筊擲15次還不能有5次成功的機率是？ 0.693

1 - pnbinom(10, 5, fit$par$prob)
[1] 0.69309
par(cex=0.5)
dnbinom(0:40, 5, fit$par$prob) %>% 
  barplot(names=0:40, xlab="No. Failures before 10 Successes",
          col=c(rep('green',11), rep('red', 30)))

💡 學習重點：
  ■ 即使是很簡單的一個案例，我們也需要綜合使用機率、程式和商業知識，才能做出合適的決策！
  ■ 機率分布的商業應用套路：
    1. 寫出各事件(變數區間)的報償
    2. 求出變數的機率分布函數
    3. 求出各事件的機率
    4. 求出不同策略的期望報償
    5. 依公司的策略目標做選擇
  ■ 如何求出機率分布函數？
    A. 將實證資料fit進已知的理論分布、並求出參數
    B. 直接就實證資料做成(平滑)機率密度函數