標準偏差の大きさは夢の大きさ―宝くじをJuliaで考える

手のひらに消えゆく白く冷ややかな刹那に、私は今年も儚い夢を重ねて見る。夢――、年末ジャンボ宝くじである。

問題提起

本記事では、
「連番10枚 / バラ10枚 / バラ単券10枚」の3つの買い方を確率論とJuliaで比較し、
一瞬のうちにただの紙切れへと堕ちる宝くじを、せめてもの最大限楽しむ方法について考察する。

---

宝くじの買い方には大きく、

• 連番
• バラ

の2種類があり、また、そこそこ当たる当たりとして、

• 下1桁一致：300円
• 下2桁一致：3000円

の2種がある。

大きな夢とは裏腹に、さしあたりこの2つの当選でいかに期待値の罠を裏切って元金割れを軽減するか、が私にとっての実際的な楽しみとなり果てる。

その点で、

• 連番10枚組は1枚みただけで10枚すべての下2桁が確定し、3000円の当たりの有無が分かってしまう
• バラでも連番でも10枚組で買うと300円の当たりは決まって1つ（10枚組の下1桁には必ず0〜9が1枚ずつ入るとされている -> 宝くじホームページ）

ため、元金回収具合のドキドキを最大限味わえていない。

結論

「宝くじを買うときは、バラ券を10枚組ではなく1枚ずつ別々に買うことで、当たりの不確実性が最大になり、くじとしての楽しみを最大限謳歌できる」

項目	連番10枚組	バラ10枚組	バラ単券10枚
期待値	600.0	600.0	600.0
標準偏差	900.0	900.0	985.9
3000円以上の当たりが出る確率	10.0%	10.0%	9.6%

(小数点以下第2位を四捨五入)

【図】連番10枚組、バラ10枚組、バラ単枚10枚を購入したときの、合計当たり額の確率分布

ちなみに最大合計当選額の33000円までのグラフ(logスケール)

33000円となる確率は、10の30乗分の1である。

---

• 連番10枚組/バラ10枚組/バラ単券10枚はいずれも期待値は600円である。
• 連番10枚組とバラ10枚組は、合計当たり額の確率分布が完全に同じである。
• これだけ不均一な分布を標準偏差で比較すべきかは微妙だが、バラ単券10枚の標準偏差が最も大きく、そのためバラ単券10枚の方が合計当たり額の不確実性が最大となる。
• (意外にも、)3000円以上の当たりが出る確率(元金を回収する確率)はバラ単券10枚が最も低い。

期待値は同じ600円でも、300円か3300円の2択しか出ないくじと、色んな可能性があるくじなら、後者の方がきっと楽しい。

問題設定

• 当たりは次の2種のみ考慮し、当たり番号は独立に決まる。（重複当選もあり得る、公式ルール）
  • 下1桁一致：300円
  • 下2桁一致：3000円
• バラ券を10枚組で買うと、組の中の1の位は0~9が1枚ずつ入る。
• バラ券を10枚別々に買うと、くじ券の数字は完全にランダムになる。

連番10枚組/バラ10枚組/バラ単券10枚の3つのパターンで10枚のくじを購入したときの、合計当たり額の分布を考える。

考察

バラ券単枚

モデル化

くじが各回それぞれ独立に以下の報酬が決まると読み替える。

• 報酬 $R_1=300$ 円が確率 $p_1=1/10$ で当たる。
• 報酬 $R_2=3000$ 円が確率 $p_2=1/100$ で当たる

本来は、例えば下1桁が1のくじの総数は限られているので、購入して手元に来るくじ番号は独立ではなく、厳密には独立に報酬は決まらないはずである。
しかし実際のところ、くじ番号が全6桁であることや組番号やユニットの違いも踏まえれば、同じ下1桁のくじの総数は購入するくじの枚数に比べ大量にあるので、購入して手元に来るくじ番号の下1,2桁は独立とみなしてよいだろう。

また、本来は当たり番号がすべてのくじで共通して当たりごとに1つ決まっていて、それとくじ番号を照らし合わせるわけだが、上記の独立性を踏まえて、各々のくじが各々の当たりについて独立に報酬が決まるものとしてモデル化する。

上記より、

$N$ 回くじを引くとき、

• $Y_1$：$R_1$ の当たり回数
• $Y_2$：$R_2$ の当たり回数

とすると、

$$
Y_1 \sim \mathrm{Binomial}(N, p_1),\qquad
Y_2 \sim \mathrm{Binomial}(N, p_2),
$$

と独立な二項分布として表し、合計報酬額 $S_N$ は

$$
S_N = R_1 Y_1 + R_2 Y_2.
$$

と表す。この $S_N$ の確率分布を知りたいということになる。

確率分布の導出・計算

目標は、 $S_N$ が具体的な値 $s$ になる確率 $P(S_N = s)$ を全ての $s$ について知ることである。

$Y_1, Y_2$ の取りうる値は、

$$
y_1, y_2 \in {0,1,2,\dots,N}
$$

したがって $S_N$ の取りうる値は、

$$
s = R_1 y_1 + R_2 y_2 \quad \left( y_1, y_2 \in {0,1,2,\dots,N} \right)
$$

