統計学を利用して、ギャンブルで勝てる賭け方をpythonで予測してみよう

分析屋

2024年10月17日 18:37

はじめに

分析屋の小泉と申します。

これまで統計学の基礎的な概念や考え方を重点的に解説してきました。
でも、統計検定に向けて勉強するのも大事ですが、
数式を覚えたり理論を覚えても、実際にどんなことに使えるかが分からないと面白くないですよね。

そこで今回は、これまでとはちょっと趣向を変えて、
以前解説したベルヌーイ試行を前提とするギャンブルを数理モデルに落とし込み、pythonで最適化することで、
どのような状況ならどれくらい所持金を増やせるかを予測してみます。

本記事では、前提の設定から統計学を用いて数理モデル化し、pythonでのコーディングと結果の解釈まで行います。

作成したコードは本記事の下層に残しますので、
ライブラリさえインストールしていればすぐ再現できますし、自由にいじってシミュレーションで遊んでみてください。

今後も内容は深まっていきますが、統計学の理解やモチベーション向上の一助となれば幸いです。

モデルの前提

今回検討するのは、ケリー基準（Kelly Criterion）と呼ばれる理論をアレンジしたものです。

これは、投資において勝率が分かっている状況で、繰り返し試行することで試算の最大化を狙うときの投資額を求める理論です。
実際はかなり複雑な計算が必要ですが、今回はこれまでに解説した内容が使えるように内容を簡略化して数理モデルに置き換えます。

1回ごとに勝つか負けるかの二択で結果が得られるため、ベルヌーイ試行ともいえますね。

ギャンブルの勝てる確率は事前処理やギャンブル中の手続きに問わず一定とし、想定できるものとします。
1回の試行が次回以降の確率に影響を及ぼさない、という意味では統計学的には復元抽出と呼ばれますね。

数式でモデルを表す

モデルの設計にあたり、各変数を定義します。

$${p}$$・・・勝てる確率

$${n}$$・・・ギャンブルを行う回数（試行数）

$${n_W}$$・・・勝つ回数

$${f}$$・・・1回の賭け金（所持金に対する割合）

$${M_0}$$・・・最初の所持金

$${M_n}$$・・・n回試行後の所持金

数が多いので、例を考えてみます。
例えば、最初の所持金$${M_0}$$を10000円、各回のベット（賭け金）を500円･･･$${f=500\div10000=0.02}$$としてみましょう。

ゲームに勝った場合は、

$${10000-500=9500}$$･･･ベットの支払い
$${M_1=9500+(500\times2)=10500}$$･･･ベットの2倍を賞金として入手

となります。
負けた場合は賞金はないため、

$${10000-500=9500}$$･･･ベットの支払い
$${M_1=9500}$$･･･ベットは帰ってこないのでそのまま

となります。

さて、変数を定義しましたが、
賭け金を割合で表すことで、所持金の推移を数式で簡単に表すことが出来ます。
最初の賭け金に対して、勝てば所持金はfの割合だけ増えるし、負ければ所持金はfの割合だけ減ってしまうからです。

先の例から具体的に考えると、$${f=0.05}$$としたとき、
勝った場合は$${M_0=10000}$$、$${M_1=10500}$$ですから、

$${M_1=M_0+M_{0}f=M_0\times(1+f)}$$

となります。
一方で、負けた場合は$${M_0=10000}$$、$${M_1=9500}$$ですから、

$${M_1=M_0-M_{0}f=M_0\times(1-f)}$$

となり、まとめて次のような数式で表されます。

$${M_1=M_0\times(1{\pm}f)}$$

同様に、2回目の試行後の所持金は、
1回目の所持金の割合から賭け金・取り分が得られるため、以下のように表されます。

$${M_2=M_1\times(1{\pm}f)}$$

また、$${M_1}$$は$${M_1=M_0\times(1\pm{f})}$$であったので、

$${M_2=M_0\times(1{\pm}f)\times(1{\pm}f)}$$

と言えます。

これを続けていくと、n回目の試行後の所持金は、最初の所持金から以下のように表せます。

$${M_n=M_{n-1}\times(1{\pm}f)}$$
$${M_n=M_0\times(1{\pm}f)^{n}}$$

なお、$${(1{\pm}f)^n}$$という表し方ですが、n回の試行には勝ちも負けも含まれています。

勝った回数$${n_W}$$ならば、負けた回数は$${n-n_W}$$で表されるため、
$${n=n_W+(n-n_W)}$$
といえます。

