ヒューリスティックコンテストでベイス推定に入門しよう

ヒューリスティックコンテストでは、一部の数値が与えられないタイプの問題が出題されることがあります。多くの場合はインタラクティブ形式の問題で、少しずつ与えられる情報で数値を予測していくことになります。
(例：AHC003, HTTF2022予選, HTTF2023本選, AHC018)

こうした推定タイプの問題では、ベイズ推定が強力な手法となることがあります。
本記事はベイズ推定を使って問題を解いてみたい方の第一歩となることを目的としています(C++での実装例も載せました)。題材はAHC003としました。問題の概要は以下です。

• 30×30のグリッド(上下左右が辺で結ばれている)があり、各辺の長さは与えられない。
• 始点と終点が与えられるので経路を出力すると、その経路の長さ(±10%の誤差を含む)が返ってくる。これを1000回繰り返す。
• 出力した経路が最短経路に近いほど得点が高くなる。回数を重ねるごとに得点の比重が上がる(つまり後半が大事)。
• (入力で与えられない)真の辺の長さは以下の二つのパターンがある。
  • 縦横60本の同一直線上の辺は大体同じ長さ ($M=1$)
  • 各直線がランダムな位置で二分割されており、それぞれの側の辺は大体同じ長さ ($M=2$)

ということで、辺の長さを正確に推測できればダイクストラ法によって良い経路を導けそうです。

この記事は優勝者のeivourさんによる解説記事を参考にしています。
qiita.com
今回は行列演算の高速化やMCMCといった高度な内容は扱いませんので、線形代数と統計学の基礎知識があれば読めるのではないかと思います。より高い点数を目指したい方はeivourさんの記事をぜひご覧ください。
(筆者も勉強したばかりの内容も含んでいますので、厳密さに欠ける部分もあるかもしれません。不備がありましたらご指摘いただけますと幸いです。)

実際に問題に取り組む前にベイズ推定について説明します。既に理解している方は飛ばしてください。

ベイズ推定

ベイズ推定の基本となるベイズの定理について調べると、

$$P(A|B) = \frac{P(B|A)P(A)}{P(B)}$$

という数式を見ると思います。この式の意味は既に多くの解説があるのでここでは詳しく説明せず、早速具体例を見ていきましょう。
いきなりAHC003に挑戦しても良いのですが、その前に一つ簡単な例で感覚をつかむことにします。

• 表が出る確率が未知のコインがある。このコインを $3$ 回投げたところ全て表だった。表が出る確率を推定しなさい。

普通は表が出た回数の割合としたいですが、そうすると表が出る確率 $p=1$ となってしまい、さすがにおかしい気がしませんか？
でもなぜおかしいと感じるのでしょうか。それはおそらく、現実のコインは普通は表と裏が出る確率にそれほど大きな差はないという前提知識があり、たまたま偏っただけだと考えるからだと思います。この感覚を数学的に表現できるのがベイズ推定です。
表が出る確率 $p$ をピンポイントで推定しても良いのですが、確率分布とすることもできます。そうすると推定の自信の度合いを分散という数学的な値で表現できて便利です。

では、表が出る確率(分布)を $P(A)$ としましょう。これが $B$ というコイントスの結果を反映し、条件付き確率 $P(A|B)$ に変化したと考えます。このように観測結果をもとに確率(分布)を更新する、というのがベイズ推定の基本的な考え方です。$P(A)$ は観測前の確率分布なので「事前分布」、$P(A|B)$ は観測後なので「事後分布」と呼びます。
ベイズの定理の式を眺めると、事後分布 $P(A|B)$ を求めるには $P(B|A)$ と $P(B)$ がわかれば良さそうです。ただ $P(B)$ は $A$ とは無関係(すなわちこのコインとは無関係)な定数*1なので後回しにして $P(B|A)$ を考えます。これは $A$ という条件のもとでの $B$ の確率、すなわち「表が出る確率が $p$ とわかっている場合に $B$ というコイントスの結果が得られる確率」のことですから、

$$P(B|A) = {}_3 \mathrm{C}_3\ p^3 (1-p)^0$$

です。一般にコインを $n$ 回投げて $k$ 回表が出たのならば、

