ランダムな推しグッズへの抵抗

これは ISer Advent Calendar 2022の12/11の記事として書かれました。

ISer Advent Calendar 2022 - Adventar

はじめに

3年間も推し活をしていると次のような問題に出会うことがとても多いです。

❓

今回発売されたグッズは全部で10種類あります。このうち特定の4種類が欲しいとします。グッズはすべてランダムで手に入るとしたとき、その4種類を手に入れるためには平均何個買うことになりますか？

アイドルの写真集などはランダムで1枚ポストカードが封入されていたりしますし、中身が見えない袋に入った状態でメンバーの人数分のグッズが展開されたりすることも日常茶飯事です。特にこの前の日向坂とココイチのコラボキャンペーンは地獄でした。詳しくはこちらの記事を参照してください。

ココイチキャンペーン第1弾前半戦

ココイチキャンペーン第1弾後半戦

この記事では、グッズを買うときに便利な確率質量関数を作ろうと思います。具体的には、グッズが $n$ 種類あって、欲しいグッズが特定の $m$ 種類のとき、 $x$ 回目で欲しいグッズが揃う確率 $G_{n,m}(x)$ を求めます。できるだけお金を使わないようにしたいので、グッズは一つずつ買って中身を確認するという戦略を採用することにします(実際は、生写真など5枚組だったりするものもあるので、一概にこの戦略が使えるわけではありません)。もし、たとえば5個一気に買うような場合は、5回目までに揃っている確率( $\sum_{k=1}^{5}G_{n,m}(k)$ )を計算すれば、5個買って揃う確率が出ると思います。

受験数学では場合の数・確率はそこまで得意じゃなかったので、間違いがあったら申し訳ないです。

例題を解いてみる

まずは上に挙げた例題を解いてみましょう。

4回で試行が終わる場合

これは簡単ですね。4回引いて、欲しい4種類が出るということなので、確率は $4!/10^4$ です。4回で欲しい4種が出るほど世の中は甘くありません。

5回で試行が終わる場合

まず5回目で初めて出るグッズを選びます( $_4C_1$ 通り)。そして最初の4個については、欲しいグッズだけで構成されている場合と、欲しくないグッズが含まれている場合があります。

欲しいグッズ(3種類)だけで構成されている場合は、1種類がダブっているので、何がダブるかで $_3C_1$ 通りあり、順列が $4!/2!$ 通りあるので、合計 $_3C_1\cdot4!/2!$ 通りあります。

欲しくないグッズが混ざっている場合は、欲しくないグッズの選び方が $10-4=6$ 通りあり、欲しいグッズ3つと欲しくないグッズ1つの順列が $4!$ あるので、合計 $6\cdot4!$ 通りあります。

したがって確率は、

\frac{_4C_1(_3C_1\cdot4!/2! + 6\cdot4!)}{10^5} = \frac{720}{10^5}

❌ Unsupported block (unsupported by Notion API)

となります。ここまではなんとか計算できます。720は6!ですね。では次は8!でしょうか。残念ながら違います。

6回で試行が終わる場合

6回目で初めて出るグッズは $_4C_1$ 通りです。そして最初の5個については、前回と同様に欲しいグッズだけで構成されている場合と、そうでない場合があります。が、そこはあまり問題ではありません。問題なのは、 $k$ 種類の欲しいグッズを重複ありで $n$ 個並べる場合の数を求める必要があることです。ただの重複順列ではなく、 $n$ 個のうち、欲しい $k$ 種類それぞれが少なくとも一つ以上ないとダメなので、少し厄介です。

例として、3種類の欲しいグッズを重複ありで5個並べる場合の数を求めてみましょう。求める場合の数は、1,2,3,4,5と書かれたボールがあり、それらをA,B,Cと書かれた箱に、空箱を作らないように入れる場合の数とも言えます。これは $3^5$ という重複順列から、空箱ができてしまうパターンを除けば良いです。

空箱が出るパターンは、2つの箱にボールが偏るパターン(どの2つの箱かで $_3C_2$ 通り、2つの箱に空箱を作らず5つのボールを入れるので $2^5-2$ 通り)と、1つの箱にボールが偏るパターン(3通り)があります。

よって合計、 $3^5-{}_3C_2(2^5-2)-3=150$ 通りあります。

このような場合の数は、第2種スターリング数に近いです。第2種スターリング数 $S(n,k)$ はWikipediaによると、

📖

番号付けされた

n

個の要素をグループ

k

個に分割する組み合わせの数

と定義されているので、グループに番号を振ってあげれば、上で扱った場合の数に一致します。つまり、上の例でいうと、 $3!S(5,3)=150$ ということです。

今回は、グループを区別する場合の数のほうをよく使うので、

S'(n,k) = k!S(n,k)

