AHC040推定パートで、長方形0を毎回入れると推定結果が悪化する件¶
概要¶
AHC040の各長方形の幅を推定するときに、長方形0を毎回入れた観測をすると推定結果が悪化する件について、調べたことのメモ。
背景¶
- 推定アプローチとして、長方形0を常に左上に配置して、そこから横方向・縦方向に伸ばすという「L字配置」をして、長方形の幅の和の情報を得る方法が考えられた
- しかし、このアプローチで推定した場合、あまり推定結果がよくならない現象が起きていた
簡単な実験¶
簡単のため、N個の幅に対して、N/2個の和(誤差あり)の情報をT回得られる場合で実験。
#include <bits/stdc++.h>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
class Game {
public:
Game(int N, int sigma) : N(N), sigma(sigma), engine(seed_gen()), normal_rand(0, sigma) {
const int U = 100000;
const int L = U / 10 + engine() % (U / 2 - U / 10 + 1);
for (int i = 0; i < N; i++) {
rects.push_back(L + engine() % (U - L + 1));
}
}
int get(const vector<int>& idx) {
int ret = 0;
for (int i : idx) ret += rects[i];
return ret + normal_rand(engine);
}
int debug(const vector<int>& idx) {
int ret = 0;
for (int i : idx) ret += rects[i];
return ret;
}
private:
const int N, sigma;
random_device seed_gen;
mt19937_64 engine;
normal_distribution<> normal_rand;
vector<int> rects;
};
int main() {
const int N = 10;
const int sigma = 5000;
const int T = 100;
Game game(N, sigma);
random_device seed_gen;
mt19937_64 engine(seed_gen());
VectorXd mu(N);
MatrixXd var = MatrixXd::Zero(N, N);
for (int i = 0; i < N; i++) {
mu(i) = 65000;
var(i, i) = 20000.0 * 20000.0;
}
for (int t = 0; t < T; t++) {
set<int> s;
s.insert(0); // !!! 常に0を入れる !!!
while (s.size() < N / 2) s.insert(engine() % N);
vector<int> idx(s.begin(), s.end());
int res = game.get(idx);
// カルマンフィルタで更新
MatrixXd H = MatrixXd::Zero(1, N);
for (int i : idx) H(0, i) = 1;
MatrixXd z(1, 1);
z(0, 0) = res;
MatrixXd R(1, 1);
R(0, 0) = sigma * (double)sigma;
MatrixXd y = z - H * mu;
MatrixXd S = H * var * H.transpose() + R;
MatrixXd K = var * H.transpose() * S.inverse();
mu += K * y;
var -= K * H * var;
}
for (int i = 0; i < N; i++) {
cout << i << "\t";
cout << game.debug({i}) << "\t";
cout << mu(i) << "\t";
cout << sqrt(var(i, i)) << endl;
}
return 0;
}
常にindex0を入れた場合
# index 実際の幅 推定幅 推定標準偏差
0 52999 77140.8 16001.1
1 47006 39634.8 4111.54
2 84034 78201.9 4100.78
3 89342 82800.4 4107.2
4 95932 90371.2 4109.58
5 54637 49830.1 4102.79
6 61122 56394.2 4106.21
7 98155 93883.8 4100.35
8 84121 77498 4103.31
9 72824 64949 4114.36
すべてランダムに選んだ場合
0 50084 50967.6 952.936
1 71539 71651.8 918.389
2 95361 96866.1 925.838
3 75503 76197.4 959.127
4 88144 88641.6 938.597
5 46058 45540.7 920.785
6 86431 85634.8 947.211
7 76880 74930.8 925.045
8 33301 33224.4 923.575
9 66842 67633.7 927.25
観測行列Hのランク¶
- 観測行列Hを、「各観測を行に並べたもの」(T x N行列)で考える
- 行列Hのrankを計算してNになっていない場合、数値計算の精度低下や、観測データが状態変数の組み合わせでしか識別できず各変数を個別に識別できない問題が発生してしまう
- 特異値分解してみると、最小の特異値が0にかなり近く、ほぼランク落ちしてしまっていることがわかる
- wikipediaにもあるように、LU分解などでランク計算するのは精度落ちなどがありえ、用いられないみたい
- AHC040の実際ので見てみると、特異値の小さい方の2つが、ゼロに近いほどではないものの、他と比べて小さめになってるっぽい
MatrixXd H = MatrixXd::Zero(T, N);
MatrixXd z(T, 1);
MatrixXd R = MatrixXd::Zero(T, T);
for (int t = 0; t < T; t++) {
// ランダムに選ぶ
set<int> s;
s.insert(0); // !!! 常に0を入れる !!!
while (s.size() < N / 2) s.insert(engine() % N);
vector<int> idx(s.begin(), s.end());
int res = game.get(idx);
for (int i : idx) H(t, i) = 1;
z(t, 0) = res;
R(t, t) = sigma * (double)sigma;
}
JacobiSVD<Eigen::MatrixXd> svd(H, Eigen::ComputeThinU | Eigen::ComputeThinV);
cout << svd.singularValues() << endl;
常にindex0を入れた場合
16.716
6.34501
5.79118
5.59349
5.31227
5.10989
4.86206
4.53603
4.11643
9.06548e-16 // !!! 0ではないがかなり小さい !!!
すべてランダムに選んだ場合
15.9283
6.08912
5.92107
5.73654
5.51335
5.44418
4.95822
4.75776
4.33734
3.8955
結論¶
- 長方形0を常に入れてしまうと、観測行列Hがランク落ち(かそれに近い状況になって)しまい、推定精度が悪化してしまうため、うまくバラけるようにHを設計する必要があった
- ランク落ちしているかどうかの確認は、LU分解ではなく、特異値分解などで確認したほうがよい