仮説検定でパチスロの設定判別してみる

概要

基本的にデータというものは母集団すべてを認識するのは難しい．だからこそ母集団の平均(母平均)と分散(母分散)を知りたいときは推定するわけなのだが，ふと思ったことがある．

パチスロって母平均わかってるじゃん．それなら小役確率から今打ってる台の母平均推定できるじゃん．それなら設定判別に使えるんじゃね？オレ頭いい～！

なんて馬鹿なことを思いついたわけです．というわけで実際のデータも使いながら，Pythonで統計的に設定判別してみる．

概要
題材
確認したいこと
小役確率から推定される，台の設定であり得るものはどれか?
設定1だとして、この小役確率になるのってあり得るのか?
- ド・モアブル=ラプラスの定理を使う
おまけで設定1と設定5は3000Gで統計的に判別可能なのか?
最後に

題材

この前打ったニューゲッターマウスのユニメモで取れたデータをもとにやってみる．ユニメモは以下のようになっていた．まだリーチ目33個達成していないのでとりあえずここまで．RB中はサブロー推しで4回だった．

また設定別の各小役確率は以下の通り:

設定	オレンジA	オレンジB	スイカ	チェリー
1	1/20.5	1/20.4	1/59.9	1/14.0
2	1/20.3	1/20.2	1/59.0	1/13.8
5	1/19.6	1/19.5	1/56.5	1/13.5
6	1/18.6	1/18.9	1/54.9	1/13.3

※2022/7/29に下記を参考にしています 1geki.jp

実戦で得られた小役確率と設定別確率を見比べると，オレンジA以外は高設定っぽい挙動になっている．RBを引いた回数は少ないが，もともとRBに寄っていた台だったのでやめるころにはBB:RB=1:1くらいになっていた．
そしてBB中は奇数っぽい感じ．これはもしや5だったのでは…?となっているわけです．
実際にはオレンジAはそこまで良くないし，それ以外の強引きなんじゃないの?とも思えてしまうんですが，よく打ってる台や店ならまだ判断する材料はありますが今回はそうはいかない．
というわけで統計学の出番です．

確認したいこと

小役確率から推定される，台の設定であり得るものはどれか?
設定1だとして、この小役確率になるのってあり得るのか?
おまけで設定1と設定5は3000Gで統計的に判別可能なのか?

小役確率から推定される，台の設定であり得るものはどれか?

とりあえずオレンジBを例にとって考えてみます．

ここで使うのは95%信頼区間です．
信頼区間とは何か？それは母集団から取ってきた一部のデータを使って，母平均や母比率等がどれくらいなのかを示す範囲です．
また95%信頼区間の場合，この信頼区間を求める作業を100回行ったときに95回はその区間の中に母平均や母比率が含まれるということになります．

これを利用して，今回の2839G回した時のオレンジBの出現率1/17.9の時の二項分布を考え，それに対する95%信頼区間を求めてみます．
例えばこの信頼区間の外に設定1の出現確率があるようなら，95%の確率で設定1を否定できるということになります．

統計的に判断するために，そしてこの題材なので二項分布を導入して考えてみます．

二項分布は成功確率 $p$ の事象を $n$ 回行ったときに，成功が $k$ 回でる確率を表す分布です．式で書くと以下のような感じ．

\begin{equation} Bin(n, p) = {}_n C_k p^k (1-p )^{n-k} \end{equation}

コイントスの例だと $p=1/2$ となります．

では今回のパチスロの例だとどうだろうか．
成功をオレンジBが出た，失敗を出なかったと考えればよいわけです．そして設定別にオレンジBの出現確率は公開されています．
今回は2839GでオレンジBが159回出現していました．そして仮に設定1だとするとオレンジBの出現確率は $1/20.4$ です．

なので $n=2839$ ， $p=1/20.4$ ， $k=159$ となるわけですね．

計算してみます．実際に計算すると面倒なのでpythonでやってしまいましょう．今回はpython3.9.6を利用しています．
以下必要なライブラリを読み込んでいます．

import numpy as np
from scipy.stats import binom, norm #二項分布,正規分布
import pandas as pd
import matplotlib
from statsmodels.stats.proportion import proportions_ztest

pythonでやるとこんな感じ．

game = 2839
m = 159
p = m / game

rv = binom(game, p)
lower, upper = rv.interval(alpha=0.95)
print(f'lower={lower}, upper={upper}')
# lower=135.0, upper=183.0

大体135~183回くらいはオレンジBが出ると信頼区間が出ました．これを確率にしてみます．

print(f'lower=1/{game/lower}, upper=1/{game/upper}')
#lower=1/21.02962962962963, upper=1/15.513661202185792

