乙Py先生のプログラミング教室
初学者のためのプログラミング学習サイト

プログラミング

今回はブラック-ショールズ方程式で
オプション価格を計算する方法についてです


解説動画はこちら




ブラック-ショールズ方程式とは

デリバティブの価格づけに現れる
偏微分方程式のことです。

デリバティブは株式や債券、通貨、商品などの
原資産から派生した金融商品の総称のことです。

ブラックショールズ方程式は主に
ヨーロピアン・オプション(満期時のみ行使可) の
理論価格を求めるためのモデルになっています。


オプション取引とは

金融商品のデリバティブの一種で
ある原資産について、あらかじめ決められた
将来の一定の日または期間において
事前に定めた権利行使価格で取引できる
「権利」のことです

原資産を
買う権利についてのオプションを「コールオプション」
売る権利についてのオプションを「プットオプション」
と呼んでいます。



日経225オプション

日経225オプションの価格情報があります。


オプション取引の価格計算ツール(日本取引所)

価格計算ツールもついているので
今回はこの計算方法をPythonで再現します。


ブラックショールズの計算式について

モデルの前提条件

•    価格は幾何ブラウン運動に従う
•    市場は完全で、裁定取引がない
•    無リスク金利は一定
•    ボラティリティ(変動率)は一定
•    配当は考慮しない(配当あり版も考慮は可)


使用する変数

S 現在の価格(Spot Price)
K 権利行使価格(Strike Price)
T 満期までの残存期間(年単位)
r 無リスク金利(年利)
σ 価格のボラティリティ(年率標準偏差)

ブラック-ショールズの数式

コールオプション価格(買う権利):
スクリーンショット 2025-07-12 17.06.54


プットオプション価格(売る権利):
スクリーンショット 2025-07-12 17.06.59

こんな感じの計算式になっていますが
途中、d1,d2というものが必要になります。

d_1, d_2 の定義
スクリーンショット 2025-07-12 17.08.28

なお、N(x) は 標準正規分布の累積分布関数(CDF)になります。


オプション計算のPythonコード

下記の関数で、オプション価格を計算できます。
import numpy as np
from scipy.stats import norm