$$P(B|A) = {}_n \mathrm{C}_k \ p^k (1-p)^{n-k}$$

となります。
今まで $P(A)$ の具体的な形は考えていませんでしたが、$P(B|A)$ と同じ形、すなわち

$$P(A) = Cp^{\alpha} (1-p)^{\beta}$$

のような形($C$ は定数で、この形の確率分布をベータ分布と呼びます)をしていれば、

$$P(A)P(B|A) = C'p^{\alpha + k} (1-p)^{\beta + n - k}$$

といったように $P(B|A)$ をかけた後も同じ形になりますね。事後分布が事前分布と同じ形の確率分布になるような事前分布のことを特に「共役事前分布」と呼びます。計算がとても楽になるので、可能なら共役事前分布を選びたいです。
なぜ計算が楽なのかと言うと、定数 $C'$ が簡単に求まるからです。事後分布もベータ分布だとわかっているなら、途中の定数の計算は省略して最後に求めれば良いです*2。つまり、共役事前分布ならば後回しにしていた $P(B)$ の計算は実は不要でした。

これで事後分布が求められるようになりました。最後に事前分布における $\alpha, \beta$ の値を具体的に設定してみましょう。実はこれは特に決まりはなく、設定する人の「感覚」が反映されます。「普通のコインなら、平均は $0.5$ で多少のばらつきがありそう」みたいな感覚です。例えば $\alpha = 50, \beta = 50$ とすると、平均が $0.5$、標準偏差が約 $0.050$ のベータ分布になります。
コインを3回投げて全て表だったという結果を反映すれば、$\alpha = 53, \beta = 50$ と更新されます。これは平均が約 $0.515$、標準偏差が約 $0.049$ のベータ分布です。たった3回の試行だったので平均はそれほど大きく動きませんでした(このコインは少し表が出やすいだろう程度でしょうか)。試しにコイントスの結果($n, k$ の値)や、事前分布の設定($\alpha, \beta$ の値)を色々変えて結果を確認してみてください。

AHC003に挑戦

それでは実際の問題に取り組んでみましょう。この記事では簡単のため $M=1$ とします(すなわち同一直線上の辺は全て大体同じ長さとします)。上のベイズ推定の項目での議論から、

• 推定したい値の確率分布を設定する
• 観測結果をもとに確率分布を更新する

という二つの作業が必要だとわかると思います。逆に言えばこの二つが適切にできれば高い性能を期待できます。

確率分布の設定

まず最初の確率分布を設定しましょう。どんな確率分布を選ぶべきかは問題設定によりますし、センスが問われるところでもあって難しいのですが、今回は正規分布を用いることにします*3。というのは、返ってくる経路の長さの誤差も正規分布に従うとしてしまえば、正規分布は共役事前分布になる(事前分布が正規分布なら事後分布も正規分布になる)から嬉しいわけです。
推定したい辺は全部で $1740$ 本ありますから、多変量正規分布になります。つまりそれぞれの辺の分散に加えて共分散も必要です。ただ変数が $1740$ 個もあると大変なので、とりあえずそれぞれの直線だけを考えて $60$ 変数としてしまいましょう。こうするとそれぞれの変数は独立になりますから、共分散は全て $0$ です。
辺の長さ $X$ は平均(ベクトル) $\boldsymbol{\mu}$、分散共分散行列 $V$ の正規分布に従うとします。各辺の平均は $5000$ としておくと分散 $v$ は、

$$v = \frac{1}{8000-2D} \int_{1000+D}^{9000-D} (x-5000)^2dx = \frac{1}{3} (4000-D)^2$$

と計算できます。$D$ の値がわからないですが、例えば {$250, 750, 1250, 1750$} の $4$ 通りを設定して全て計算しておき、後でどれが近そうかを推測することにしましょう。これで $\mu, V$ を具体的に設定できました(今回は共分散が全て $0$ ですから $V$ は対角行列です)。

事後分布の導出

次に、得られた観測結果から事後分布を求める方法を考えます。出力した経路においてその辺を通る回数*4を要素としたベクトル $\boldsymbol{c}$ を考えると、経路の長さの期待値 は $c^T \mu$ とベクトルの内積で表せます($c$ や $\mu$ は列ベクトルで、転置の記号 $T$ を用いて行ベクトルとしています)。そして、返ってきた経路の長さ $y$ は誤差項 $\varepsilon$ が加わった