というわけで大体1/21.03 ~ 1/15/51くらいの確率までは95%信頼区間に入っているようです．これはニューゲッターマウスの設定1から6まですべて入っていますね．1を否定はできないようです．
ゲーム数が3000Gもいってないからかもしれないので，仮に1日がんばれば回せるであろう10000Gの時に同じ確率だったとしたらどうかも見てみましょう．

game = 10000

rv = binom(game, p)
lower, upper = rv.interval(alpha=0.95)
print(f'lower={lower}, upper={upper}')
print(f'lower=1/{game/lower}, upper=1/{game/upper}')
#lower=515.0, upper=606.0
#lower=1/19.41747572815534, upper=1/16.501650165016503

これだと設定1どころか設定5まで否定してくれました．1日あればオレンジBだけでも十分に設定判別できそうです．
本音を言えば2~3000Gくらいで判別できるといいんですがね…

設定1だとして、この小役確率になるのってあり得るのか?

今回は2839GでオレンジBが159回出現していました．
ちなみに設定1だと期待値は約139回です．仮に設定1だとすると20回多く引いていることになります．

では実際に引けるものなのか?

まずは $n=2839$ ， $p=1/20.4$ ， $k=159$ をpython内で定義．

game = 2839
m = 159
p = 1/20.4

そして二項分布を定義し，オレンジBが0回から2839回出現する確率をすべて求めてみます．

xs = [i for i in range(game+1)]

rv = binom(game, p)
ys = rv.pmf(xs)

いったんグラフ化してみましょう．全範囲書くとわかりにくいので75回~200回くらいまでの確率で描画してみます．

pd.DataFrame({'game':xs, 'get':ys}).pipe(
    lambda df: df[df['game']>=75]
).pipe(
    lambda df: df[df['game']<=200]
).plot.line(x='game', y='get')

グラフを見ても大体の確率は見えなくもないですが，設定1を2389G回した時に159回オレンジBを取れる確率を求めてみましょう．

rv.pmf(159)
# 0.007838770272102068

0.78%くらいですね．ちなみに期待値である139回出現する確率は

rv.pmf(139)
# 0.03467343446793342

となり3.47%くらいになります．4倍以上の差がありますね．
とはいえどれくらい起こりえないのか数値が小さすぎてよくわかりません．

ド・モアブル=ラプラスの定理を使う

ド・モアブル=ラプラスの定理とは何か? 成功確率 $p$ ，試行回数を $n$ 回とします．
この時， $n$ が十分大きいときにこの2項分布は平均 $np$ ，分散 $np(1-p)$ の正規分布に近似できる，という内容のものです．

正規分布とは以下のような式であらわされる確率分布となっています．

\begin{equation} N(x; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right) \end{equation}

ここで $\mu$ は平均， $\sigma^2$ は分散を表す．

これを $\mu=0$ ， $\sigma^2=1$ 描画してみると以下のようになる．

xs = [i for i in np.arange(-4, 4.1, 0.1)]
norm_pdf = norm.pdf(x=xs, loc=0, scale=1)
pd.DataFrame({'game':xs, 'get':norm_pdf}).plot.line(x='game', y='get')

上で書いた二項分布のグラフによく似ていますね．
正規分布は見ての通り平均値を境に左右対称なグラフになり，釣り鐘のような形をしています．

また平均0，分散1にするようなデータの変換方法があり，特に標準化と呼ばれる方法があります．変換は以下のように行います．

\begin{equation} Z = \frac{X - \mu}{\sigma} \end{equation}

これは後ほど使います．

ここから上記の特性を利用した平均値の検定を行います．
具体的には，データから得られたオレンジBの出現確率 $\bar{X}$ と設定1の出現確率 $\mu$ が同等かどうかを検定します．

要は仮説検定を行うわけなのですが，仮説検定は帰無仮説 $H_0$ ，対立仮設 $H_1$ をそれぞれ立てます．
今回の例だと帰無仮説は $\bar{X} = \mu$ ，対立仮設は $\bar{X} \neq \mu$ となります．

仮説検定は帰無仮説がそれらしいかを検証し，そうでなければ帰無仮説を棄却し対立仮設を採用します．
その時の判断基準となる優位水準 $\alpha$ を定めるのですが，よく使われる $\alpha=0.05$ としておきます．

この優位水準のイメージですが，正規分布の外側5%の位置，今回の場合左右に2.5%ずつの領域を取るようなイメージになります．

グラフの縦軸は確率になるので，この赤い箇所は確率が低い場所であることがわかります．
つまりは低確率というのがどのあたりの領域のことなのかを定めているものとなります．
この領域に先ほどの標準化を行った後の数値が入っていれば，めったに起こらないことが起きたということになり，帰無仮説は棄却されます．逆に入っていなければ帰無仮説を否定する材料なしということになります．

