アークナイツガチャにおける確率分布

# 概要
<a href="https://qiita.com/suiton/items/3d8a58cd2d7fb97ab539">前回の記事</a>からの続きです。前回の記事では、お目当てを引くまでにどれぐらいの回数が必要かという観点で累積分布の比較を行いましたが、今回は確率分布を比較してみたいと思います。具体的には、<b>N回のガチャで最高レアリティを引ける合計枚数の確率分布</b>ですね。一般的なガチャとアークナイツのガチャでそれぞれ比較してみます。

細かい計算過程とかどうでもいいからアークナイツのガチャって結局どんな感じなの？ってことだけ知りたいって人は、最後のまとめだけ見てください。

今回も Statistics and Machine Learning Toolbox の関数をガンガン使っていきます。

# 一般的なガチャの確率分布
一般的なガチャは二項分布に従います。二項分布とはN回の試行中に何回成功するかという離散確率分布です。最高レアリティ排出率は、前回の記事にて計算したアークナイツガチャの実質的な排出率と同じ値としています。（なので中途半端な値です）

```matlab 
% パラメータ設定 (必要に応じて変更)
%  - 試行回数 N
%  - モンテカルロシミュレーションの試行回数 N_mc
%      [注意] N * N_mc の数だけループを回すため、設定値には注意
%  - 最高レアリティ排出率 ssr_prob

N = 5000
N_mc = 1000000
ssr_prob = 0.0289

% 二項分布の確率分布オブジェクトの作成
pd_gen = makedist('Binomial', 'N', N, 'p', ssr_prob);

% 近似した正規分布 (X軸の範囲を絞るために先に計算)
mu_gen      = N * ssr_prob;
sigma_gen   = sqrt(N * ssr_prob * (1-ssr_prob));
pd_gen_norm = makedist('Normal', 'mu', mu_gen, 'sigma', sigma_gen);

% ±3sigma までのX軸定義
x_axis = max(0, floor(mu_gen-3*sigma_gen)) : min(N, ceil(mu_gen+3*sigma_gen));

% 確率分布プロット
figure(1)
clf
hold on;
grid on;
plot(x_axis, pdf(pd_gen, x_axis));
title('二項分布の確率密度関数')
xlabel('成功回数 X')
xlim(minmax(x_axis))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112872/1c898d97-9743-ded7-9a75-519f152c1006.png)

二項分布は以下の条件の時に正規分布に近似できます。（<a href="https://jp.mathworks.com/help/stats/binomial-distribution.html">二項分布 - Mathworks</a> より引用）
> N が大きく、p が大きすぎも小さすぎもしない場合、パラメーター N と p をもつ二項分布は、平均 N*p と分散 N*p*(1–p) をもつ正規分布で近似できます。

パラメータ設定にて、適切な N と p を設定しているため、以下のプロットで正規分布に近似できていることがわかります。

```matlab
% 確率分布の比較プロット
figure(2)
clf
hold on;
grid on;
plot(x_axis, pdf(pd_gen, x_axis));
plot(x_axis, pdf(pd_gen_norm, x_axis));
title('確率分布の比較')
xlabel('成功回数 X')
legend('二項分布', '正規分布近似')
xlim(minmax(x_axis))
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112872/fc1d657f-11e2-5931-8009-f0699d53afef.png)

正規分布の特徴としては、以下があります。
* 最頻値、平均値、中央値が一致する
* 平均値を中心にして左右対称
* ±1σの範囲で68%, ±3σの範囲で99.7%

プロットの表示範囲は ±3σ としているため、この範囲に99.7%が入ります。つまり、同じ試行を行うと 1000 回中 997 回はこの範囲に収まるということになります。

# アークナイツガチャの確率分布
## モンテカルロシミュレーション
さて、次はアークナイツガチャの確率分布です。前回の記事でも書いたように、「徐々に最高レアリティが出る確率が上昇し、最終的には確定になる」というタイプのガチャです。
この方式では、二項分布の前提である「確率 p は常に一定である」を満たさないため、二項分布は使えません。そのため、前回同様にモンテカルロ法にて、確率分布を求めたいと思います。

```matlab
% 結果格納用配列
res_mc = zeros(1, N_mc);

% Monte-Carlo Sim.
for i = 1:N_mc
    local_cnt = 0;    % S6がN連続出現していない回数
    s6_cnt    = 0;    % S6を引いた回数
    
    % ピックアップが出るまで回す
    for j = 1:N
        local_cnt = local_cnt + 1;
        
        % 現在の local_cnt に応じたS6排出率
        if local_cnt > 50
            s6_rate = 0.02 + (local_cnt-50)*0.02;
        else
            s6_rate = 0.02;
        end
        
        % S6排出判定
        if rand < s6_rate
            s6_cnt = s6_cnt + 1;
            local_cnt = 0;
        end
    end
    
    % 結果配列へ格納
    res_mc(i) = s6_cnt;
end
```

## ヒストグラムの計算と確率分布プロット

hist_mc 変数には、N_mc 回の試行結果があります。これに対してヒストグラムを取り、各ビンのカウントを試行数で割ると、モンテカルロシミュレーションから求めた確率分布となります。