$$y = c^T X + \varepsilon$$

という確率変数であるとします。今回の誤差は±10%でしたから、誤差項 $\varepsilon$ は $X$ とは独立で、平均 $0$、分散 $R=(0.1y)^2$ (標準偏差 $0.1y$) の正規分布に従うとしてみます。
この時の事後分布の導出は、こちらの記事(分散のみですがWoodburyの公式*5についても言及があります)やこちらの記事(一般の場合について導出を追えます)が大変参考になると思います。
ここでは結論だけを書きます。$y$ の分散を $v_y = c^T Vc + R$ とおくと、事後分布の平均 $\mu '$ および分散共分散行列 $V'$ は、

$$\mu ' = \mu + Vc\ v_{y}^{-1} (y-c^T \mu)$$
$$V' = V - Vc\ v_{y}^{-1}c^T V$$

となります。すなわち $X$ の事後分布は平均 $\mu '$、分散共分散行列 $V'$ の正規分布です。

$D$ の値はどれが適切なのか？

ここまでは勝手に $D$ の値を $4$ 通り設定して全てについて計算を進めていましたが、その計算結果からどの値が適切かを推測するにはどうしたら良いでしょうか。
$X$ を既知とした時に実際の観測結果 $y$ が得られる確率を比較すれば良さそうです(その確率があまりにも低いなら、そもそも仮定した $D$ の値が不適切だったと言えそうです)。この確率とは何かと言うと、「パラメータがわかっている時にその観測結果が得られる条件付き確率」ですから、上のベイズ推定の議論で言う $P(B|A)$ に他なりません。これを統計学の用語で「尤度」と呼びます。
確率変数 $y$ は平均 $c^T \mu$、分散 $v_y$ の正規分布に従うので、尤度は

$$p(y|X) = \frac{1}{\sqrt{2\pi v_y}} \exp(-\frac{(y-c^T \mu)^2}{2v_y})$$

です(意味合いは異なりますが、形は確率密度関数と同じです)。全ての観測結果の同時確率について尤度を考えたいので、これまでの尤度を全てかけ合わせます。しかしかけ算を繰り返すと値がほとんど $0$ になってしまって扱いにくい(例えば $e^{-100}$ とかになったりします)ので、実装する際には自然対数を取って足し合わせることにしましょう。そして、全ての $D$ の値について尤度を計算してその比で「もっともらしさ」を判断することにします。

経路の選択

ここまでで、$4$ 種類のモデルについて尤度と事後分布を求められるようになりました。あとはダイクストラ法で経路を決定したいのですが、そのためには各辺の長さを固定する必要があります。
モデルごとに考えます。平均ベクトル $\mu$ をそのまま長さとしても良いかもしれませんが、せっかく分散も計算してあるのでそれも活用したいです。分散が大きい場合はその辺についての推定が弱いということですから、積極的に経路に組み込んで情報を得たくなります。そこで、分散が大きいほど辺の長さが小さくなるようにしてみたくなります。例えば、平均から標準偏差の定数倍を引いた

$$w = \mu - \gamma \sqrt{v_x}$$

などが考えられそうです(分散 $v_x$ は $V$ の対角成分です)。分散の影響をどの程度反映するかは調整が可能なので、実際に色々な計算式を試してみてください。モデルごとに辺の長さを確定させたら、尤度の比で重みをつけて最終的な辺の長さとします。

C++での実装例

ここで示したのはあくまで例であり、機能としても最低限のものです。ここまでの議論が理解できた方はぜひ自力で実装してみてください。
前提として行列演算のライブラリが必要です。C++の場合はEigenという便利なライブラリがありますが、コンテストで使うために自分で実装してみます。

struct VectorXd {
    std::vector<double> vec;
    bool transposed = false;