と定義しておきます。

これを使うと、6回で試行が終わる場合の確率が出せます。

\frac{_4C_1(S'(5,3)+(10-4) _5C_1\cdot S'(4,3)+(10-4)^2{}_5C_2\cdot S'(3,3))}{10^6}=\frac{13560}{10^6}

欲しくないグッズが含まれている場合は、どの欲しくないグッズを選ぶかとそれをどこに配置するかを決めた後、残ったスペースに欲しいグッズを上の要領で配置します。

7回以降で試行が終わる場合

6回目で試行が終わる場合と同様に計算できます。実は4,5回目で試行が終わる場合に関しても6回目と同じ議論で求められます。

以上の考察により、 $x$ 回グッズを買ったときに試行が終わる確率 $G_{10,4}(x)$ が求まります。

G_{10,4}(x) = \frac{4\sum_{i=0}^{k-4} {}_{k-1}C_i\cdot (10-4)^i \cdot S'(k-1-i,3)}{10^x}

ためしに以下のようなプログラムで期待値を求めてみましょう。

プログラム

import math

def Sprime(n, m):
    res = 0
    for k in range(1, m+1):
        if (m - k) % 2 == 0:
            res += math.comb(m, k) * k ** n
        else:
            res -= math.comb(m, k) * k ** n

    return res


def G(n, m, x):
    if n == 1:
        return 1
    if m == 1:
        return 1/n * (1-1/n) ** (x-1) 
    
    target = 0
    for i in range(x - m + 1):
        target += math.comb(x-1, i) * (n-m) ** i * Sprime(x-1-i, m-1) 

    target *= m
    res = target / n ** x
    return res

def E(n,m):
    res = 0
    for x in range(m, n * 10):
        res += x * G(n,m,x)

    return res


if __name__ == "__main__":
	print(E(10,4))

スターリング数はたぶん漸化式で計算したほうがいいんだろうけど、

$ python goods.py
20.820465970094187

だいたい21個です。想像よりかなり多いですね。

一般化

上のプログラムで少し出てきてしまったのですが、一応、一般の場合についてまとめておきます。

グッズが $n$ 種類あり、特定の $m$ 種類が欲しいときに、 $x$ 回目で試行が終わる確率 $G_{n,m}(x)$ は

G_{n,m}(x) = \left\{ \begin{array}{ll} 1 & (n = 1)\\ (1-\frac{1}{n})^{x-1}\frac{1}{n} & (m = 1)\\ \frac{m}{n^k}\sum_{i=0}^{k-m}{}_{k-1}C_i\cdot(n-m)^i\cdot S'(k-1-i,m-1) & (otherwise) \end{array} \right.

となります。 $n=1$ のときは、グッズが1種類しかないので、1回引けば当たります(2回以降も便宜的に確率は1とします)。 $m=1$ のときは、幾何分布の確率質量関数とほぼ同じです。

また第2種スターリング数の一般項は以下のとおりです。

S(n,k)=\frac{1}{k!}\sum_{m=1}^{k} (-1)^{k-m}{}_mC_k m^n

これは、

S'(n+1,k)= (n+1)^k - \sum_{m=1}^{n} {}_{n+1}C_m S'(m,k)

を用いて、数学的帰納法で証明できます(当初、スターリング数の存在を知らず、自力でスターリング数を見つけて、自力で一般項を導出、証明した私を褒めてください)。

できれば期待値や分散も他の確率分布のようにすっきりした形で書きたいのですが、

E_{n,m}(X) = \sum_x x\cdot G_{n,m}(x)

のような無限級数を求めることは僕にはできませんでした。スターリングの公式とかを使ったらいけるのかもしません。コンピュータで計算することにします。

確率分布

さて、私達はいまグッズを買う際に、欲しいグッズが手に入る確率を確認したいのでした。一目でだいたい何個ぐらい買うと全部揃うのかを確認できるようなグラフを作ってみます。

プログラムはこちら

スターリング数はたぶん漸化式で計算したほうがいいんだろうけど、めんどくさいので一般項で計算してます。多分結構大きい引数を与えるとオーバーフローします。そこまでグッズ数多いことはあんまりないので、これでいいのです。

import math
import random
import matplotlib.pyplot as plt
import numpy as np
import sys


def Sprime(n, m):
    res = 0
    for k in range(1, m+1):
        if (m - k) % 2 == 0:
            res += math.comb(m, k) * k ** n
        else:
            res -= math.comb(m, k) * k ** n

    return res 

def G(n, m, x):
    if n == 1:
        return 1
    if m == 1:
        return 1/n * (1-1/n) ** (x-1) 
    
    target = 0
    for i in range(x - m + 1):
        target += math.comb(x-1, i) * (n-m) ** i * Sprime(x-1-i, m-1) 

    target *= m
    res = target / n ** x
    return res
        

def E(n, m):
    res = 0
    for x in range(m, n * 10):
        res += x * G(n, m, x)

    return res

def V(n,m):
    E_X = 0
    E_Xsqr = 0

    for x in range(m, n * 10):
        E_X += x * G(n, m, x)


    for x in range(m, n * 10):
        E_Xsqr += x**2 * G(n, m, x)
    
    return E_Xsqr - E_X ** 2


if __name__ == "__main__":
    args = sys.argv

    if len(args) < 3:
        print("Two arguments are required")
        exit()

    
    n = int(args[1])
    m = int(args[2])

    limit = int(n) * 10
    xs = list(range(1,limit))
    
    ys = []
    cums = []
    cum = 0

    for x in xs:
        cums.append(cum)
        y = G(n,m,x)
        cum += y
        ys.append(y)
    
    fig, ax1 = plt.subplots()

    ax1.bar(xs, height=ys, width=1.0, label='pmf')
    ax1.minorticks_on()
    
    Ex = E(n,m) 
    Vx = V(n,m) 
    Sx = math.sqrt(Vx)

    ax1.vlines(x=Ex, ymin=0.0, ymax=np.max(ys), color='orange', linestyle='--', label='$E[X]$')  
    ax1.vlines(x=Ex - Sx, ymin=0.0, ymax=np.max(ys), color='orange', linestyle=':', label='$E[X] - \\sqrt{V[X]}$')  
    ax1.vlines(x=Ex + Sx, ymin=0.0, ymax=np.max(ys), color='orange', linestyle=':', label='$E[X] + \\sqrt{V[X]}$')  

    ax1.grid(axis='x', which="major", color="black", alpha=0.5)
    ax1.grid(axis='x', which="minor", color="gray", linestyle=":")


    ax2 = ax1.twinx()
    ax2.plot(xs,cums,color="red", label="cdf")
    ax2.minorticks_on()
    
    
    ax2.set_yticks([0.0, 0.25, 0.50, 0.75, 1.0])
    ax2.grid(which="major", color="black", alpha=0.5)
    ax2.grid(which="minor", color="gray", linestyle=":")
    
    
    h1, l1 = ax1.get_legend_handles_labels()
    h3, l3 = ax2.get_legend_handles_labels()
    ax2.legend(h1 + h3, l1 + l3, loc='lower right')

    plt.savefig("goods.png")

上のグラフは $n=10,m=2$ としたときの確率分布と累積確率分布です。他にも期待値や1シグマ範囲を表示しています。

$n=10,m=2$ としたのには理由があります。12月17日から始まる日向坂とココイチのコラボキャンペーンで、チキン4つセットを一箱買うと1枚ポストカードがもらえるというキャンペーンがあるのですが、そのポストカードが全部で10種類あり、欲しいポストカードが2種類なのです(きょんことひなちゃんです)。

そして私はいまそのチキンを5箱予約しているので、上のグラフから5個で揃う確率は15%しかないことがわかります。予想より低くて軽く絶望しています。

じゃあもうせめてきょんこだけでもいいので5個で当たる確率が知りたい、となったら、 $n=10,m=1$ とすればよいです。

おお！40%の確率で当たるぞ！少し希望が持てます。

では、12月20日に発売する金村美玖の写真集のランダムポストカードはどうでしょうか。全部で6種類あり、どうせみんな全種類ほしいと思うので、 $n=6,m=6$ として画像を生成してみます。

平均で14冊ほど買う必要がありそうですね。これは、各店舗ごとのポストカード特典と限定カバー特典を集めてたら多分このぐらいに到達すると思うので、ちょうどいいですね。

私は今、非常に金欠なので、1冊で十分です。

上のグラフたちを眺めていたら、形がガンマ分布に似ていることに気づきました。なので、ガンマ分布のフィッティングをした結果、かなりの精度でフィットすることがわかりましたが、ガンマ分布の引数と $n,m$ がどう関係しているかがいまいちわからないので、今後の課題とします(しない)。もしかしたら $G_{m,n}(x)$ がガンマ分布っぽく書けるのかもしれませんが、どう式変形するのか検討もつかないのでやめておきます。

最後に

常に金欠気味な私にとっては、どのくらいグッズを買うことになるのかを予測することは非常に重要です。今回、求めた確率分布を手軽に計算できる形で表現することはできませんでしたが、多少なりとも使える確率質量関数が作れたのではないかと思います。

今後もランダム性のあるグッズを大量に買うことになると思うので、その都度確率を見積もってから買う癖をつけようと思います。推しが出るかどうかが結果を見るまでわからないというスリルがなんだかんだ好きです。

参考文献

← Go home