Search on the blog

2017年4月30日日曜日

はじパタ第5章 k最近傍法(kNN法)

 平井有三先生の「はじめてのパターン認識」第5章を読んだ。
kNN自体はもちろん知っていたが、理論的な考察についてはほとんど知らなかったので勉強になった。

面白かったところ
  • 式を使ったボロノイ図の定義
  • kNNとベイズ誤り率の関係
  • kNNにおける次元の呪いの話
理解度チェックリスト
  • 最近傍法とは?
  • k最近傍法とは?
  • ボロノイ図とは?
  • kNNを使うときに問題となる次元の呪いについて説明せよ。
    • 予測データと学習データの距離はどうなる?
  • kNNの計算量は?
  • kNNの計算量を減らすための工夫をいくつか説明せよ。

2017年4月27日木曜日

英会話フレーズ: get back to

「あとで〜するね」というときに使える便利なフレーズ。

ミーティングなどで「あとで確認して連絡します。」という場合は、
I'll get back to you later.

「金曜までに具体的な数字を教えてくれますか?」とお願いするときは
Can you get back to me with some figures by Friday?

何かを注文しようとして、その場で即決できず「ちょっと考えてあとで連絡するよ」と言いたいときは
Thank you. Let me think about it and get back to you.

2017年4月22日土曜日

単位超球から一様サンプリング

 誰もが一度はモンテカルロ法を使って円周率を計算してみたことがあると思います。二次元の場合は、点が円の中に入ったり、入らなかったりするのですが、次元が大きくなると円の中に入る頻度が極端に少なくなります。
 
 じゃあ、もし高次元の単位超球内に一様分布する点をサンプリングしたくなったらどうすればよいでしょうか?

 以下自分がやってみたこと。

分散を計算するときの桁落ち対策

分散 = 二乗平均 - 期待値の二乗

という式を使って分散を計算すると、右辺の2つの項の値が非常に近い場合、桁落ちが発生する。

たとえば以下のようなデータを考える。

x1 = 1010
x2 = 1010 + 1
x3 = 1010 + 2
x4 = 1010 + 3
x5 = 1010 + 4

この場合、二乗平均も期待値の二乗も1020程度である。
両者の差は2だが、倍精度浮動小数点の有効桁数は15桁程度なので、この差の情報は失われてしまう。

分散は平行移動しても変わらないので、平均が0になるように中心化を行う。

x1 = -2
x2 = -1
x3 = 0
x4 = 1
x5 = 2

すると、二乗平均は2、期待値の二乗は0となり、差の情報は失われない。

以下のプログラムを実行してみると、桁落ちの様子を確認できる。

#include <iostream>
#include <vector>
#include <iomanip>

using namespace std;

// データを中心化する
vector<double> centerize(const vector<double> &v) {
  int n = v.size();
  double avg = 0.0;
  for (auto &x : v)
    avg += x;
  avg /= n;
  vector<double> w(n);
  for (int i = 0; i < n; i++)
    w[i] = v[i] - avg;
  return w;
}

// 分散を計算する
double calc(const vector<double> &v) {
  int n = v.size();
  double avg_sq = 0.0;
  double avg = 0.0;
  for (auto &x : v) {
    avg += x;
    avg_sq += x * x;
  }
  avg_sq /= n;
  avg /= n;

  return avg_sq - avg * avg;
}

// データ生成
vector<double> gen() {
  vector<double> v;
  for (int i = 0; i < 5; i++)
    v.push_back(1e10 + i);
  return v;
}

int main(int argc, char *argv[])
{
  vector<double> v = gen();
  cout << fixed << setprecision(15);
  cout << calc(v) << endl;
  cout << calc(centerize(v)) << endl;
  return 0;
}

以下実行結果。
$ ./main       
0.000000000000000
2.000000000000000

SRM 712 Div1 600 AverageVarianceSubtree

問題
ノードに重みがつけられた木が与えられる。
木のすべての部分木について重みの分散を計算し、この分散の平均値を求めよ。

解法
DFSしてノード上でDPする頻出テクニックを使う問題だが、分散を計算するために一工夫必要。すべての部分木の分散を効率よく計算するためには(ルート, 部分木のサイズ)ごとに"重みの二乗和"と"重みの和の二乗"を保持すればよい。

例として、ノード(a, b)からなる木とノード(c, d)からなる木をマージする作業を考えてみる。
木(a, b)の二乗和と木(c, d)の二乗和がわかっていたとすると、マージした木の二乗和はそのままこの二つを足せばいいだけなので簡単。

問題は、和の二乗の方で、
(a + b + c + d)^2 = a^2 + b^2 + c^2 + d^2 + 2ab + 2ac + 2ad + 2bc + 2bd + 2cd
(a + b)^2 = a^2 + 2ab + b^2
(c + d)^2 = c^2 + 2cd + d^2
となるので、単純に足すわけにはいかない。

和の二乗 = 二乗の和 + 2 * 異なる項の積和