    VectorXd() {}
    VectorXd(size_t __new_size) {
        vec.resize(__new_size);
    }
    VectorXd(size_t __new_size, double val) {
        vec.resize(__new_size, val);
    }
    void resize(size_t __new_size) {
        vec.resize(__new_size);
    }
    void resize(size_t __new_size, double val) {
        vec.resize(__new_size, val);
    }
    size_t size() const {
        return vec.size();
    }
    void emplace_back(const double val) {
        vec.emplace_back(val);
    }

    double& operator [] (const int num) {
        assert(num >= 0 && num < (int)vec.size());
        return vec[num];
    }
    double operator [] (const int num) const {
        assert(num >= 0 && num < (int)vec.size());
        return vec[num];        
    }
    VectorXd& operator += (const VectorXd& A) {
        assert(vec.size() == A.vec.size());
        for(int i = 0; i < (int)vec.size(); i++) {
            vec[i] += A.vec[i];
        }
        return *this;
    }
    VectorXd& operator -= (const VectorXd& A) {
        assert(vec.size() == A.vec.size());
        for(int i = 0; i < (int)vec.size(); i++) {
            vec[i] -= A.vec[i];
        }
        return *this;
    }
    VectorXd& operator *= (const double a) {
        for(auto& val : vec) {
            val *= a;
        }
        return *this;
    }
    VectorXd operator + (const VectorXd& A) const {
        assert(vec.size() == A.vec.size());
        VectorXd B = *this;
        return B += A;
    }
    VectorXd operator - (const VectorXd& A) const {
        assert(vec.size() == A.vec.size());
        VectorXd B = *this;
        return B -= A;
    }
    VectorXd operator * (const double a) const {
        VectorXd B = *this;
        return B *= a;
    }

    double dot (const VectorXd& A) const {
        assert(vec.size() == A.vec.size());
        double res = 0.0;
        for(int i = 0; i < (int)vec.size(); i++) {
            res += vec[i] * A.vec[i];
        }
        return res;
    }
    double squaredNorm() {
        return this->dot(*this);      
    }
    double norm() {
        return sqrt(this->squaredNorm());
    }

    void transpose() {
        transposed = !transposed;
    }
};
std::ostream& operator << (std::ostream& os, const VectorXd& A) {
    for(auto& val : A.vec) {
        os << val << " ";
    }
    return os;
}

struct MatrixXd {
    int row_size, column_size;
    std::vector<VectorXd> mat;
    std::vector<VectorXd> L, U;

    MatrixXd() {}
    MatrixXd(int r, int c) : row_size(r), column_size(c) {
        mat = std::vector(r, VectorXd(c));
    }
    MatrixXd(int r, int c, double val) : row_size(r), column_size(c) {
        mat = std::vector(r, VectorXd(c, val));
    }
    void resize(int r, int c) {
        row_size = r; column_size = c;
        mat.resize(r, VectorXd(c));
    }
    void resize(int r, int c, double val) {
        row_size = r; column_size = c;
        mat.resize(r, VectorXd(c, val));
    }

    int get_row_size() const {
        return row_size;
    }
    int get_column_size() const {
        return column_size;
    }

