楽しい！！！！！！！！！すごく勉強になったん！！！！！
大きく 2 つのやり方があって、高速ゼータ変換を用いた O(3ⁿ) から O(n2ⁿ) への高速化や、期待値に関する重要な考察をするなど色々やれるん。

問題概要

$n$ 人がテストを行う。毎秒ごとに $i$ 番目の人がテストを終える確率は $p_i/q_i$ で与えられる。全員がテストを終えるまでに要する時間の期待値を求めよ。

制約

$1 \le n \le 20$

解法 1: O(3ⁿ) な期待値 DP を O(n2ⁿ) に高速化する

まず、1 人の場合は有名問題で、毎ステップごとに確率 $p$ で終わる事象が終わるまでのステップ数の期待値は $\frac{1}{p}$ になる。これは期待値の方程式を立てればわかる:

$E = (1-p)(E+1) + p$

これを解いて、 $E = \frac{1}{p}$ が導出される。
それでは $m$ 人の場合 ( $S$ で表される集合の人たち) の場合の式を書いてみる。ここで簡単のために

$P_{S} = \prod_{i \in S} p_{i}$
$Q_{S} = \prod_{i \in S} (1 - p_{i})$

とおいておく。

$E_{S} = Q_{S}(E_{S} + 1) + \sum_{T \subsetneq S} ( Q_{T}P_{S-T}(E_{T} + 1))$

という感じの式になる。この式で $T$ は、 $S$ のうち 1 秒後も終わらずに残ったものたちを表す。この式を $E_{S}$ について解くと、