となっているので、異なる項の積和をうまく部分問題から計算できればよいことが分かる。

木(a, b)の異なる項の積和 = ab
木(c, d)の異なる項の積和 = cd
木(a, b, c, d)の異なる項の積和 = ab + ac + ad + bc + bd + cd = ab + cd + (a + b)(c + d)

となるので、部分問題の異なる項の積和部分問題の和を使えば計算できることが分かる。

以上をふまえて、
  • sum[v][i] = ノードvを根とするサイズiの木の重み和
  • sum_sq[v][i] = ノードvを根とするサイズiの木の重み二乗和
  • sum_prd[v][i] = ノードvを根とするサイズiの木の重みの異なる項の積和
を葉から根方向にDPさせればよい。木をマージするときは、マージする相手方のサイズの分だけ自身が繰り返し現れるので、その分を掛けるのを忘れないようにする。

あとこの問題は数値計算の精度が厳しいらしく、計算誤差を小さくするため以下の工夫が必要。
  • あらかじめ重みを中心化しておく。
  • long doubleを使う。

ソースコード

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

int n;
long long sz[55][55];
double long sum[55][55];
double long sum_sq[55][55];
double long sum_prd[55][55];
double long w[55];
vector<int> child[55];

void dfs(int v) {
  REP (i, n+1) {
    sz[v][i] = 0;
    sum[v][i] = sum_sq[v][i] = sum_prd[v][i] = 0.0;
  }

  sz[v][1] = 1;
  sum[v][1] = w[v];
  sum_sq[v][1] = w[v] * w[v];

  for (auto &c : child[v]) {
    dfs(c);
    for (int i = n; i > 1; i--) {
      for (int j = 1; j < i; j++) {
        int k = i - j;
        sz[v][i] += sz[v][j] * sz[c][k];
        sum[v][i] += sum[v][j] * sz[c][k] + sum[c][k] * sz[v][j];
        sum_sq[v][i] += sum_sq[v][j] * sz[c][k] + sum_sq[c][k] * sz[v][j];
        sum_prd[v][i] += sum_prd[v][j] * sz[c][k] + sum_prd[c][k] * sz[v][j] + sum[v][j] * sum[c][k];
      }
    }
  }
}

class AverageVarianceSubtree {
  public:
  double average(vector<int> p, vector<int> weight) {
    n = weight.size();
    double long avg = 0.0;
    REP (i, n) avg += weight[i];
    avg /= n;
    REP (i, n) w[i] = weight[i] - avg;
    REP (i, n) child[i].clear();
    REP (i, p.size()) child[p[i]].push_back(i+1);

    dfs(0);

    long long tot = 0;
    long double ret = 0.0;

    REP (v, n) {
      for (int i = 1; i <= n; i++) {
        tot += sz[v][i];
        ret += sum_sq[v][i] / i - (sum_sq[v][i] + 2 * sum_prd[v][i]) / i / i;
      }
    }
    
    return ret / tot;
  }
};

2017年4月10日月曜日

はじパタ第4章 確率モデルと識別関数

 平井有三先生の「はじめてのパターン認識」を読んでいる。今日は第4章を読んだ。


悩んだところ
  • ピマ・インディアンデータのROCカーブの説明が間違ってるのでは?と思ったら、病気の人は陰性、健康な人は陽性というデータ定義だったらしい(つまり病気の検査をするときの一般的な陽性・陰性の定義と逆だった)。よく見ると注釈に書いていた。
面白かったところ
  • 共分散行列の固有ベクトル方向にデータを回転して、固有値でスケーリングする技は知っていたが、「白色化」という名前を初めて知った。
  • 独立であることと、相関がないことは同値ではないらしい。独立なら相関はないが、相関がないからといって独立とは限らない。
  • マハラノビス距離は知っていたが多次元正規分布の指数の部分に現れているのは気づかなかった。確かに指数の部分はxが現れる確率密度を表すものなのでマハラノビス距離が遠いものは現れにくいということに対応している。
  • 最尤推定を使った正規表現のパラメータ推定は1回は自分でちゃんと計算しておいた方がよい。
  • 正規分布を使った識別関数が、条件を加えることで、二次局面、線形、最近傍法と変化する話が面白かった。
  • 2つのクラスの共分散行列は同じと仮定しないといけない場合、一般にはそうでないことが多い。そのような場合は、クラスをまとめたものを共通の共分散行列として使うといいらしい。
試してみたこと

Pythonで白色化をしてみた。

import numpy as np
import pandas
from matplotlib import pyplot as plt

mean = [0.0, 0.0]
sigma = [[1, 2], [2, 5]]
data = np.random.multivariate_normal(mean, sigma, 300)
x, y = data.T
plt.plot(x, y, 'o')


df = pandas.DataFrame(data)
u = df.mean()
sigma = df.cov()
lam, S = np.linalg.eig(sigma)