    VectorXd& operator [] (const int num) {
        assert(num >= 0 && num < (int)mat.size());
        return mat[num];
    }
    MatrixXd& operator += (const MatrixXd& A) {
        assert(row_size == A.row_size && column_size == A.column_size);
        for(int i = 0; i < row_size; i++) {
            for(int j = 0; j < column_size; j++) {
                mat[i][j] += A.mat[i][j];
            }
        }
        return *this;
    }
    MatrixXd& operator -= (const MatrixXd& A) {
        assert(row_size == A.row_size && column_size == A.column_size);
        for(int i = 0; i < row_size; i++) {
            for(int j = 0; j < column_size; j++) {
                mat[i][j] -= A.mat[i][j];
            }
        }
        return *this;
    }
    MatrixXd& operator *= (const MatrixXd& A) {
        assert(column_size == A.row_size);
        MatrixXd B(row_size, A.column_size);
        for(int i = 0; i < row_size; i++) {
            for(int k = 0; k < column_size; k++) {
                for(int j = 0; j < A.column_size; j++) {
                    B.mat[i][j] += mat[i][k] * A.mat[k][j];
                }
            }
        }
        return *this = B;
    }
    MatrixXd& operator *= (const double a) {
        for(int i = 0; i < row_size; i++) {
            for(int j = 0; j < column_size; j++) {
                mat[i][j] *= a;
            }
        }
        return *this;
    }
    MatrixXd operator + (const MatrixXd& A) const {
        assert(row_size == A.row_size && column_size == A.column_size);
        MatrixXd B = *this;
        return B += A;
    }
    MatrixXd operator - (const MatrixXd& A) const {
        assert(row_size == A.row_size && column_size == A.column_size);
        MatrixXd B = *this;
        return B -= A;
    }
    MatrixXd operator * (const MatrixXd& A) const {
        assert(column_size == A.row_size);
        MatrixXd B = *this;
        return B *= A;
    }
    MatrixXd operator * (const double k) const {
        MatrixXd B = *this;
        return B *= k;
    }
    VectorXd operator * (const VectorXd& A) const {
        assert(column_size == (int)A.size());
        VectorXd B(row_size);
        for(int i = 0; i < row_size; i++) {
            for(int j = 0; j < column_size; j++) {
                B[i] += mat[i][j] * A[j];
            }
        }
        return B;
    }
    MatrixXd pow(const long long n) const {
        assert(row_size == column_size);
        if(n == 0){
            MatrixXd E(row_size, column_size);
            for(int i = 0; i < row_size; i++) {
                E.mat[i][i] = 1.0;
            }
            return E;
        }
        MatrixXd A = pow(n >> 1);
        A *= A;
        if(n & 1) A *= *this;
        return A;
    }

    VectorXd row(const int num) const {
        assert(num >= 0 && num < row_size);
        return mat[num];
    }
    VectorXd col(const int num) const {
        assert(num >= 0 && num < column_size);
        VectorXd A(row_size);
        for(int i = 0; i < row_size; i++) {
            A[i] = mat[i][num];
        }
        return A;
    }
    void row_apply(const int num, const VectorXd& A) {
        assert(num >= 0 && num < row_size);
        mat[num] = A;
    }
    void col_apply(const int num, const VectorXd& A) {
        assert(num >= 0 && num < column_size);
        for(int i = 0; i < row_size; i++) {
            mat[i][num] = A[i];
        }
    }
    void row_append(const VectorXd& A) {
        assert(column_size == (int)A.size());
        row_size++;
        mat.emplace_back(A);
    }
    void column_append(const VectorXd& A) {
        assert(row_size == (int)A.size());
        column_size++;
        for(int i = 0; i < row_size; i++) {
            mat[i].emplace_back(A[i]);
        }
    }
    VectorXd diag() {
        assert(row_size == column_size);
        VectorXd A(row_size);
        for(int i = 0; i < row_size; i++) {
            A[i] = mat[i][i];
        }
        return A;
    }
    MatrixXd transpose() {
        MatrixXd A(column_size, row_size);
        for(int i = 0; i < column_size; i++) {
            for(int j = 0; j < row_size; j++) {
                A[i][j] = mat[j][i];
            }
        }
        return A;
    }

    void LU_decomposition() {
        assert(row_size == column_size);
        int N = row_size;
        U = mat;
        L.resize(N, VectorXd(N));
        for(int i = 0; i < N; i++) {
            L[i][i] = 1.0;
        }
        for(int i = 0; i < N; i++) {
            for(int j = i + 1; j < N; j++) {
                double val = U[j][i] / U[i][i];
                for(int k = i; k < N; k++) {
                    U[j][k] -= U[i][k] * val;
                }
                L[j][i] += val;
            }
        }
    }
    MatrixXd inverse() {
        assert(row_size == column_size);
        int N = row_size;
        MatrixXd A = *this;
        MatrixXd E(N, N);
        for(int i = 0; i < N; i++) {
            E[i][i] = 1.0;
        }
        for(int i = 0; i < N; i++) {
            double inv = 1.0 / A[i][i];
            for(int j = 0; j < N; j++) {
                A[i][j] *= inv;
                E[i][j] *= inv;
            }
            for(int j = 0; j < N; j++) {
                if(i == j) continue;
                double val = A[j][i];
                for(int k = 0; k < N; k++) {
                    A[j][k] -= A[i][k] * val;
                    E[j][k] -= E[i][k] * val;
                }
            }
        }
        return E;
    }
    VectorXd inverse_product(const VectorXd& A) const {
        assert(row_size == column_size);
        int N = row_size;
        assert((int)A.size() == N);
        VectorXd Y(N);
        for(int i = 0; i < N; i++) {
            double sum = 0.0;
            for(int j = 0; j < i; j++) {
                sum += L[i][j] * Y[j];
            }
            double res = (A[i] - sum) / L[i][i];
            Y[i] = res;
        }
        VectorXd X(N);
        for(int i = N - 1; i >= 0; i--) {
            double sum = 0.0;
            for(int j = N - 1; j > i; j--) {
                sum += U[i][j] * X[j];
            }
            double res = (Y[i] - sum) / U[i][i];
            X[i] = res;
        }
        return X;
    }
};
std::ostream& operator << (std::ostream& os, const MatrixXd& A) {
    for(int i = 0; i < A.row_size; i++) {
        os << A.row(i) << std::endl;
    }
    return os;
}

