Zenn
📦

食玩問題〜複数個パッケージ販売の場合〜

2024/09/18に公開
Haskell
統計学
tech

食玩問題とは

ある1機のガシャポンが全n種のキーホルダーを売っているとします。このとき何回ガシャを回せば全部コンプリートできるか？などを考える問題です。コンプリートできる確率を求めたり、あるコンプ確率以上を出すための購入数の期待値であったりという観点から分析されています。

背景

私は「KAMITSUBAKI STUDIO」というバーチャル系音楽レーベルのファンであり、そこのグッズ展開ではこのステッカーの様に、全18種類を内容ランダム3個入りパック(パック内重複なし)で販売するなどという販売方式が多く取られています。これを「n種類s個入りパック販売する場合の、パック購入数とn種コンプリートする確率」の問題と捉え解いてみたいと思いました。なお、この疑問が頭に浮かび上がった当日にzenn記事：食玩問題についてが公開されたため少々運命を感じました。素晴らしい記事であり、この私の記事のベースとなっていますのでこれも読まれることをお勧めします。

Webアプリとしての実装

このページ以降では確率の式の導出を行いました。その式を使用してこのURLに実際の確率計算を行うWebアプリを作りました。先にこれでイメージを掴むと以降の説明がわかりやすいかも知れません。

問題の設定

上記の状況の元、全n種類s個入りパック販売のものを、rパック購入してn種類の内k種類が揃う確率a_{r,k}を求めます。

式の導出

n=4種類s=2個入りパックを例とします。
まず以下のようなことが言えます。(r:パック購入数,k:揃っている種類数)

\begin{align*} a_{r=0,k=0}&=1.0\ \text{(なにも買わない場合、0枚が必ず揃っている。)}\newline a_{r=0,k>0}&=0.0\ \text{(なにも買わない場合、1枚も揃わない)}\newline a_{r=1,k=2}&=1.0\ \text{(1枚買う場合、2枚は必ず揃う)} \end{align*}

最終的に4枚を揃える経路は以下のパターンなどがあります(実際の物販では1パックずつ購入で揃えることは状況的に難しいのですが説明しやすい様に状況設定しています)。

  • 1パック購入で2種揃える。次の1パック購入で、一回目と重複しない2種を引当ててコンプリート
  • 1パック購入。次の購入で1回目と一つ重複し3種そろう。次の購入で最後1種を引き当てる。
  • 1パック購入。もう1回パック購入で全4枚揃う。でも予備にほしいのでもう1パック購入。

このことから、以下のような漸化式が成立します。(r:パック購入数,k:揃っている種類数)

\begin{align*} a_{r=2,k=4}&=[2枚持っている時に重複しない2枚を引く確率]a_{1,2} \\ a_{r=3,k=4}&=[3枚持っている時に重複しない1枚を引く確率]a_{2,3} \\ &+[4枚持っている時に重複する2枚を引く確率]a_{2,4} \\ a_{r=2,k=3}&=[2枚持っている時に重複しない1枚を引く確率]a_{1,2} \end{align*}

これをr,kに関して一般化すると以下が成立します。

\begin{align*} a_{r+1,k}&=[k-2枚持っている時に重複しない2枚を引く確率]a_{r,k-2} \newline &+[k-1枚持っている時に重複しない1枚を引く確率]a_{r,k-1} \newline &+[k枚持っている時に重複しない0枚(逆に言えば既に持っている2枚)を引く確率]a_{r,k} \\ a_{r+1,k}&=\sum^{2}_{i=0}[k-i枚持っている時に重複しないi枚を引く確率]a_{r,k-i} \end{align*}

ここで式内の[k-i枚持っている時に重複しないi枚を引く確率]について考えていきます。
まずk=3,i=1として[2枚持っている時に重複しない1枚を引く確率]について考えていきます。全4種類をそれぞれ、A,B,C,Dと名前をつけると、1パックに含まれる2枚の組み合わせは以下のC(4,2)=6パターンです(Cは組み合わせの総数のCです。zennで綺麗に掛けませんでしたので、、、)。なおABは「AとBが封入されている1つのパック」の表記です。