となる。

ただし、異なる $(y_1, y_2)$ の組み合わせで同じ $s$ となりうることに注意して、 $P(S_N = s)$ は、

$$
P(S_N = s)
= \sum_{\substack{s = R_1 y_1 + R_2 y_2 \ y_1, y_2 \in {0,1,2,\dots,N} }}
P( Y_1 = y_1, Y_2 = y_2 ).
$$

と表せる。

たとえば、 $s=3000$ を与える $(y_1, y_2)$ は $(10,0), (0,1)$ の2通りあるので、
$$
P(S_N = 3000)
= P( Y_1=10, Y_2=0) + P( Y_1=0, Y_2=1)
$$
となる。

さらに、 $Y_1, Y_2$ は独立なので、

$$
P(Y_1=y_1, Y_2=y_2)
= P(Y_1=y_1) \times P(Y_2=y_2).
$$

また、二項分布の PMF (probability mass function) は

$$
\begin{aligned}
P(Y_1=y_1) &= \binom{N}{y_1} p_1^{y_1} (1-p_1)^{N-y_1} \\
P(Y_2=y_2) &= \binom{N}{y_2} p_2^{y_2} (1-p_2)^{N-y_2}
\end{aligned}
$$

より、

$$
P(S_N = s)
= \sum_{\substack{s = R_1 y_1 + R_2 y_2 \ y_1, y_2 \in {0,1,2,\dots,N} }}
\binom{N}{y_1} p_1^{y_1} (1-p_1)^{N-y_1}
\binom{N}{y_2} p_2^{y_2} (1-p_2)^{N-y_2}
$$

と表せる。

これを解析的に評価して期待値や分散を計算するのは難しそうなので、今回は実際に上式を計算により展開して、

$$
\begin{align*}
E[S_N] &= \sum_i s_i p_i \\
\mathrm{Var}(S_N) &= \sum_i (s_i - E[S_N])^2 p_i \\
\mathrm{SD}(S_N) &= \sqrt{\mathrm{Var}(S_N)}
\end{align*}
$$

期待値 $E$ 、 分散 $\mathrm{Var}$ 、 標準偏差 $\mathrm{SD}$ をそれぞれ求める。

Juliaスクリプト

using DataFrames
using Distributions
using Statistics
using Random
using Plots

plotlyjs()

if Sys.iswindows()
  Plots.default(fontfamily="MS Gothic")
end

# 理論
function bara_single_theory(N, p1, p2, R1, R2)
  bin1 = Binomial(N, p1)  # Y_1 の二項分布
  bin2 = Binomial(N, p2)  # Y_2 の二項分布

  pmf = Dict{Int,Float64}() # S_N の確率質量関数

  for y1 in 0:N
    for y2 in 0:N
      s = R1 * y1 + R2 * y2
      p = pdf(bin1, y1) * pdf(bin2, y2) # Julia は pmf を関数で取得できる
      pmf[s] = get(pmf, s, 0.0) + p # 同じsが既にあれば足し合わせる
    end
  end
  return pmf
end

# 実験
function bara_single_simulate_once(R1::Int, R2::Int)
  # 当たり番号を独立に生成（1回の試行に対して1回だけ）
  winning_1digit = rand(0:9)  # 300円の当たり：下1桁
  winning_2digit = rand(0:99)  # 3000円の当たり：下2桁

  # ランダムな2桁を生成
  last_digit = rand(0:9)
  second_digit = rand(0:9)
  two_digit = last_digit + 10 * second_digit

  prize = 0
  # 下1桁が一致するかチェック
  if last_digit == winning_1digit
    prize += R1
  end
  # 下2桁が一致するかチェック
  if two_digit == winning_2digit
    prize += R2
  end

  return prize
end

function bara_single_simulate_set10(R1, R2)
  total_prize = 0
  for _ in 1:10
    total_prize += bara_single_simulate_once(R1, R2)
  end
  return total_prize
end

function bara_single_simulate(set_N::Int, R1, R2)
  result = Dict{Int,Float64}()

  for _ in 0:set_N
    prize = bara_single_simulate_set10(R1, R2)
    result[prize] = get(result, prize, 0) + 1
  end
  return result
end

kuji_N = 10 # くじを引く回数
kuji_p1 = 1 / 10 # 当たり1確率:下1桁一致
kuji_p2 = 1 / 100 # 当たり2確率:下2桁一致
kuji_R1 = 300 # 当たり1報酬
kuji_R2 = 3000 # 当たり2報酬

set_N = Int(1e7) # 実験回数

bara_single_theory_pmf = bara_single_theory(kuji_N, kuji_p1, kuji_p2, kuji_R1, kuji_R2);
bara_single_simulate_result = bara_single_simulate(set_N, kuji_R1, kuji_R2);

xticks = sort(collect(keys(bara_single_theory_pmf)))
bara_single_theory_prob = [bara_single_theory_pmf[s] for s in xticks]

