FGOのスター集中の確率・期待値の計算

#注意
この記事はFGOことFate/Grand Orderというソーシャルゲームのあるアルゴリズムについての記事です。

#はじめに
FGOではクリティカルが戦術の大きな部分を占めますが、クリティカルスターの集中の計算はちょっと面倒です。

> [サーヴァント/隠しステータス - Fate/Grand Order @wiki 【FGO】- アットウィキ](https://w.atwiki.jp/f_go/pages/304.html)

仕様は

* サーヴァント毎に集中度（SW）が設定されている。(10～200)
* 集中度アップ・ダウンの効果によって増減する。
* サーヴァント3騎から5枚のカードをランダムに配布する。
* それぞれのカードについて、ランダムに集中度にボーナス(+0x2, +20x2, +50x1)が加算される。
* ボーナス加算後の集中度に応じて保有しているクリティカルスターを振り分ける。
* 10個振り分けられたカードにはそれ以上振り分けられない。

#計算例
以下のようなサーヴァントで考えます。

* SW
 * サーヴァント1：200
 * サーヴァント2：50
 * サーヴァント3：10

サーヴァント1から2枚、2から2枚、3から1枚配布され、ボーナスを+0, +0, +50, +20, +20で考えると…

* いずれのカードもスター10個未満の時

| 配布カード | サーヴァント1<br>カードA | サーヴァント1<br>カードB | サーヴァント<br>カードC | サーヴァント2<br>カードD | サーヴァント3<br>カードE |
|:-:|:-:|:-:|:-:|:-:|:-:|
| 集中度 | 200 | 200 | 50 | 50 | 10 |
| ボーナス | +0  | +0 | +50 | +20 | +20 |
| 合計 | 200 | 200 | 100 | 70 | 30 |
| 比率(%) | 33.3 | 33.3 | 16.7 | 11.7 | 5 |

* カードAがスター10個になった時

| 配布カード | サーヴァント1<br>カードA | サーヴァント1<br>カードB | サーヴァント<br>カードC | サーヴァント2<br>カードD | サーヴァント3<br>カードE |
|:-:|:-:|:-:|:-:|:-:|:-:|
| 集中度 | - | 200 | 50 | 50 | 10 |
| ボーナス | -  | +0 | +50 | +20 | +20 |
| 合計 | - | 200 | 100 | 70 | 30 |
| 比率(%) | - | 50 | 25 | 17.5 | 7.5 |

* カードAとBがスター10個になった時

| 配布カード | サーヴァント1<br>カードA | サーヴァント1<br>カードB | サーヴァント<br>カードC | サーヴァント2<br>カードD | サーヴァント3<br>カードE |
|:-:|:-:|:-:|:-:|:-:|:-:|
| 集中度 | - | - | 50 | 50 | 10 |
| ボーナス | -  | - | +50 | +20 | +20 |
| 合計 | - | - | 100 | 70 | 30 |
| 比率(%) | - | - | 50 | 35 | 15 |

………と、このように確率が変化していきます。

#計算
ボーナスは振り分けの過程で変化しないので無視して、一般化して考えます。
カード$n$に$y$個のスターが振り分けられる確率は、

* カード$n$に$y-1$個スターがあって、カード$n$に1個スターが振り分けられる確率
* カード$n$に$y$個スターがあって、カード$n$以外に1個スターが振り分けられる確率

の合計です。
これをPythonプログラムで実装します。

#実装
```python:fgostarcalc.py

import math
#各カードのSW（ボーナス込み）：カード1～カード5
SW = [200, 200, 100, 70, 30]

#各カードのスター0～10個の状態の確率を保持
Probabilities = [None for _ in range(11**6)]
Probabilities[0] = 1

#任意の桁の数字の取得(各カードのスターの個数)
def digitnum(decim, base, ind):
    return(math.floor((decim % (base**ind))/(base**(ind-1))))

#各桁の数字の合計(スターの合計)
def sumdigit(decim, base):
    if decim == 0:
        return 0
    else:
        return(sum([digitnum(decim, base, i) for i in range(1, 2+math.ceil(math.log(decim, base)))]))

#有効性の検証
def is_valid(decim, base, totalstar, card, holdstar):
    if digitnum(decim, base, card) == holdstar:
        if sumdigit(decim, base) == totalstar:
            return True
    return False
        
#各状態の確率計算の再帰関数
def calc_starprob(decim, base, totalstar, card, holdstar):
    #print(decim, base, totalstar, card, holdstar)
    starweights = [0 if digitnum(decim, base, x+1)
                   == base-1 else SW[x] for x in range(5)]
    starprob = [SW[x]/(sum(starweights)-starweights[x]+SW[x])
                for x in range(5)]
    starweights[card] = SW[card]
    if decim < 0 or totalstar < 0 or holdstar < 0:
        print("Error: ", decim, base, totalstar, card, holdstar)
        raise
    if Probabilities[decim] != None:
        return(Probabilities[decim])
    else:
        tmp = [0]
        for x in range(5):
            if digitnum(decim, base, x+1) == 0:
                continue
            else:
                if x+1 == card and holdstar > 0:
                    tmp += [calc_starprob(decim-base**x, base,
                                          totalstar-1, card, holdstar-1)*starprob[x]]
                elif x+1 != card:
                    tmp += [calc_starprob(decim-base**x, base,
                                          totalstar-1, card, holdstar)*starprob[x]]
        Probabilities[decim] = sum(tmp)
        return(Probabilities[decim])

#totalstar個スターがある時のcardにholdstar個のスターが集まる確率
def calc_probability(totalstar,card,holdstar):
    return(sum([calc_starprob(x, 11, totalstar, card, holdstar)
                 for x in range(11**5) if is_valid(x, 11, totalstar, card, holdstar)]))

#totalstar個スターがある時のcardにholdstar個以上のスターが集まる確率
def calc_moreprobability(totalstar,card,holdstar):
    tmp=0
    for star in range(holdstar, 11):
        tmp+=sum([calc_starprob(x, 11, totalstar, card, star)
                for x in range(11**5) if is_valid(x, 11, totalstar, card, star)])
    return(tmp)

#totalstar個スターがある時のcardにholdstar個以下のスターが集まる確率
def calc_lessprobability(totalstar,card,holdstar):
    tmp=0
    for star in range(0,holdstar+1):
        tmp+=sum([calc_starprob(x, 11, totalstar, card, star)
                for x in range(11**5) if is_valid(x, 11, totalstar, card, star)])
    return(tmp)

#totalstar個スターがある時のcardに集まるスターの個数の期待値
def calc_expectation(totalstar,card):
    tmp=0
    for star in range(0, 11):
        tmp+=calc_probability(totalstar, card, star)*star
    return(tmp)

#cardに集まるスターの個数の期待値がstar個となるのに必要なスターの個数(整数値)
def calc_totalstar_expectation(card,star):
    if star>10:
        star=10
    for totalstar in range(50):
        if calc_expectation(totalstar,card)>star:
            return(totalstar)
```

#テスト

```python

#サンプル(25個スターがある時に8個がカード1に集まる確率)
print(calc_probability(25, 1, 8))
#0.1693947010897705

#サンプル(20個スターがある時にカード1に集まるスターの個数の期待値)
print(calc_expectation(20, 1))
#6.862519232471602

#サンプル(カード1に期待値4個以上集めるのに必要なスターの個数)
print(calc_totalstar_expectation(1,4))
#12
```

#おわりに
あまり早い実装とは言えませんが、一応計算できました。
確率、期待値、必要な個数なども計算できます。
ランダムなボーナスも$\frac{5!}{2!2!}=30$通りなのでこれを応用すればすべてのパターンを計算することが出来ると思います。