AB,AC,AD,BC,BD,CD

仮にABが1パック目だったとすると、以下を引くことが[2枚持っている時に重複しない1枚を引く]の全パターンになります。

AC,AD　(既に持っているAを含む),
BC,BD　(既に持っているBを含む)

この数はC(2,1)\cdot C(2,1)=4となります。先のCは既に持っているA,Bの内1枚を選ぶ場合の数であり、後のCは全4種類の内で持っていない2枚の内1枚を引く場合の数です。この意味を踏まえると、n=4,s=2,k=3,i=1としてC(k-i,s-i) \cdot C(n-(k-i),i)となることが分かります。

確認のため、k=4,i=2として[2枚持っている時に重複しない2枚を引く確率]を見てみましょう。1パック目がABとすると残りはCDのみですので場合の数は1となります。これはn=4,s=2,k=4,i=2としたC(k-i,s-i) \cdot C(n-(k-i),i)=C(2-2,2-2) \cdot C(4-(4-2),2)=C(0,0)\cdot C(2,2)と一致します。他にもn=5の場合なども試しましたが、同様に言えるということが確認できています。場合の数はこれで分かりましたので、各1パック購入のパターン数も合わせると[k-i枚持っている時に重複しないi枚を引く確率]は以下となります。

\frac{C(k-i,s-i) \cdot C(n-(k-i),i)}{C(n,s)}

結果の漸化式

これまでの結果を踏まえると以下となります。

a_{r+1,k}=\frac{1}{C(n,s)}\sum^{s}_{i=0}{C(k-i,s-i)C(n-(k-i),i)a_{r,k-i}}

プログラムでの再帰実行のために購入数rに関して変形した式が以下です。

a_{r,k}=\frac{1}{C(n,s)}\sum^{s}_{i=0}{C(k-i,s-i)C(n-(k-i),i)a_{r-1,k-i}}

これは以下のようなHaskellコードとなります。高速化のためメモ化もしてます。

Syokugan.hs
module Syokugan where

import Data.MemoTrie (memo2)
import Data.Ratio ((%))

factorial :: Integer -> Integer
factorial n = product [1..n]

combi :: Integer -> Integer -> Integer
combi n r = fn `div` (fr * fm)
  where
    fn = factorial n
    fr = factorial r
    fm = factorial $ n - r

solve :: Integer -> Integer -> Integer -> Integer -> Rational
solve nn ss rr kk = memoizedSolve rr kk
  where
    memoizedSolve = memo2 solve'
    rcombi a b = toRational $ combi a b
    solve' :: Integer -> Integer -> Rational
    solve' 0 0 = 1 % 1 
    solve' r 0 = 0 % 1 
    solve' 0 k = 0 % 1 
    solve' r k   
      | (r * ss) < k = 0 % 1
      | ss > k = 0 % 1
      | otherwise = result
          where
            c = sum $ map ft [0..ss]
            ft i = memoizedSolve (r-1) (k-i)
                   * rcombi (k-i) (ss-i)
                   * rcombi (nn - (k - i)) i
            result = c / (rcombi nn ss)

以下の様に実行出来ます。

import Syokugan (solve)

main = do
  let n = 11  -- All goods number.
  let s = 3  -- Set number.
  let r = 10 -- Buy try number.
  let ans = fromRational $ solve n s r n
  print ans

結果：

0.6056587736776045

この結果から、全11種類の3個入り1パックの商品を10パック購入した場合のコンプリート確率は61%程度であることが分かります。

モンテカルロ法との結果比較