$E_{S} = \frac{1}{1 - Q_{S}} (Q_{S} + \sum_{T \subsetneq S} ( Q_{T}P_{S-T}(E_{T} + 1))$

となる。これはちょうど bit DP の枠組みで解けるように見える。状態量が $2^{N}$ な bit DP で、状態遷移が submask 全体に行き渡るタイプの bit DP なので計算量は $O(3^{N})$ になる。このままでは間に合わないので頑張って高速化する。

$E_{S}$ の計算をするのに毎回 $\sum_{T \subsetneq S}$ の形の足し算をするのが時間かかるので、ここをなんとかしたい。もしこれの中身が $S$ には依存せずに $T$ のみに依存する式だったら、

$t_{S} = \sum_{T \subsetneq S} ()$

みたいな補助変数を用意して、これを更新しながら計算すればよさそうである。DP を累積和を用いて高速化する話によく似ている。これができるようにするために、式変形をする。具体的には $P_{S-T}$ の部分をくくり出す。ここで

$R_{S} = \frac{Q_{S}}{P_{S}} = \prod_{i \in S} (\frac{1}{p_{S}} - 1)$

としておく。そうすると

$E_{S} = \frac{P_{S}}{1 - Q_{S}} (R_{S} + \sum_{T \subsetneq S} ( R_{T}(E_{T} + 1))$

という風に式変形できる。これでめでたく sum のところの中身が $T$ のみに依存する式になったので、アイディアとしては

$t_{S} = \sum_{T \subsetneq S} ( R_{T}(E_{T} + 1))$

とおくことで、これを毎回計算しながら頑張ればよい。 $t_{S}$ の更新は高速ゼータ変換の亜種 (naoya_t さんのブログの (2) バージョン) の出番である！

ここでいう高速ゼータ変換は

集合に対する関数 $f(S)$ があって

$g(S) = \sum_{T \subseteq S} f(T)$

を計算したい

というものを考える ( $T$ の向きが逆で $T \supseteq S$ の方を高速ゼータ変換と呼ぶことの方が多いかも)。

これは以下のようなシンプルなコードが知られているが

for (int i = 0; i < n; ++i)
  for (int bit = 0; bit < (1<<n); ++bit)
    if (bit & (1<<i))
      f[bit] += f[bit ^ (1<<i)];

これは配列再利用をしているため少し見えづらくなっている。結局以下のような DP をしている:

t[ i ][ bit ] := bit のうち、下から i 桁については bit の subbit 全体を動かし、それ以外の桁については bit と同じ値に固定したもの全体の総和

とすると

初期化: t[ 0 ][ bit ] = f(bit)
更新:

t[ i + 1 ][ bit ] =
t[ i ][ bit ]　( i が bit に含まれないとき)
t[ i ][ bit ] + t[ i ][ bit ^ (1<<i) ]　( i が bit に含まれるとき)

答え: t[ n ][ bit ]

なお、こうして t を二次元配列として定義してあげると、

for (int i = 0; i < n; ++i)
  for (int bit = 0; bit < (1<<n); ++bit)
    if (bit & (1<<i))
      f[bit] += f[bit ^ (1<<i)];

の i と bit のループの順番を変更しても大丈夫になる。これを利用して、dp[bit] を求めながら、t[n][bit] も求めていく。もう 1 つの注意点として、本来は t の初期化は

t[ 0 ][ bit ] = f(bit)

という感じで行うのだが、f(bit) の値が予めわからないので、後で

for (int i = 0; i <= n; ++i) t[i][bit] += R*(dp[bit] + 1);

という感じで足している。

#include <iostream>
#include <vector>
using namespace std;

class TestProctoring {
public:
    double expectedTime(vector <int> a, vector <int> b) {
        int n = (int)a.size();
        vector<double> p(n);
        for (int i = 0; i < n; ++i) p[i] = (double)(a[i]) / b[i];
        vector<double> dp(1<<n, 0.0);
        vector<vector<double> > t(n+1, vector<double>(1<<n, 0.0));
        for (int bit = 0; bit < (1<<n); ++bit) {
            double P = 1.0, Q = 1.0, R = 1.0;
            for (int i = 0; i < n; ++i) {
                t[i+1][bit] = t[i][bit];
                if (bit & (1<<i)) {
                    P *= p[i], Q *= (1.0 - p[i]), R *= (1.0/p[i] - 1);
                    t[i+1][bit] += t[i][bit ^ (1<<i)];
                }
            }
            if (bit) dp[bit] = P / (1 - Q) * (R + t[n][bit]);
            for (int i = 0; i <= n; ++i) t[i][bit] += R*(dp[bit] + 1);
        }
        return dp[(1<<n)-1];
    }
};

解法 2: 試行回数の期待値 = (k 回以上かかる確率) の和

てんぷらさんに教えてもらった解法です。kmjp さんの解法とも一致しています。今回はちょうど $k$ 回で終わる確率を $P_k$ として、期待値の定義から

$\sum_{k = 1}^{\infty} k P_{k}$

を計算したいです。しかし「ちょうど $k$ 回になる確率」というのはとても計算しづらいので、代わりに

$Q_{k} = k$ 回以上かかる確率

というのを代わりに考えます。このとき $P_{k} = Q_{k} - Q_{k+1}$ になるので、

$\sum_{k = 1}^{\infty} k P_{k}$

$= \sum_{k = 1}^{\infty} k (Q_{k} - Q_{k +1})$

$= (Q_1 - Q_2) + 2(Q_2 - Q_3) + 3(Q_3 - Q_4) + \dots$

$= Q_1 + Q_2 + Q_3 + \dots$

$= \sum_{k = 1}^{\infty} Q_k$

という感じになります (収束性についてはきちんと議論しないとですが...)。標語的に言えば

試行回数の期待値 = ( $k$ 回以上かかる確率) の $k$ についての無限和

ということになります。これは「ヒストグラムを書いて縦に集計するものを横に集計し直したもの」というようなものの見方ができます。結構汎用的なものの見方な気がします。AGC 023 C - Painting Machines でも出て来たとか。

さて、各 $k$ について $Q_k$ を計算します。 $k$ 回以上かかるというのはすなわち「 $k-1$ 秒たってもまだ誰かしらが終わってない確率」とういことになります。これは余事象をとると「 $k-1$ 秒で全員が終わる確率」です。 $r_{i} = 1 - \frac{a_{i}}{b_{i}}$ とおくと、 $i$ さんが $k$ 秒たっても終わらない確率は $r_{i}^{k-1}$ となることから、 $i$ さんが $k$ 秒で終わる確率は $1 - r_{i}^{k-1}$ になります。したがって、

$Q_{k} = 1 - (1 - r_{1}^{k-1})(1 - r_{2}^{k-1})\dots(1 - r_{n}^{k-1})$

になります。これを展開すると、次のように $2^{n} - 1$ 個の項が出て来ます。 $S$ を { $1, 2, \dots, n$ } の空でない部分集合を動くとして、

$Q_{k} = \sum_{S \in 2^{n}, S \neq \emptyset} (-1)^{|S|+1} \prod_{i \in S}r_{i}^{k-1}$

という感じになります。この式は包除原理で導いたものと考えることもできます (kmjp さんのブログでは包除原理)。

この項それぞれについて、 $k = 1, 2, \dots,$ についての無限級数の総和をとったものを求めればよいです。ここで

$R = \prod_{i \in S}r_{i}$

とおくと、無限級数の和は

$1 + R + R^{2} + R^{3} + \dots = \frac{1}{1 - R}$

になります。あとはこれを各 $S \in 2^{n}, S \neq \emptyset$ について足し引きすればよいです。

#include <iostream>
#include <vector>
using namespace std;

class TestProctoring {
public:
    double expectedTime(vector <int> a, vector <int> b) {
        double res = 0.0;
        int n = (int)a.size();
        for (int bit = 1; bit < (1<<n); ++bit) {
            double r = 1.0;
            int bitcount = 0;
            for (int i = 0; i < n; ++i) if (bit & (1<<i)) r *= 1.0 - (double)(a[i])/b[i], ++bitcount;
            double tmp = 1.0 / (1.0 - r);
            if (bitcount & 1) res += tmp;
            else res -= tmp;
        }
        return res;
    }
};

けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

TopCoder TCO 2018 Round 3B Medium - TestProctoring

問題概要

制約

解法 1: O(3ⁿ) な期待値 DP を O(n2ⁿ) に高速化する

解法 2: 試行回数の期待値 = (k 回以上かかる確率) の和

問題概要

制約

解法 1: O(3n) な期待値 DP を O(n2n) に高速化する

解法 2: 試行回数の期待値 = (k 回以上かかる確率) の和

解法 1: O(3ⁿ) な期待値 DP を O(n2ⁿ) に高速化する