# ヒストグラム
p = bar(xticks, bara_single_theory_prob,
  label="バラ単券10枚（理論）",
  xlabel="合計当たり額",
  ylabel="確率",
  legend=:right,
  bar_width=200,
  alpha=0.5,
  # xticks=collect(0:300:maximum(x)),
  xformatter=x -> string(Int(x)),
)

bar!(xticks, [get(bara_single_simulate_result, s, 0.0) / set_N for s in xticks],
  label="バラ単券10枚（実験）",
  xlabel="合計当たり額",
  ylabel="確率",
  legend=:right,
  bar_width=200,
  alpha=0.5,
  color=:red,
)
savefig(p, "imgs/bara_single_comparison.png")

theory_E = sum([s * p for (s, p) in bara_single_theory_pmf])
theory_Var = sum([(s - theory_E)^2 * p for (s, p) in bara_single_theory_pmf])
theory_SD = sqrt(theory_Var)

simulate_E = sum([s * count / set_N for (s, count) in bara_single_simulate_result])
simulate_Var = sum([(s - simulate_E)^2 * count / set_N for (s, count) in bara_single_simulate_result]) # 母分散
simulate_SD = sqrt(simulate_Var)

# md表として出力
println("| 項目 | 理論 | 実験 |")
println("| :----: | ----: | ----: |")
for row in eachrow(
  DataFrame(
    "項目" => ["期待値", "分散", "標準偏差"],
    "理論" => [theory_E, theory_Var, theory_SD],
    "実験" => [simulate_E, simulate_Var, simulate_SD],
  )
)
  println("| $(row.項目) | $(round(row.理論, digits=3)) | $(round(row.実験, digits=3)) |")
end

【図】Juliaによる理論と実験(シミュレーション)の結果の比較

項目	理論	実験
期待値	600.0	600.658
分散	972000.0	974096.243
標準偏差	985.901	986.963

10の7乗回、10枚の合計当たり額をシミュレーションした結果、ぴったり理論と重なった。

バラ券10枚組

バラ券10枚組の場合、単券の場合と違い10枚のくじの下1桁は0〜9が1枚ずつ入るため、以下のような報酬の決まり方となる。

10枚組をまとめて考える。

• 下1桁の300円の当たり

  10枚のくじの下1桁は0〜9が1枚ずつ入るため、必ず1枚が当たり、残りの9枚ははずれとなる。

• 下2桁の3000円の当たり

  1/100の確率で当たるくじが10枚あるので、1/100 x 10 = 1/10の確率で1枚が当たる。
  このとき、他の9枚は下1桁の数字が違うので全てはずれとなることに注意。
  したがって、1/10の確率で3000円の報酬、9/10の確率で0円の報酬となる。

したがって、バラ券10枚組を1セット購入したときの合計報酬額 $S$ は、300円か3300円の2通りで、それぞれ確率(分布)は、

$$
\begin{align*}
P(S = 300) &= \frac{9}{10} \\
P(S = 3300) &= \frac{1}{10}
\end{align*}
$$

であり、期待値 $E$ 、 分散 $\mathrm{Var}$ 、 標準偏差 $\mathrm{SD}$ は、

$$
\begin{align*}
E[S] &= \sum_i s_i p_i \\
&= 300 \times \frac{9}{10} + 3300 \times \frac{1}{10} \\
&= 600 \\
\mathrm{Var}(S) &= \sum_i (s_i - E[S])^2 p_i \\
&= (300 - 600)^2 \times \frac{9}{10} + (3300 - 600)^2 \times \frac{1}{10} \\
&= 90,000 \times \frac{9}{10} + 7,290,000 \times \frac{1}{10} \\
&= 810,000 \\
\mathrm{SD}(S) &= \sqrt{\mathrm{Var}(S)} \\
&= \sqrt{810,000} \\
&= 900
\end{align*}
$$

となる。

連番10枚組

連番10枚組の場合、以下のような報酬の決まり方となる。

• 下1桁の300円の当たり

  10枚のくじの下1桁は0〜9が1枚ずつ入るため、必ず1枚が当たり、残りの9枚ははずれとなる。

• 下2桁の3000円の当たり

  10の位が10枚のくじで共通で、1の位は10枚のくじで0〜9まで網羅しているため、10の位の数字が当たり番号の10の位の数字と一致するかどうか報酬が決まることになり、1/10の確率で3000円の報酬、9/10の確率で0円の報酬となる。

すなわち、10枚単位で考えれば、バラ券10枚組と連番10枚組の合計当たり額の確率分布は全く同じと言える。

(今回のように下数桁の当たりのみ考えている場合に限る、高額当選にあるような前後賞を考えると確率分布は変わるはず。)

まとめ

確率とは上手くできているもので、買い方を変えようと期待値は変わらない。
だが、バラ単券10枚の購入は、合計当たり額が0円となるリスクを生む一方、収支+30000円というロマンもくれる。
300円か3300円の2つに未来を限るのではなく、300区切りの沢山の可能性を見出そうではないか。

きっと標準偏差の大きさは夢の大きさと言っても過言ではない――。

※ 対面販売所でバラ単券をバラバラに注文した際に生じる迷惑・トラブル等には一切責任を負いません。ネット注文 とかにしておきましょう。