```matlab
% ヒストグラムの各ビンのエッジを取得
%   各ビンの中央 = min : 1 : max
%   各ビンのエッジ = min-0.5 : 1 : max+0.5
res_mc_minmax = minmax(res_mc);
mids  = res_mc_minmax(1):res_mc_minmax(2);
edges = (res_mc_minmax(1)-0.5):(res_mc_minmax(2)+0.5);

% ヒストグラムを計算し、モンテカルロ試行数 N_mc で割る
hist_mc = histcounts(res_mc, 'BinEdges', edges);
hist_mc = hist_mc / N_mc;

% 確率分布プロット
figure(3)
clf
hold on;
grid on;
plot(mids, hist_mc)
title('モンテカルロシミュレーションで求めた確率分布')
xlabel('成功回数 X')
xlim(minmax(x_axis))
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112872/fc8698da-4098-0dfd-9d4d-aa4394c58bd2.png)


# 一般的なガチャとアークナイツガチャの確率分布を比較
ここまでで、それぞれのガチャの確率分布を求めてきました。それでは、二つの確率分布グラフを重ねてみます。

```matlab
% ガチャの確率分布比較
figure(4)
clf
hold on;
grid on;
plot(x_axis, pdf(pd_gen, x_axis));
plot(mids, hist_mc)
title('ガチャの確率分布比較')
xlabel('成功回数 X')
xlim(minmax(x_axis))
str = sprintf('一般的なガチャ (SSR率=%.3G%%)', ssr_prob*100);
legend(str, 'アークナイツガチャ')
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112872/e3c2aed2-eabd-2dad-9442-e4aba0a18443.png)

一般的なガチャと比較すると、分散が明らかに小さいですね。<b>上振れも起きにくいですが、下振れも起きにくいのが特徴</b>といえるでしょう。

# アークナイツガチャの確率分布は正規分布に近似するか？
プロットを目視で見る限り、確率分布は正規分布に近似しそうだなという印象は受けます。ですが、本当に近似するのでしょうか？

実際に正規分布にフィッティングさせて確かめてみましょう。

```matlab
% 正規分布へのフィッティング
pd_ak_norm = fitdist(res_mc, 'Normal');

% プロット
figure(5)
clf
hold on;
grid on;
plot(mids, hist_mc, 'Color', '#cccccc', 'LineWidth', 3)
plot(mids, pdf(pd_ak_norm, mids));
title_str = sprintf('正規分布近似 (N=%d, N\\_mc=%d)', N, N_mc);
title(title_str)
xlabel('成功回数 X')
xlim(minmax(x_axis))
legend('モンテカルロ', '正規分布近似')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112872/137ae35b-cbee-00bf-e87b-51f3386f6ab0.png)

僅かにですが、以下の点でフィットしきれていないように見えます。N = 5000 だとかなりわかりにくいですが、N = 1000 にするとかなり明確です。

* 中央値がズレている
* 左側の確率密度が正規分布のほうが高く、右側が逆の傾向

正規分布のフィッティングにしっくりこなかったので、いくつか他の確率分布へのフィッティングを試行していたところ、綺麗にフィットする確率分布がありました。対数正規分布です。

```matlab
% 対数正規分布へのフィッティング
pd_ak_lognorm = fitdist(res_mc, 'Lognormal');

% プロット
figure(6)
clf
hold on;
grid on;
plot(mids, hist_mc, 'Color', '#cccccc', 'LineWidth', 3)
plot(mids, pdf(pd_ak_lognorm, mids));
title_str = sprintf('対数正規分布近似 (N=%d, N\\_mc=%d)', N, N_mc);
title(title_str)
xlabel('成功回数 X')
xlim(minmax(x_axis))
legend('モンテカルロ', '対数正規分布近似')
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/112872/c3a6380a-1193-6e6b-d008-12d274444821.png)

正規分布よりも綺麗にフィッティングしています。こちらのほうが近似としては適切にみえます。

## 対数正規分布
調べてみたところ、どうやら以下の特徴があるみたいです。

* Xを対数軸としたときの正規分布
* 小さい側が多く、大きい側が少ない
    * 所得の分布とかで使われるらしい
* 最頻値、平均値、中央値が異なり、平均値 > 中央値 > 最頻値 となる
    * N = 1000 のケースでは、以下となった
        * 平均値：28.6
        * 中央値：28.4
        * 最頻値：28.0

一般のガチャでは平均値以下のことを下振れと呼んだりしますが、アークナイツガチャでは平均値以下でも下振れではない区間があります。全体をみると50％の境界は中央値ですので、上の例でいうと28.4~28.6の区間は、平均値以下でも50％よりは上（＝上振れ）と言えるでしょう。

# まとめ

2回の記事に分けてアークナイツのガチャについて、モンテカルロシミュレーションによる確率分布導出で分析してきました。ここまでの結果より、アークナイツのガチャは以下のような特徴があるといえるでしょう。

*  ★6の実質的な確率
    *  約 2.89 %
*  ピックアップが引けるまでの確率
    *  恒常追加時なら、約140連で90%
    *  限定なら、約220連で90％
        *  120連以降でピックアップ1.04%の一般的なガチャと同じような確率に漸近
*  N連でSSRが引ける枚数の確率分布
    *  一般的なガチャに比べると上振れ/下振れが少ない傾向
    *  平均値以下になる確率が高い
        *  平均値より下側のほうが確率密度が高いため

# 参考

https://sci-fx.net/math-log-norm-dist/