高速ゼータ変換(FMT)とMin-Max反転による収集期待値の計算

ある$n$個の要素からなる集合において、各ステップで特定の部分集合が選ばれる確率分布が与えられているとします。全ての要素が少なくとも1回選ばれる状態になるまでに必要なステップ数の期待値を求めることを考えます。この種の問題は、確率論における「クーポン収集問題」の一般化と見なすことができます。

数学的背景と期待値の導出

まず、確率$p$で成功する試行を繰り返す場合、初めて成功するまでの期待回数が$\frac{1}{p}$であることを確認します。これは幾何分布の期待値として知られています。確率変数$X$を成功までの回数とすると、$X=k$となる確率は$(1-p)^{k-1}p$です。したがって、期待値$E[X]$は以下のように計算できます。

\[ E[X] = \sum_{k=1}^{\infty} k(1-p)^{k-1}p \]

ここで、幾何級数の和$\sum_{k=0}^{\infty} x^k = \frac{1}{1-x} (|x|<1)$を両辺微分して$\sum_{k=1}^{\infty} kx^{k-1} = \frac{1}{(1-x)^2}$を得ます。$x=1-p$を代入すると、

\[ E[X] = p \times \frac{1}{(1-(1-p))^2} = p \times \frac{1}{p^2} = \frac{1}{p} \]

となることが導けます。

Min-Max反転とアルゴリズム

本問題では「全ての要素が揃うまでの時間」を求めます。要素$i$が初めて選ばれる時間を$T_i$とすると、求める期待値は$\max(T_1, T_2, \dots, T_n)$の期待値です。ここで、Min-Max反転(Max-Minの包除原理)を適用します。この変換により、最大値の期待値は、特定の部分集合に関する最小値の期待値の線形和として表現できます。

集合$S$に対し、その要素のうち「いずれか1つでも選ばれる」までの期待時間を$E[\min_{i \in S} T_i]$とします。先ほどの導出より、これは1ステップで集合$S$のいずれかの要素が選ばれる確率$p_S$の逆数$\frac{1}{p_S}$に等しいです。したがって、全要素が揃うまでの期待値$E_{all}$は、空でない全ての部分集合$S$について以下の式で計算できます。

\[ E_{all} = \sum_{S \neq \emptyset} (-1)^{|S|+1} \frac{1}{p_S} \]

ここで、$p_S = 1 - q_S$であり、$q_S$は「集合$S$の要素が一つも選ばれない」確率です。これは「補集合$\bar{S}$の要素のみからなる部分集合が選ばれる」確率と等しくなります。

各部分集合が選ばれる生の確率が与えられている場合、任意の集合$U$について、「$U$の部分集合が選ばれる確率の総和」を高速に計算する必要があります。これには高速ゼータ変換(FMT)または高速モービウス変換、競技プログラミングではしばしばSOS DP(Sum Over Subsets Dynamic Programming)と呼ばれる手法を用います。これは、各ビットに対して操作を行うことで、$O(n2^n)$で全ての部分集合の累積和を計算するアルゴリズムです。

実装と例題

典型的な例として、各要素が独立にビットごとに選ばれる確率が与えられる問題(Luogu P3175 [HAOI2015] 按位或など)を考えます。以下のC++コードは、FMTを用いて各部分集合の確率を計算し、Min-Max反転によって期待値を求める実装例です。コードの可読性を高めるため、変数名やループ構成を変更しています。


#include <iostream>
#include <vector>
#include <cmath>
#include <bitset>
#include <iomanip>

using namespace std;

int main() {
    // 高速化のための標準入出力
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int num_elements;
    cin >> num_elements;

    const int total_subsets = 1 << num_elements;
    vector<double> probability(total_subsets);
    vector<int> bit_count(total_subsets);

    // 各部分集合が選ばれる確率を入力
    for (int mask = 0; mask < total_subsets; ++mask) {
        cin >> probability[mask];
        // 集合のサイズ(ビット数)を事前計算(Dinicのアプローチなど)
        bit_count[mask] = bit_count[mask >> 1] + (mask & 1);
    }

    // 高速ゼータ変換 (FMT)
    // prob[mask] に maskの部分集合である全ての集合の確率の和を格納する
    for (int bit = 0; bit < num_elements; ++bit) {
        for (int mask = 0; mask < total_subsets; ++mask) {
            if (mask & (1 << bit)) {
                probability[mask] += probability[mask ^ (1 << bit)];
            }
        }
    }

    double expected_time = 0.0;
    const double EPS = 1e-10;

    // Min-Max反転を用いた期待値の計算
    for (int mask = 1; mask < total_subsets; ++mask) {
        // maskに含まれない要素のみからなる部分集合が選ばれる確率を求める
        int complement_mask = (~mask) & (total_subsets - 1);
        double prob_none = probability[complement_mask];
        
        // 全く選ばれない確率が1に近い場合、収束しない(無限)
        if (1.0 - prob_none < EPS) {
            cout << "INF" << endl;
            return 0;
        }

        // 集合maskのいずれかが選ばれる確率
        double prob_any = 1.0 - prob_none;
        
        // 符号反転:要素数が奇数なら加算、偶数なら減算
        double term = 1.0 / prob_any;
        if (bit_count[mask] % 2 == 1) {
            expected_time += term;
        } else {
            expected_time -= term;
        }
    }

    cout << fixed << setprecision(10) << expected_time << endl;

    return 0;
}

7月21日 00:57 投稿