さらに、勝てる確率を$${p}$$と定義しているため、$${n_W}$$は$${p}$$を用いて、

$${n_W=np}$$

と表せます。

以上のことから、
$${\pm}$$のうち勝ったときは$${+}$$、負けたときは$${-}$$となるため、

$${M_n=M_0\times(1{+}f)^{n_W}\times(1{-}f)^{n-n_W}}$$
$${M_n=M_0\times(1{+}f)^{np}\times(1{-}f)^{n(1-p)}}$$

と表せました。

この2択の考え方は、統計学におけるベルヌーイ試行の理解があることで簡単に導き出すことが出来るはずです。
ベルヌーイ分布・二項分布の確率質量関数と似た構造をとっていることからも、統計学の理解がモデル構築の大前提となることをご理解頂けると思います。

なお、これはケリー基準の勉強をする方に向けた余談ですが、
勝てる確率$${p}$$に対して、$${1-p}$$で表される確率（=負ける確率）をオッズといいます。
ギャンブルなどで「オッズが〇〇倍」というのは、勝てそうな確率に対して負けそうな確率が何倍あるかを示しているんですね。
ですので、「このゲームはオッズが4倍」と言われたら、

$${\frac{1-p}{p}=4}$$
$${1-p=4p}$$
$${p=0.2}$$

となります。
実際は他の参加者がどこに賭けているかによって計算されますが、
少なくとも自分が勝つ確率$${p}$$に対してオッズ$${1-p}$$（自分が負ける確率）もギャンブルの世界ではよく用いられるようです。

このあとは$${p}$$ばかり出てきますが、オッズという言葉でイメージがつきやすい方は、こちらも頭に入れておくと良いかも知れません。

指数関数的成長率の導入

さて、ここまででn回の試行後の所持金を求める数理モデルが作成できました。
今回の結論として求めたいのはギャンブルでの最終的な勝敗ですので、
一般化すると、最終的な$${M_n}$$が$${M_0}$$より大きいことを勝利、小さいことを敗北とすると定義できます。

ということは、$${M_n-M_0>0}$$を満たすようにすればよい・・・と考えられますが、
実はそれでは不十分です。

というのも、fは割合を表すため、ベット額は試行前の所持金によって毎回変わってしまいます。
すなわち、$${M_n}$$の値は$${M_0}$$の値によって変動してしまうのです。

今回の議題では初期の所持金額を問わず、ギャンブル前の所持金に関わらずに利益を最大化することを目的としているため、
単純な$${M_n}$$と$${M_0}$$の絶対比較は望ましくありません。

では、どのようにして数理モデルを評価するべきでしょうか。
ここで導入するのが、平均変化率（成長率）という考え方です。
基準となる値を$${A}$$とし、$${n}$$期後の値を$${B}$$としたとき、以下の等式が成立します。

$${A\times(1+r)^n=B}$$

基本的に$${n}$$は1であることが多く、この場合は$${r}$$が増減の割合を示すこととなります。

$${n}$$は基準となるステップ数であり、$${A}$$と&&{B}&&の間を等間隔に分割すれば任意に定めることができます。
例えば1年後の売り上げなら年単位だと$${n=1}$$、四半期単位なら$${n=4}$$、月単位なら$${12}$$とステップ幅に応じてnが定まり、
それが何ステップあったかで売上の粒度を変えることが出来ます。

このため、$${r}$$は各ステップの平均変化率と呼ばれます。

平均変化率の利用のメリットとしては、
$${r}$$が0以上ならば$${B{\geq}A}$$が成立し、
$${r}$$が0未満であれば$${B{<}A}$$が成立することから、
$${r}$$のみで初期値に関わらず変化量を評価できることがあります。

今回の場合、初期値Aは$${M_0}$$となり、試行後の値Bは$${M_n}$$となります。
先ほどの成長率の等式を変形し、

$${(1+r)^n=\frac{M_n}{M_0}}$$
$${1+r=(\frac{M_n}{M_0})^{\frac{1}{n}}}$$

自然対数をとることで$${\frac{1}{n}}$$を外に出し、

$${\ln{(1+r)}=\frac{1}{n}\ln{(\frac{M_n }{M_0})}}$$