では設定1でもあり得るのか，以下の手順で行います．

二項分布を正規分布で近似
得られた正規分布の平均値を標準化
標準化後の値が赤の領域にあるか確認

ではまず正規分布で近似します．やることは正規分布で近似したときの分散を計算します．

game = 2839 #ゲーム数
xbar = 159 #オレンジB出現回数
p = 1/20.4 #設定1のオレンジB出現確率
mu = 139 #設定1でのオレンジB出現回数の期待値

sig = (game * p * (1-p)) ** 0.5

次に標準化します．
$\bar{X}$ はオレンジBの出現回数になので159， $\mu$ はオレンジBの出現回数の期待値なので139，[tex \sigma]は上でも止めているので標準化後の値が計算できます．

z = (xbar - mu) / sig
print(z)
# 1.7385076331917815

結果約1.74ほどとなりました．

最後に標準化後の値が赤い領域にあるかを確認します．
赤の領域は片側で2.5%分あるので，2.5%，97.5%になる横軸の値をそれぞれ求めます．

lower = norm.ppf(q=0.025, loc=0, scale=1)
upper = norm.ppf(q=0.975, loc=0, scale=1)
print(f'lower={lower}, upper={upper}')
lower=-1.9599639845400545, upper=1.959963984540054

領域としては-1.96以下，もしくは1.96以上なら設定1を否定する材料になりますが，明らかに $-1.96<1.74<1.96$ ですね．というわけで設定1は否定できません．
先ほどのグラフに1.74のところに線を引いてみるとこんな感じ．

じゃあ10000Gまで回すとどうなるか．
10000Gまでで同じ確率でオレンジBを引き続けた場合だと以下のようにできる．

game = 10000
xbar = int(159 * 10000 / 2839)
p = 1/20.4
mu = 139

sig = (game * p * (1-p)) ** 0.5

z = (xbar - mu) / sig
print(z)
# 19.498957751686675

というわけでばっちり範囲外になったので設定1はほぼ否定できました．同じくグラフも書いてみるとこんな感じ．

ほぼ確率0ですねこれは．

おまけで設定1と設定5は3000Gで統計的に判別可能なのか?

普段設定判別する時は，本来であれば一つの小役だけでなく設定差のある小役は一通りみるし，他にもボーナス確率，お店の状況等色々と使います．
なので小役1つだけで設定判別するのはどうなのよという話になりますが，複数行うとまた処理が面倒なので(?)，またもやオレンジBで考えます．

とりあえず3000Gで設定1と設定5で差があると判断できるかを確認してみましょう．

この場合，2標本の比率の検定を行いますが、pythonだと以下のようにできる．

p1 = 1/20.4
p5 = 1/19.5
game = 3000

m1 = int(game * p1)
m5 = int(game * p5)

zvalue, pvalue = proportions_ztest([m1, m5], [game, game], alternative='smaller')
print(f'z-value={zvalue}, p-value={pvalue}')
# z-value=-0.3554093266554545, p-value=0.3611414813580965

やっているのは，出現率から求められるオレンジBの出現期待値が設定5のほうが大きいと判断できるか？という検定を行っています．結果を見るに，p値は0.05よりも大きいので3000G程度では大小関係はつけられない場合が多い，つまり設定1でも設定5と同じくらいには小役を引いてしまう可能性が高く，また逆も然りという感じでしょうか．

じゃあ10000Gだとどうでしょう.

p1 = 1/20.4
p5 = 1/19.5
game = 10000

m1 = int(game * p1)
m5 = int(game * p5)

zvalue, pvalue = proportions_ztest([m1, m5], [game, game], alternative='smaller')
print(f'z-value={zvalue}, p-value={pvalue}')
# z-value=-0.7130990801206974, p-value=0.23789222434418184

10000Gでも厳しそうですね．

最後に

最後の検定は検出力まで求めたかったけどよくわからなかったので断念．statsmodels.stats.power.NormalIndPowerが使えそうなんだけどね．

後，パチスロは分散を最大化するように，つまり荒くなるように抽選部分を作っていると聞いたことがあります．
それもそのはず，荒くないようにされるとAタイプで1/150でボーナス引ける台でも100G以内の連荘に期待ができなくなりますからね．これは結構厳しいです．
なので今回のように判定できるかは不安が残りますが，設定判別の一つの指標として使っていけたらなとは考えています．(ホールでの計算をどうするかは知らない)

数学的にどうなのさ?

大学時代にちょっと長く数学を勉強した人の雑記。数学のこと(主に統計)や趣味、メモなどが多くなります