もっと便利にできるとは思いますが、ひとまずこれで挑みます。(LU分解や逆行列の計算は今回は使いません)
準備ができたので、ベイズ推定のクラスを作ってみます。(これもかなり洗練の余地が残っていると思います)

struct distribution {
    double mu, var;
    distribution(double mu, double var) : mu(mu), var(var) {}
};
struct multi_distribution {
    VectorXd mu;
    MatrixXd cov;
    multi_distribution(VectorXd& mu, MatrixXd& cov) {
        this->mu = mu;
        this->cov = cov;
    }
};
class Bayesian {
private:
    multi_distribution md;
    double log_likelihood;
    double l_sum;

public:
    Bayesian(multi_distribution& md) : md(md) {
        log_likelihood = 0.0;
        l_sum = 0.0;
    }

    void update(const VectorXd& c, double y, const distribution& err) {
        double Syy = c.dot(md.cov * c) + err.var;
        VectorXd Sxy = md.cov * c;
        int siz = c.size();
        MatrixXd res(siz, siz);
        for(int i = 0; i < siz; i++) {
            for(int j = 0; j < siz; j++) {
                res[i][j] = Sxy[i] * Sxy[j];
            }
        }
        res *= (1.0 / Syy);
        double dif = y - c.dot(md.mu);

        md.cov -= res;
        md.mu += Sxy * dif * (1.0 / Syy);
        log_likelihood = -0.5 * (log(2.0 * M_PI * Syy) + dif * dif / Syy);
        l_sum += log_likelihood;
    }

    VectorXd get_mu() const {
        return md.mu;
    }
    MatrixXd get_cov() const {
        return md.cov;
    }
    double get_likelihood() const {
        return log_likelihood;
    }
    double get_sum_likelihood() const {
        return l_sum;
    }
};

平均と分散共分散行列を指定した多変量分布でオブジェクトを作り、update 関数で確率分布や対数尤度を更新します。続いて、以下の calc_cost 関数で各モデルの事後分布から各辺の重みを計算します(INFは適当な巨大な数としてください)。辺の番号は、横の辺を上から順に $0-29$、縦の辺を左から順に $30-59$ としています。

VectorXd calc_cost(const std::vector<Bayesian>& bayes, int turn) {
    VectorXd w(29 * 60);
    int siz = bayes.size();
    double tot_l = 0.0;
    double max_likelihood = -INF;
    for(int m = 0; m < siz; m++) {
        if(bayes[m].get_sum_likelihood() > max_likelihood) {
            max_likelihood = bayes[m].get_sum_likelihood();
        }
    }
    for(int m = 0; m < siz; m++) {
        tot_l += exp(bayes[m].get_sum_likelihood() - max_likelihood);
    }
    for(int m = 0; m < siz; m++) {
        VectorXd cost = bayes[m].get_mu();
        VectorXd dev = bayes[m].get_cov().diag();
        for(int i = 0; i < (int)dev.size(); i++) {
            dev[i] = sqrt(dev[i]);
        }
        double weight = 1.0 - (double)turn / 1000.0;
        cost -= dev * weight;
        cost *= exp(bayes[m].get_sum_likelihood() - max_likelihood) / tot_l;
        for(int i = 0; i < 29 * 60; i++) {
            w[i] += cost[i/29];
        }
    }
    for(int i = 0; i < 29 * 60; i++) {
        w[i] = std::max(w[i], 1000.0);
        w[i] = std::min(w[i], 9000.0);
    }
    return w;
}