assert np.linalg.norm(np.dot(np.dot(S.T, sigma), S) - np.diag(lam)) < 1e-9

whitening = np.dot(np.diag(lam**(-1/2)), S.T)
data_whitened = np.dot(whitening, (df - u).T).T
x, y = data_whitened.T
plt.plot(x, y, 'o')
Out[88]:


理解度チェックリスト
  • データの前処理についてそれぞれの手法を説明せよ。
    • 標準化
    • 無相関化
    • 白色化
  • パラメトリックモデルとノンパラメトリックモデルの違いは?
  • 多次元正規分布に従うn個のデータを使って、分布のパラメータを最尤推定せよ。

2017年4月9日日曜日

Google Code Jam Qualification Round 2017

 Code Jam 2017の予選があった。ABCの3問を解いて、通過ラインを超えることができた。Dの問題を復習しておく。

問題
Fashion Show

考察
  • 各モデル達はチェスのルーク、ビショップ、クイーンのどれかに対応する
  • 8クイーン問題的に考えることができる
  • ビショップとルークの問題を独立して考えることができる
  • クイーンのポイントは2点なので、ビショップとルークの問題を解いてそのままマージすればOK
  • 各配置問題は(行, 列)または(左斜め軸、右斜め軸)の二部グラフのマッチングをすればOK

ソースコード

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

int n;
int m;
int bd[100][100];
int be[100][100];

class BipartiteMatching {
  int V;
  vector<vector<int> >G;
  vector<int> match;
  vector<bool> used;
  
  bool dfs(int v) {
    used[v] = true;
    for (int i = 0; i < (int)G[v].size(); i++) {
      int u = G[v][i];
      int w = match[u];
      if (w < 0 || (!used[w] && dfs(w))) {
        match[v] = u;
        match[u] = v;
        return true;
      }
    }
    return false;
  }
  
public:
  BipartiteMatching(int v_size) : V(v_size), G(V), match(V), used(V) {}
  
  void add_edge(int u, int v) {
    G[u].push_back(v);
    G[v].push_back(u);
  }
  
  int count() {
    int ret = 0;
    fill(match.begin(), match.end(), -1);
    for (int v = 0; v < V; v++) {
      if (match[v] < 0) {
        fill(used.begin(), used.end(), false);
        if (dfs(v))
          ++ret;
      }
    }
    return ret;
  }

  int getPair(int v) {
    return match[v];
  }
};

void solve() {
  cin >> n >> m;
  memset(bd, 0, sizeof(bd));
  REP (i, m) {
    char t;
    int y, x;
    cin >> t >> y >> x;
    --x, --y;
    if (t == 'x') bd[y][x] = 1;
    else if (t == '+') bd[y][x] = 2;
    else if (t == 'o') bd[y][x] = 3;
  }

  REP (i, n) REP (j, n) be[i][j] = bd[i][j];

  // rook
  BipartiteMatching rook(2 * n);
  set<int> used;
  REP (i, n) REP (j, n) if (bd[i][j] & 1) {
    used.insert(i);
    used.insert(j + n);
  }
  REP (i, n) REP (j, n) if (!used.count(i) && !used.count(j+n)) rook.add_edge(i, j+n);
  rook.count();
  REP (i, n) {
    int j  = rook.getPair(i);
    if (j != -1) {
      be[i][j-n] |= 1;
    }
  }
  
  // bishop
  BipartiteMatching bishop(4*n);
  used.clear();
  REP (i, n) REP (j, n) if (bd[i][j] & 2) {
    used.insert(n-1+i-j);
    used.insert(2*n-1+i+j);
  }
  REP (i, n) REP (j, n) {
    if (!used.count(n-1+i-j) && !used.count(2*n-1+i+j))
      bishop.add_edge(n-1+i-j, 2*n-1+i+j); 
  }
  bishop.count();
  REP (i, n) REP (j, n) {
    int d1 = n-1+i-j;
    int d2 = 2*n-1+i+j;
    if (bishop.getPair(d1) == d2)
      be[i][j] |= 2;
  }

  // output score
  int score = 0;
  REP (i, n) REP (j, n) {
    if (be[i][j] & 1) ++score;
    if (be[i][j] & 2) ++score;
  }
  cout << score << " ";

  // output changed arrangement
  int cnt = 0;
  REP (i, n) REP (j, n) cnt += bd[i][j] != be[i][j];
  cout << cnt << endl;

  REP (i, n) REP (j, n) {
    if (be[i][j] != bd[i][j]) {
      char t = be[i][j] == 3 ? 'o' : (be[i][j] == 1 ? 'x' : '+');
      cout << t << " " << i+1 << " " << j+1 << endl;
    }
  }
}

int main() {
    ios_base::sync_with_stdio(0);
    int T;
    cin >> T;
    REP (i, T) {
        cerr << "Case #" << i+1 << ": " << endl;
        cout << "Case #" << i+1 << ": ";
        solve();
    }

    return 0;
}