モンテカルロ法はランダムシミュレーションとも言われる通り、現実のランダム要素の生成過程をランダム関数で模倣して大まかな確率や確率分布を算出する方法です。例えば、1~6の数値をランダム関数で生成すれば現実の6面サイコロを模倣したことになります。これを今回複数個パッケージ販売の例に適用して先述の数式と比べて確からしさを検証しました。例えば全n=4種類s=2個入りパック販売の場合、1パックの内容は以下の何れかです。これらからランダムに10個を復元抽出すればr=10パックを購入したこととなります。

AB,AC,AD,BC,BD,CD

このr=10パック購入を10000回試行し、各試行で全コンプしたかの統計を取れば、全4種類のコンプ確率であるa_{10,4}のある程度の精度を持った値が分かります。

以下のコードで確認しました。全種類数n=10,15,20、パック内個数s=1,2,3、購入数r=5,10,15,30の全(n,s,r)の組み合わせで全て誤差\pm 1%程度となっており、先の式の正しさがある程度保証出来たかと思います。

Simu.hs
module Simu where

import Control.Parallel.Strategies (Eval, runEval, rpar)
import Data.Set as S
import Data.List (subsequences, permutations)
import System.Random (randomRIO)
import Control.Monad (replicateM)
import Control.Parallel.Strategies (runEval, parMap, rpar)

combinations :: Int -> [a] -> [[a]]
combinations m xs = Prelude.filter ((== m) . length) (subsequences xs)

combinationsWithDup :: Int -> [a] -> [[a]]
combinationsWithDup 0 _  = [[]]
combinationsWithDup m xs = [y:ys | y <- xs, ys <- combinationsWithDup (m-1) xs]

pick :: Int -> Int -> Int -> [[[Int]]]
pick kind_num set_num pick_num =
  let
    k_all = [1..kind_num]
    combi = combinations set_num k_all
  in
    combinationsWithDup pick_num combi

matchall:: Int -> [[Int]] -> Bool
matchall num packs
  = num == (length . S.unions $ Prelude.map S.fromList packs)


pickNRandom :: Int -> [a] -> IO [a]
pickNRandom n xs = do
    indices <- replicateM n $ randomRIO (0, length xs - 1)
    return [xs !! i | i <- indices]

trueRate :: [Bool] -> Double
trueRate xs = trueCount / totalCount
  where
    trueCount = fromIntegral $ length (Prelude.filter id xs)
    totalCount = fromIntegral $ length xs


calcSimuProba :: Int -> Int -> Int -> Int -> IO Double 
calcSimuProba kind_num set_num buy_num proba_try_num  = do
  let all_combi = pick kind_num set_num buy_num
  let combi = combinations set_num [1..kind_num]
  ps <- replicateM proba_try_num
          $ fmap (matchall kind_num)
            $ pickNRandom buy_num combi
  return $ trueRate ps
compare.hs
import Simu (calcSimuProba)
import Syokugan (solve)
import Numeric (showFFloat)


main :: IO ()
main = do
  let proba_try_num = 10000

  let params = [
                (k, s, r) |
                  k<-[10, 15, 20],
                  s<-[1, 2, 3],
                  r<-[5, 10, 15, 30]
               ]
  let xs_sim = map (\(k, s, r) -> calcSimuProba k s r 10000) params
  let xs_calc = map (\(k, s, r) -> solve' k s r) params
                  where
                    solve' k s r = fromRational $ solve nn ss rr nn
                      where
                        nn = toInteger k
                        ss = toInteger s
                        rr = toInteger r
  xxs_sim <- sequence xs_sim
  -- print xxs_sim
  -- print $ xs_calc

  let diff = map abs $ zipWith (-) xxs_sim xs_calc
  --print $ diff
  mapM_ pput $ diff
    where pput x = putStrLn $ showFFloat (Just 6) x ""

おわりに

高校数学レベルの統計知識をイジって、趣味と結びついたある程度使える数式が立式出来たのは感動でした。本当は参考としたzenn記事：食玩問題についての様に行列形式にしてからの対角化、一般項の導出までやりたく挑戦しましたが、固有ベクトルの計算で心が折れました、、、数学に強かったらもう少し進めると思うので勉強を続けます。

Discussion