def black_scholes_option_price(S, K, T, r, sigma, option_type='call'):
    """
    ブラック-ショールズ方程式によるオプション価格の計算
    S: 現在の株価
    K: 権利行使価格
    T: 残存期間(年単位)
    r: 無リスク金利(年利)
    sigma: ボラティリティ(年率標準偏差)
    option_type: 'call' か 'put'
    """
    d1 = (np.log(S / K) + (r + sigma**2 / 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    if option_type == 'call':
        price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    elif option_type == 'put':
        price = K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    else:
        raise ValueError("option_type must be 'call' or 'put'")
    return price

実際に値を入れて計算する場合はこのようになります。
日数で計算したい場合は年単位なのでうまく合わせます。
# パラメータ例
S = 39425       # 現在の価格
K = 40000       # 権利行使価格
T = 62/365      # 満期まで62日
r = 0.0         # 無リスク金利 0%
sigma = 0.1883  # ボラティリティ 18.83%

call_price = black_scholes_option_price(S, K, T, r, sigma, option_type='call')
put_price = black_scholes_option_price(S, K, T, r, sigma, option_type='put')

print(f"Call Option Price: {call_price:.4f}")
print(f"Put Option Price : {put_price:.4f}")
Call Option Price: 963.0367
Put Option Price : 1538.0367




まとめ

この計算結果をどう応用するかですが
オプション取引の戦略を決定するために使います。

市場価格 > 理論価格

オプションが「割高」
プットまたはコールを「売る(ショート)」

市場価格 < 理論価格

オプションが「割安」
プットまたはコールを「買う(ロング)」


ということになります。

この辺りを実際の売買手順フローに落とし込むと
① 必要データ収集(S, K, T, r, σ)
        ↓
② BSモデルで理論価格算出(call/put)
        ↓
③ 実際の市場価格と比較(割高・割安を判断)
        ↓
④ 戦略選定(買う/売る or 戦略組み合わせ)
        ↓
⑤ 証券口座で売買、モニタリング

という感じになるので
自動取引ができるんじゃ無いかと
画策しているところです。

今回はブラックショールズ方程式による
オプション価格の計算方法についてでした。

今回はここまでです
それでは。






今回は地震データの可視化についてです
最近地震が多いので、plotlyで可視化してみました。


解説動画はこちら








データの入手先


気象庁の地震データベース

コードを動かしてみたい方は
CSVがダウンロードできるようなので
手元にダウンロードしてみてください。



データの読み込み

こちらのコードはGoogle Colabで動くようになっています。
動かしたい場合は Colabの画面左メニューから
フォルダマークをクリックすると
ファイル置き場が見えるので、そこにCSVファイルを
ドラッグなどで配置します。

ファイルを読み込みするには下記のコードです。

# 必要なライブラリのインポート
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
import numpy as np
import warnings
warnings.filterwarnings('ignore')

file_path = "/content/地震リスト.csv"
df = pd.read_csv(file_path)
df.head()

df変数にデータが読み込まれると思います。



データの整形


CSVそのままのデータでは
うまく可視化が行えないため
可視化用にデータを加工する必要があります。

下記のコードを実行すると
可視化用のカラムなどが追加されます。


# データの前処理
s_order = ['震度1', '震度2', '震度3', '震度4', '震度5弱', '震度5強', '震度6弱', '震度6強', '震度7']
target = "トカラ列島近海"

def preprocess_data(df):
    # 緯度・経度の変換(度分秒から度へ)
    def convert_coordinate(coord_str):
        # 例: "29°28.2′N" -> 29.47
        if '°' in coord_str and '′' in coord_str:
            parts = coord_str.replace('N', '').replace('E', '').replace('S', '').replace('W', '')
            degree_part = parts.split('°')[0]
            minute_part = parts.split('°')[1].replace('′', '')
            return float(degree_part) + float(minute_part) / 60
        return float(coord_str)
    
    df['緯度_数値'] = df['緯度'].apply(convert_coordinate)
    df['経度_数値'] = df['経度'].apply(convert_coordinate)
    
    # 深さの数値化
    df['深さ_数値'] = df['深さ'].str.replace(' km', '').astype(float)
    
    # 震度の数値化
    s_mapping = {
        '震度1': 1, '震度2': 2, '震度3': 3, '震度4': 4, '震度5弱': 5,
        '震度5強': 5.5, '震度6弱': 6, '震度6強': 6.5, '震度7': 7
    }
    df['最大震度_数値'] = df['最大震度'].map(s_mapping)
    df['最大震度'] = pd.Categorical(df['最大震度'], categories=s_order, ordered=True)
    return df

# データの前処理を実行
df = preprocess_data(df)
df = df[df["震央地名"]==target]
df.head()



時系列の可視化


最初は時系列で
どれだけ地震が発生しているのかを
見てみましょう

日別、震度別で時系列で地震の回数を
表示してみます。

# 時系列分析(日付別)
df['日付'] = pd.to_datetime(df['地震の発生日'])
daily_shindo_counts = df.groupby(['日付', '最大震度']).size().reset_index(name='地震回数')

# 震度の順序を維持
daily_shindo_counts['最大震度'] = pd.Categorical(
    daily_shindo_counts['最大震度'], 
    categories=s_order, 
    ordered=True
)

# 震度別横並び棒グラフ
fig_time = px.bar(
    daily_shindo_counts,
    x='日付',
    y='地震回数',
    color='最大震度',
    title='日付別地震発生回数(震度別)',
    labels={'地震回数': '地震発生回数', '日付': '発生日'},
    category_orders={'最大震度': s_order},
    barmode='group'  # 横並びに表示
)

# レイアウトの調整
fig_time.update_layout(
    xaxis_tickangle=-45,
    height=600,
    width=1000,
    legend=dict(
        orientation="v",
        yanchor="top",
        y=1,
        xanchor="left",
        x=1.01
    )
)

fig_time.show()
スクリーンショット 2025-07-05 17.27.01

データは直近の1週間分しか無いようです。
7月前後で急増しているのが分かります。



箱ひげ図の可視化


次は箱ひげ図で
「トカラ列島近海」の詳細を見てみましょう。


震度ごとにマグニチュードをまとめると
このような感じになります。

# 震央地名 最大震度別での箱ひげ図

fig_box = px.box(df, 
                 x='震央地名', 
                 y='M',
                 color='最大震度',
                 title='震央地名別 マグニチュード分布(最大震度別)',
                 labels={'M': 'マグニチュード', '震央地名': '震央地名'},
                 category_orders={'最大震度': s_order})

fig_box.update_layout(
    xaxis_tickangle=-45,
    height=600,
    width=1000
)
fig_box.show()
スクリーンショット 2025-07-05 17.29.52
同じ震度でも、マグニチュードにはバラツキがありますが
マグニチュードが上がるにつれ、震度も大きくなっています。



地図表示

今後は地図で震度を表示してみましょう
OpenStreetMapを使用して、緯度軽度を用いて
地図に震度を表示させてみます。

# マップ表示用のデータ準備
df_map = df.copy()
df_map['緯度_数値'] = pd.to_numeric(df_map['緯度_数値'], errors='coerce')
df_map['経度_数値'] = pd.to_numeric(df_map['経度_数値'], errors='coerce')
df_map['M'] = pd.to_numeric(df_map['M'], errors='coerce')
df_map = df_map.dropna(subset=['緯度_数値', '経度_数値', 'M'])
df_map['最大震度'] = pd.Categorical(df_map['最大震度'], categories=s_order, ordered=True)

# 緯度経度と最大震度・Mを用いたmap表示(OpenStreetMap使用)
fig_map2 = go.Figure()

# 震度別に色分けしてプロット
震度_colors = {
    '震度1': '#90EE90',  # 薄緑
    '震度2': '#FFD700',  # 金色
    '震度3': '#FFA500',  # オレンジ
    '震度4': '#FF6347',  # トマト色
    '震度5弱': '#FF4500',  # 赤オレンジ
    '震度5強': '#FF0000',  # 赤
    '震度6弱': '#8B0000',  # 濃い赤
    '震度6強': '#4B0082',  # インディゴ
    '震度7': '#000000'   # 黒
}

for 震度 in s_order:
    if 震度 in df_map['最大震度'].values:
        subset = df_map[df_map['最大震度'] == 震度]
        fig_map2.add_trace(go.Scattermapbox(
            lat=subset['緯度_数値'],
            lon=subset['経度_数値'],
            mode='markers',
            marker=dict(
                size=subset['最大震度_数値'] ** 2 ,  # サイズ調整
                color=震度_colors.get(震度, '#000000'),
                sizemode='diameter'
            ),
            text=subset['震央地名'],
            hovertemplate='%{text}
' + '緯度: %{lat:.2f}
' + '経度: %{lon:.2f}
' + 'マグニチュード: %{customdata[0]}
' + '最大震度: %{customdata[1]}
' + '', customdata=list(zip(subset['M'], subset['最大震度'])), name=震度 )) fig_map2.update_layout( title='地震発生位置(最大震度・マグニチュード)- OpenStreetMap', mapbox=dict( style="open-street-map", center=dict(lat=df_map['緯度_数値'].mean(), lon=df_map['経度_数値'].mean()), zoom=8 ), height=700, width=1000 ) fig_map2.show()
スクリーンショット 2025-07-05 17.34.15

今回の地震は、かなり局所的に起きていることが分かります。
特に震度5以上が発生した箇所がかなり近く
この近辺は注意が必要な地域に見えます。


散布図表示


最後に
マグニチュード、震度、深さを
散布図で見てみましょう

# 深さとマグニチュードの関係
fig_scatter = px.scatter(
    df,
    x='深さ_数値',
    y='M',
    color='最大震度',
    size='最大震度_数値',
    hover_name='震央地名',
    hover_data={
        '地震の発生日': True,
        '地震の発生時刻': True
    },
    title='深さとマグニチュードの関係(最大震度別)',
    labels={'深さ_数値': '深さ (km)', 'M': 'マグニチュード'},
    category_orders={'最大震度': s_order}
)

fig_scatter.show()
スクリーンショット 2025-07-05 17.36.27


どの深さでも、満遍なく起きているように見えますね
地震自体は、いつ何時、どこでも起きてしまいます。



まとめ

緯度経度が記載されているデータは
plotlyなどを用いれば地図表示が比較的
簡単に行うことができます。

可視化の表現の幅が広がるので
覚えておくと良いかもしれませんね


今日はここまでです
それでは




今回はドルコスト平均法による
積立投資のシミュレーションです

解説動画はこちら



ドルコスト平均法とは

価格が変動する金融商品に対して
一定金額を定期的に購入していく投資方法のことです

今回は毎月10万円積み立てるとして
それがどうなるかをシミュレーションしていきます。



eMaxisSlim米国株式S&P500

リンク先

S&P500指数(配当込み、円換算ベース)に
連動する運用成果を目指す投資信託
2018年7月3日に設定された商品です。

今回はこのデータを用いて
シミュレーションをしていきます。



eMAXIS Slim 米国株式 S&P500の価格推移

価格のデータを定義して
価格推移をだしてみます。

import numpy as np
import pandas as pd
import plotly.express as px
import seaborn as sns
import matplotlib.pyplot as plt

# データ定義
data = {
    2018 : [0,0,0,0,0,0,10038,10458,10709,11051,10188,10468],
    2019 : [8809,9854,10408,10565,10900,10031,10725,10978,10491,10890,11106,11680],
    2020 : [11873,11892,10825,9474,10655,11205,11465,11886,12707,12182,11767,12989],
    2021 : [13339,13407,14010,15217,15754,15924,16540,16712,17301,16677,18267,18005],
    2022 : [19291,18130,17601,19347,18798,18658,18052,19265,19394,18280,20283,19728],
    2023 : [17690,18714,19167,19388,20234,20666,22867,23260,23411,22911,22677,24155],
    2024 : [24154,25486,27473,28563,28573,29833,31693,29763,29789,29974,31336,32768],
    2025 : [33928,34065,32500,30512,28931,30867]
}

records = []
for year, prices in data.items():
    for month_idx, price in enumerate(prices):
        month = month_idx + 1
        if price == 0:
            continue
        date_str = f"{year}-{month:02d}-01"
        date = pd.to_datetime(date_str)
        records.append({"Date": date, "Price": price})

df = pd.DataFrame(records)
df = df.sort_values("Date")

# Plotlyでプロット
fig = px.line(df, x="Date", y="Price", title="Emaxis Slim value")
fig.show()
スクリーンショット 2025-06-21 16.43.01


良い感じに右肩上がりのようですが
この商品を毎月10万円ずつ積み立てたらどうなっていたでしょうか



2018年から2025年まで毎月積み立てた際のシミュレーション

毎月10万円
購入できる購入数量を算出して
2025年の最後の価格で
どれだけのリターンになっているかを算出します。
monthly_investment = 100000  # 毎月XX万円

# シミュレーション開始
total_units = 0  # 累計購入数量
total_invested = 0  # 累計投資金額
for year in sorted(data.keys()):
    for month_idx, price in enumerate(data[year], 1):
        if price == 0:
            continue
        units_bought = monthly_investment / price
        total_units += units_bought
        total_invested += monthly_investment

latest_year = max(data.keys())
latest_price = data[latest_year][-1]
current_value = total_units * latest_price
profit = current_value - total_invested
roi = (current_value / total_invested - 1) * 100

# 結果表示
print(f"総投資額: {total_invested:,.0f}円")
print(f"最終評価額: {current_value:,.0f}円")
print(f"損益: {profit:,.0f}円")
print(f"累計購入数量: {total_units:.4f}口")
print(f"リターン: {roi:.2f}%")
total_return = (100 + roi) / 100
years = 8
cagr = (total_return ** (1 / years)) - 1
print(f"年平均リターン: {cagr * 100:.2f}%")

総投資額: 8,400,000円
最終評価額: 16,458,895円
損益: 8,058,895円
累計購入数量: 533.2198口
リターン: 95.94%
年平均リターン: 8.77%


この価格推移のデータ上では
年平均8%超

7年で二倍近くの運用成績になっています。



ターンヒートマップ


今度は開始年から終了年までの
リターンヒートマップを出してみましょう。


開始と終了年を設定して
開始終了までのリターンを算出します。

years = list(data.keys())
monthly_invest = 100000
results = pd.DataFrame(index=years[:-1], columns=years[1:])
for start in years[:-1]:
    for end in years[years.index(start)+1:]:
        total_invest = 0
        total_units = 0
        for y in range(start, end+1):
            for price in data[y]:
                if price == 0:
                    continue
                units = monthly_invest / price
                total_units += units
                total_invest += monthly_invest
        # 最後の年の最後の価格で評価額計算
        final_price = [p for p in data[end] if p != 0][-1]
        final_value = total_units * final_price
        return_pct = (final_value - total_invest) / total_invest * 100
        results.loc[start, end] = return_pct

# ヒートマップ描画
results = results.astype(float)
plt.figure(figsize=(10, 7))
sns.heatmap(results, annot=True, fmt=".2f", cmap="coolwarm", center=0)
plt.title("dollar cost averaging method return heatmap (%)")
plt.xlabel("end year")
plt.ylabel("start year")
plt.show()
download

年単位だとどの年で買って売っても
プラスにはなっているようです。



まとめ

eMAXIS Slim 米国株式 S&P500のドルコスト平均法による積み立ては
年単位での運用であればマイナスは無いようです。

ただし、月単位だと乱高下があり
損する場合も有ったので、やはり長期積立で
見守るのが良いのではないかと思われます。

投資は自己判断で行うことが大切です。
その金融商品の良し悪しを判断し
投資するかどうかは人に言われるでなく
自分自身の判断で行うのが重要です。

投資判断を行うための材料作りとして
プログラムを用いたシミュレーションは
かなり役に立ちます。

データさえあればシミュレーションでき
結果の良し悪しから
投資判断の材料に使えるようになると思うので
プログラミングが出来るようになっていると
すぐに試すことができるのでお勧めです。

投資 x プログラミング
というテーマを今後も取り扱っていくので
両方できるようになりたい方は
要チェックしてみてください。

それでは


今回は金融庁の
NISA積立シミュレーターがやばいので
その解説です。


解説動画はこちら


 

金融庁の積立シミュレーター

これがそのシミュレーターです。

サイト

ここで
毎月の積立1万円
年利12%
40年間

で試算してみます。
スクリーンショット 2025-06-14 17.12.47


こんな感じになりました。

年利が高すぎますが
計算がわかりやすいようにして
計算してみた結果です。

40年後の資産は
11765万円になりました



Pythonでのシミュレーションコード


先程のシミュレーターと同じ計算を
Pythonでも計算してみましょう。

シミュレーションコード

monthly_rate = 0.12 / 12       # 月利(年利12%)
months = 40 * 12               # 総月数(480ヶ月)
monthly_payment = 10000       # 毎月の積立額(マイナス表記にする)

total = 0
for _ in range(months):
    total = total * (1 + monthly_rate) + monthly_payment

print(f"将来の資産総額: {total:,.0f} 円")
将来の資産総額: 117,647,725 円


金融庁のシミュレーターの結果は
11765万円 なので
ほぼ同じ結果になりました。


落とし穴

実はここに落とし穴があります!!!!

年利計算での正しいシミュレーション結果は
こうなります。
(年末に積立するパターン)

year_rate = 0.12            # 年利12%
years = 40                  # 総年数(40年)
year_payment = 120000    # 年の積立額

total = 0
for _ in range(years):
    total = total * (1 + year_rate) + year_payment

print(f"将来の資産総額: {total:,.0f} 円")
将来の資産総額: 92,050,970 円

先程の計算結果 : 117,647,725 円
とは大きく乖離してきます。


からくり

金融庁のサイトのシミュレーションは年利ではなく
年利を12ヶ月で割った月利の複利(年利 / 12)で計算しているようです。

そのため年利と月利では、計算結果が大きく異なってきます。

もし月利複利計算を正しくするとしたら
12乗して年利になるような値を月利にしないといけません。

月利 = (1 + 0.12) ** (1/12) - 1
= 0.009488792934583046 = 0.948%

年利から月利に直して計算するコードはこちら
monthly_rate = (1 + 0.12) ** (1/12) - 1  # 月利(年利12%)
months = 40 * 12                          # 総月数(480ヶ月)
monthly_payment = 10000                 # 毎月の積立額(マイナス表記にする)

total = 0
for _ in range(months):
    total = total * (1 + monthly_rate) + monthly_payment

print(f"将来の資産総額: {total:,.0f} 円")
将来の資産総額: 97,010,200 円

月利計算だとこのくらいの金額になります。




資産推移のシミュレーター

年利を月利換算しても正しく計算するようにしたものは
こんな感じになります。
import pandas as pd
import plotly.express as px

def plot_investment(rate_percent, years, payment_man):
    annual_rate = rate_percent / 100
    monthly_rate = (1 + annual_rate) ** (1/12) - 1
    months = years * 12
    monthly_payment = payment_man * 10000

    results, ci_total, am_total = [], 0, 0
    for month in range(1, months + 1):
        ci_total = ci_total * (1 + monthly_rate) + monthly_payment
        am_total += monthly_payment
        results.append({'月': month, '資産総額': ci_total, '積立金額': am_total})

    df = pd.DataFrame(results)
    fig = px.line(df, x='月', y=['資産総額', '積立金額'],
                  labels={'value': '金額(円)', 'variable': '項目'},
                  title=f'毎月{payment_man}万円積立(年利{rate_percent:.1f}%・{years}年)の資産推移')
    
    fig.update_layout(
        yaxis_tickformat=',',
        height=400
    )
    fig.show()

積立シミュレーション設定

#@title 積立シミュレーション設定
rate_percent = 5.1 # @param {type:"slider", min:1.0, max:20.0, step:0.1}
years = 30 # @param {type:"slider", min:10, max:50, step:1}
payment_man = 20 # @param {type:"slider", min:1, max:30, step:1}

plot_investment(rate_percent, years, payment_man)

colabで実行するとこのような画面が出てきます
スクリーンショット 2025-06-14 17.34.43

これを設定して実行すると
スクリーンショット 2025-06-14 17.34.53
積立のシミュレーション結果が反映されます。

colabで試せるので
コピペして試してみてください。



まとめ

金融庁のサイトのシミュレーターは
年利と月利を齟齬しているので
数十年後の金額が大きく乖離してくるようです。

単利と複利、月利と年利
この辺りの関係を正しく計算できるように
日々プログラムでシミュレーションすると安全です。

資産形成に取り組みたい方は
この機会にプログラミングも併せて
学んでみてはいかがでしょうか?


それでは。


今回は昔懐かしい
カイジの沼のシミュレーションです。

解説動画はこちら




カイジの沼のシミュレーション

沼とは

漫画『賭博破戒録カイジ』に登場する架空のパチンコ台

帝愛グループの運営するカジノに設置されている
パチンコの中でも特別な台で玉が1発あたり4000円
挑戦は最低300万円からという設定

その一方で、今までの挑戦者が消費した玉が
全てプールされる仕組みになっており
運良く大当たりを引ければ総取りで莫大な額の金を手にできるとか


仕掛け・役物

ここからは確率を計算するための
沼のギミックについてです。

1.釘の森

釘が密集した配置になっている「釘の森」
通常時の設定Cは1/100程度

2.スル

開閉する羽根のある電動役物「スルー」
釘の森を突破した玉を更にふるい落とす
確率は不明なので
ここでは通過率を1/5としておきます

3.3段クルーン

沼の本丸である「三段クルーン」
円形の皿に穴が3つ/4つ/5つ空いていて
当たりの穴に入ると次の段に進める

数学的にはクルーンに入ってさえしまえば
約1/60の確立で当たる


というような設定にはなっていますが
漫画の設定上ではこの台にはイカサマがあり
この確率では入らないようにはなっています。

しかし、今回は考慮しないで
イカサマなしだとどうなるかを見ていきたいと思います。



沼のシミュレーションコード


簡易な確率計算で大当たりしたかどうかだけ返す関数を用いて
シミュレーションしていきます。
最大回数は10万回(4億円相当)とします。

import numpy as np
import random
import pandas as pd
import matplotlib.pyplot as plt

# 確率でヒットするかどうかを判定する
def is_hit(n: int) -> bool:
    return random.randint(1, n) == 1

def numa_charenge():
  # 釘の森の通過
  mori = is_hit(100)
  if not mori:
    return False

  # スルーの通過
  throw = is_hit(5)
  if not throw:
    return False

  # クルーン1段目通過
  first = is_hit(3)
  if not first:
    return False

  # クルーン2段目通過
  second = is_hit(4)
  if not second:
    return False

  # クルーン3段目通過
  third = is_hit(5)
  if not third:
    return False
  else:
    return True


# 1000回施行して測ってみる
ball_price = 4000 # パチンコ玉の価格

def numa_result(num):
  all_result = []
  for a in range(1000):
    cumulative = 0
    for i in range(num):
      cumulative += ball_price
      is_atari = numa_charenge()
      if is_atari:
        all_result.append([i+1, cumulative,1])
        break
      if i == num-1:
        all_result.append([i+1, cumulative,0])
  return all_result

num = 100000 # 1回あたりの最大施行回数
all_result = numa_result(num)

df = pd.DataFrame(all_result,columns=["balls","cumulative","hit"])
df.head()
balls cumulative hit
0 30959 123836000 1
1 29604 118416000 1
...

plt.figure(figsize=(12, 6))
df["cumulative"].hist(bins=20)
plt.show()
download

おおよそこのような結果になります。


どれくらい注ぎ込めば当たるのか?


シミュレーション結果は使った額の累計を入れているので
それを用いれば、その金額以下の回数と
全体の回数から大体の確率を求められます。
times_300 = len(df[df["cumulative"]<3000000])
print(f"300万円で大当たりになる可能性 : {times_300/1000*100} % ")

なおこのシミュレーションの結果は
幾何分布に近似するため以下のようなコードで
回数あたりの大当たり確率の理論値が出せます。
# 幾何分布上の理論値
for i in range(10):
  p = 1 / 30000
  n = (i+1) * 10000
  prob = 1 - (1 - p)**n
  print(f"{i+1:02}万回で当たる確率 : {prob:.6f}")
01万回で当たる確率 : 0.283473
02万回で当たる確率 : 0.486589
03万回で当たる確率 : 0.632127
04万回で当たる確率 : 0.736409
05万回で当たる確率 : 0.811130
06万回で当たる確率 : 0.864669
07万回で当たる確率 : 0.903032
08万回で当たる確率 : 0.930520
09万回で当たる確率 : 0.950215
10万回で当たる確率 : 0.964328


シミュレーション結果もこの理論値に
近似してますね。


まとめ


沼の設定が玉が1発あたり4000円
3万発に1回の確率(100x5x3x4x5)
で当たるのだとしたら3億注ぎ込めば
9割くらいは大当たりするのではないか?

あくまでイカサマが無い事が前提ですが.....
イカサマさえなければ無尽蔵にお金を注ぎ込めるなら
勝てる気がしなくもない今日この頃でした。

それでは

このページのトップヘ