辺の長さの計算にあたり、標準偏差にかかる定数を $\gamma = 1.0 - t/1000$ としてみました。ここで $t$ は試行回数です。また、長さを小さく見積もりすぎて負の値になってしまうとダイクストラ法が使えなくなるので、$1000$ 以上 $9000$ 以下という制約を設けることにしました。
尤度の比を計算する際に注意点があります。対数尤度から尤度に直すには指数関数を用いれば良いですが、前述の通り尤度はほとんど $0$ であるため計算できない可能性があります。重要なのは比だけですから、あらかじめそれぞれの対数尤度に適当な値を足してから比を計算するようにしています。これで辺の長さが求まったので、ダイクストラ法で経路を決めることができます。
以上の実装で提出したところ、約94.4%の得点が得られました。コンテストの順位では上位10%より少し下くらいでしょうか。しかし変数を $60$ 個に減らしたり $M=1$ のみを考えたりと制限が多いですから、もっと高い点数を目指せそうです。

どんな場合にベイズ推定が有効なのか？

推定タイプの問題と言っても様々なパターンがあります。ベイズ推定を活用するには、繰り返しにはなりますが

• 推定したい値の確率分布を適切に設定できる
• 観測結果をもとに確率分布を更新できる

ことが重要です。正規分布が便利なことが多いですが、平均と分散共分散行列を適切に決められる必要があります。また、今回は $y = c^T X + \varepsilon$ という線形の関係式を仮定できましたが、関係式がよくわからない場合はどうするのか、という問題もありそうです。

ベイズ推定を応用したガウス過程回帰という手法があります。誤解を恐れずにすごく簡単に説明すると、新しい入力が与えられた時、今までの入力との"距離(どれくらい似ているか)"で重み付けをして出力を予測する手法です。関係式がわからなくても、新しい入力と今までの入力との"距離"を適切に設定できそうなら、ガウス過程回帰を活用できるかもしれません。実際にAHC018ではガウス過程回帰を利用した解法がありました。興味がある方はてりーさんの解説記事やyosssさんの解説スライドをぜひご覧ください。ガウス過程回帰についてはここではこれ以上触れませんが、いつか記事を書くことがあるかもしれません。

長くなってしまいましたが、最後までお読みいただきありがとうございました。数字をいじって結果を眺めて遊んでみたり、色々な問題でベイズ推定を活用できないか考えてみたりして、ベイズ推定と仲良くなりましょう。

*1:これはコインによらず一般にこの結果が得られる確率のことなので、表が出る確率を $p$ と置いて $p$ の全区間で積分して求めることになります。この積分は計算が困難な場合もあるので、計算しなくて済むと大変嬉しいです。

*2:確率分布は全区間で積分して $1$ になる必要があるので、事後分布の積分さえ計算できれば良いです。有名な確率分布の場合は積分の結果が既に知られていることが多いです。その上、そもそも定数を求めなくても確率分布の平均と分散は計算できてしまいます。

*3:問題の入力生成の項目を見ると、辺の長さは本来正規分布に従うわけではなさそうです。とはいえ正規分布を仮定しておけば、試行回数を重ねるごとに、平均は真の辺の長さに近づき分散は小さくなった正規分布が得られると期待できますから、悪くない仮定のように思えます。ただし $M=2$ の場合は簡単ではないです。詳しくはeivourさんの記事をご覧ください。

*4:今回は $60$ 変数としたため同一直線上の辺を全て同一視していますから、同じ辺を複数回通ることがあり得ます。全ての辺を考えた $1740$ 変数のモデルならば、ベクトル $c$ の各要素は $1$ または $0$ をとります。

*5:導出方法によっては $V' = (V^{-1} + cR^{-1} c^T)^{-1}$ となっているかもしれませんが、Woodburyの公式より記事中の式と等しいです。Woodburyの公式についてやどちらの式が望ましいかは記事中のリンクをご覧ください。