左辺の$${\ln{(1+r)}$$は指数関数的成長率と呼ばれ、$${G}$$として表し、

$${G=\frac{1}{n}\ln{(\frac{M_n }{M_0})}}$$

という等式が作成できました。

これにより、このギャンブルでは
$${G>0}$$のとき$${M_n>M_0}$$となり、最終的に勝つと、
$${G<0}$$のとき$${M_n<M_0}$$となり、最終的に負けると結論づけられます。

$${M_n}$$は$${M_n=M_0\times(1+f)^{np}\times(1-f)^{n(1-p)}}$$となりますから、

pを任意に設定し、
所持金のベット割合であるfを$${0<f{\leq}1}$$の範囲でGをプロットし、
最もGが大きいf(便宜上$${f_{G_{max}}}$$と呼ぶ)をモデルの最適化から求める
という手順で進めていきます。

なお、ここでいう「所持金が増える・減る」という計算結果は期待値であり、必ずこの結果になるというものではありません。
厳密に期待値の定義から考えると、同じ状況(勝てる確率がpのベルヌーイ分布確率モデルに従うゲームにおいて、所持金のうち$${f\times100%}$$)を$${n}$$回賭ける状況）が極限に多数想定するとき、
所持金の指数関数的成長率$${G}$$の平均がどのくらいになるかを表すものです。

ちょっとややこしいですが、グラフで現れる線は最も可能性が高い結果という解釈にとどめてくださいね、ということです。

スクリプトの作成

ここからはスクリプトの作成になります。
対数を利用するのでnumpy
グラフの可視化のためにmatplotlibのpyplot（お好みでjapanize-matplotlibも）
モデルの最適化を検討するためscipyのoptimize
を利用します。

Pythonスクリプトの作り方は人によってある程度流派が分かれると思いますが、
今回はここまでで作成した数理モデルを利用し、
①numpyを用いて数理モデルを関数として作成
②fの最適化をscipy.optimizeで実行
③シミュレーションを実行して結果をmatplotlibで可視化
という流れで行います。
まずは、前段階としてライブラリのインポートと変数の設定を行います。
Numpyのlinespaceを利用して、fのシュミレーションに使う値を用意しておきましょう。

最小値は0（ベットなし）、最大値は1（所持金オールイン）とし、要素数は線に見えるよう最低でも1000くらいあるとよいでしょう。
今回は大変な計算もないので、100000くらいまで欲張ります。

また、pは可視化する際に比較しやすいよう、配列で指定しておきます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# matplotlibの日本語化、お好みで
import japanize_matplotlib as jp_plt
jp_plt.japanize()

# 各変数の定義
f_val = np.linspace(0, 1, 100000)
p_val = [0.2, 0.5, 0.8] # 想定される勝率
n_val = 1

次に、数理モデルから最適化指標となるGを求めて返す関数calc_Gを作ります。
先に求めた数式をそのまま入力するだけですね。

なお、対数の表記は分野によって異なりますが、Numpyは数学的な表現で書かれているため、数式上ではlnを用いていますが、コード上ではlogを利用します。

また、対数は負の値や0でエラーとなります。
対数関数は極限まで0に近づけると負の無限大に発散するため、エラーハンドリングとして負の無限大を返すようにしておくとよいでしょう。

# G(平均変化率)の計算
def calc_G(f, p, n):
    try:
        M_n = (1 + f)**(n * p) * (1 - f)**(n * (1 - p))
        G = (1 / n) * np.log(M_n)
        return G
    # 極限値より、エラー値は-infで処理する
    except ValueError:
        return -np.inf

次に、fの最適化を行う関数f_optimalを作成します。

scipyのminimizeを用いて、先のcalc_Gを-にして実行し最小の値を探索します。
初期値は0だとエラーになるので、0以外の小さい値（今回は0.01）を設定しておきましょう。
もし見つからなければNoneを返します。

# f(ベット値)が最大（-calc_Gが最小）となるようなfの最適化
def f_optimal(p):
    rslt = minimize(lambda f: -calc_G(f, p, n_val), 0.01, bounds=[(0, 1)])
    if rslt.success:
        return rslt.x[0]
    else:
        return None

続いて、ここまでの関数の呼び出しとmatplotlibによる可視化の設定を行うrun_sim関数を作成します。

Pもnも配列にしているので、for-in文を用いて連続処理を行います。

matplotlibの設定は適当ですが、japanize-matplotlibを利用しない場合はタイトルやx,yラベルから日本語を削除してください。

# 配列ごとに最適化を実行
def run_sim(p_val, n_val, f_val):
    for p in p_val:
        G_val = [calc_G(f, p, n_val) for f in f_val]
        plt.plot(f_val, G_val, label=f'p={p}, n={n_val}')

        # 最適化したfをf_optとして計算、図示
        f_opt = f_optimal(p)
        if f_opt is not None:
            G_opt = calc_G(f_opt, p, n_val)

            # 最適化点の出力設定
            plt.scatter(f_opt, G_opt, color='black')
            plt.text(f_opt, G_opt,  f'f_opt = {f_opt:.2f}\nG_opt = {G_opt:.2f}', fontsize=8, verticalalignment='bottom')

    # グラフ作成
    plt.xlabel('f (ベット額割合；Betting Money Fraction)')
    plt.ylabel('G (対数平均変化率；Growth Rate)')
    plt.title('ケリー基準の最適化検証；Kelly Criterion Optimization')
    plt.legend()
    plt.grid(True)
    plt.show()

最後に、変数を引数に指定して実行を書けば完成です。

# シミュレーション実行
run_sim(p_val, n_val, f_val)

結果の解釈・結論

スクリプトを実行した結果が以下のようになります。

n_val=1で出力。最適化された値は黒でプロットされている

横軸は所持金に対するベット数割合f、縦軸は試行後の所持金の増加率を対数で表したGです。

グラフの見方として、
縦軸のGは、自然対数のためG=1のとき試行後に初期の所持金がネイピア数$${e=2.718･･･}$$倍に増えると解釈されます。

一方で、0の時は$${M_n/M_0=1}$$となり所持金に変動がなく、
負の値をとる場合は$${M_n/M_0<1}$$となり、試行後の方が所持金が減ってしまう（敗北）と解釈されます。

すなわち、$${p=0.2}$$や$${p=0.5}$$のときは、最も$${G}$$が大きくなる$${f_{opt}=0}$$のとき$${G_{opt}=0}$$であるため、
勝てる確率が0.5未満ならば、所持金の0%を賭けたとき最もゲーム後の所持金の期待値が多い･･･すなわち、賭けるだけ所持金は減る可能性が高いということになります。

また、$${p=0.8}$$のときは$${f_{opt}=0.6}$$で$${G_{opt}=0.19}$$のため、
勝てる確率が80%のときは、所持金の60%をベットすると、得られる期待値は$${\exp(0.19)=1.21}$$倍になる可能性が高いといえます。

「8割の確率で勝てるギャンブル」と聞くとつい食いついてしまいますが、大抵1.2倍程度にしかならないと考えると、なんだか世知辛いですね。

では、n=1000に増やして試行回数を増やしたらどうなるでしょうか。
実際はそんな回数ギャンブルにつぎ込むなよという話ではありますが、試行回数を増やせば「確率は収束する」や「負けた分も勝った分も結局取り戻せる」みたいな印象がありますし、
シミュレーションなら実際の所持金が減るわけでもないのでガンガン試してみましょう。

コードのnの部分を以下のように1000に変えて実行してみます。

n_val = 1000

n_val=1000で出力。線の形はn=1と比べて変わっているけれど･･･

先のn=1と比較して、グラフの形は結構変わりました。
n=1のときはf=1.0付近で一気にGが低下した印象ですが、
n=1000のときはもっと前（fが小さいところ）からGが低下しています。

また、最適化されたf_optにも注目してみましょう。
n=1のときのf_optと全く同じ値を示し、G_optも同じ値を得ています。

このことから分かるのは、いくら試行回数$${n}$$を増やしても、この数理モデルにおいてギャンブルに勝てる（所持金を増やす）方法は、
想定される勝利確率$${p}$$のみに依存するということです。

また、$${p\leq0.5}$$のときは、そもそも賭けたら所持金が減る可能性の方が高いということになります。

もちろん、実際のギャンブルでは勝利と敗北の2択ではないですし、
「低確率に挑むから燃えるんだ」という賭博黙示録のような人もいると思います。

とはいえ、「統計学的にはこう」という実証もでき、pythonでのモデル構築の練習にもなる良い例として紹介しました。

本記事の末尾にスクリプトを全文記載しますので、
是非ご自身の環境で確率や試行回数を変化して遊んでみてください。

ところで、$${f_opt}$$とpには相関があり、簡単な関数で求めることが出来ます。
pythonでのシミュレーションのほか、微分・二回微分によって手計算でも証明が可能です。

また、統計的仮説検定を活用すると、$${p>0.5}$$で少なくとも損をしない$${f}$$も求められるようになります。

いずれは解説することになると思いますが、興味がある方は是非遊んでみてください。

終わりに

今回は、これまでの統計学の解説を行った内容をモデル化し、pythonを用いてシミュレーションを行いました。

なお、今回の内容は次回以降に解説する大数の法則や正規分布へ大いに活用出来るものです。
nが大きくなるとカーブが大きくなるのは何故か？などという視点にも着目すると、今後の理解にも深く繋がると思います。

統計学は数式が多く、データサイエンス領域との兼ね合いが分かりづらい部分が大いにあります。
今回のように簡単なコードで可視化し、シミュレーションすることで統計学への興味・関心が増えてくれれば非常に嬉しいです。

おまけ（今回作成したpythonコード）

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# matplotlibの日本語化、お好みで
import japanize_matplotlib as jp_plt
jp_plt.japanize()

# 各変数の定義
f_val = np.linspace(0, 1, 100000)
p_val = [0.2, 0.5, 0.8]  # 想定される勝率
n_val = 1 # 試行回数

# G(平均変化率)の計算
def calc_G(f, p, n):
    try:
        M_n = (1 + f)**(n * p) * (1 - f)**(n * (1 - p))
        G = (1 / n) * np.log(M_n)
        return G
    # 極限値より、エラー値は-infで処理する
    except ValueError:
        return -np.inf

# f(ベット値)が最大（-calc_Gが最小）となるようなfの最適化
def f_optimal(p):
    rslt = minimize(lambda f: -calc_G(f, p, n_val), 0.01, bounds=[(0, 1)])
    if rslt.success:
        return rslt.x[0]
    else:
        return None


# 配列ごとに最適化を実行
def run_sim(p_val, n_val, f_val):
    for p in p_val:
        G_val = [calc_G(f, p, n_val) for f in f_val]
        plt.plot(f_val, G_val, label=f'p={p}, n={n_val}')

        # 最適化したfをf_optとして計算、図示
        f_opt = f_optimal(p)
        if f_opt is not None:
            G_opt = calc_G(f_opt, p, n_val)

            # 最適化点の出力設定
            plt.scatter(f_opt, G_opt, color='black')
            plt.text(f_opt, G_opt,  f'f_opt = {f_opt:.2f}\nG_opt = {G_opt:.2f}', fontsize=8, verticalalignment='bottom')

    # グラフ作成
    plt.xlabel('f (ベット額割合；Betting Money Fraction)')
    plt.ylabel('G (対数平均変化率；Growth Rate)')
    plt.title('ケリー基準の最適化検証；Kelly Criterion Optimization')
    plt.legend()
    plt.grid(True)
    plt.show()

# シミュレーション実行
run_sim(p_val, n_val, f_val)

ここまでお読みいただき、ありがとうございました！
この記事が少しでも参考になりましたら「スキ」を押していただけると幸いです！

これまでの記事はこちら！

株式会社分析屋について

弊社が作成を行いました分析レポートを、鎌倉市観光協会様HPに掲載いただきました。

ホームページはこちら。

noteでの会社紹介記事はこちら。

【データ分析で日本を豊かに】
分析屋はシステム分野・ライフサイエンス分野・マーケティング分野の知見を生かし、多種多様な分野の企業様のデータ分析のご支援をさせていただいております。「あなたの問題解決をする」をモットーに、お客様の抱える課題にあわせた解析・分析手法を用いて、問題解決へのお手伝いをいたします！

【マーケティング】
マーケティング戦略上の目的に向けて、各種のデータ統合及び加工ならびにPDCAサイクル運用全般を支援や高度なデータ分析技術により複雑な課題解決に向けての分析サービスを提供いたします。

【システム】
アプリケーション開発やデータベース構築、WEBサイト構築、運用保守業務などお客様の問題やご要望に沿ってご支援いたします。

【ライフサイエンス】
機械学習や各種アルゴリズムなどの解析アルゴリズム開発サービスを提供いたします。過去には医療系のバイタルデータを扱った解析が主でしたが、今後はそれらで培った経験・技術を工業など他の分野の企業様の問題解決にも役立てていく方針です。

【SES】
SESサービスも行っております。

この記事が参加している募集

#数学がすき

3,435件

#やってみた

42,500件