Page List

Search on the blog

2017年4月21日金曜日

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;
}

2017年4月8日土曜日

はじパタ第3章 ベイズの識別規則

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


 
面白かったところ
  • ベイズのMAP推定のところで、事後確率 = 修正項 × 事前確率というのは当然知っていたけど、修正項のところの概念的な意味づけが面白かった。
  • 図3.1 識別境界とベイズ誤り率、図3.6 ROC曲線と等損失直線群の図が考えさせられた。特に図3.6では最適動作点を探す方法を新しく知れたのでよかった。
  • 尤度比のところの話は単純に式変形としてみるのではなく、何を意味しているのかを考えてみると面白かった。

チェックリスト
  • ベイズの最大事後確率に基づく識別規則を式を書いて説明せよ。
    • 尤度比を使って識別規則を表現せよ。
  • ベイズの識別規則の誤り率を計算せよ。(2クラス識別で)
  • 誤りを犯すことによって発生する危険性がクラスによって異なる事例を述べよ。
    • そういう場合はどうすればよいか?
  • リジェクトとはなにか?
  • ROCカーブはどのようにプロットするか説明せよ。
  • 以下の語句の意味を説明せよ。
    • false positive rate
    • true positive rate
    • precision
    • recall
    • accuracy
    • F-measure
  • AUCとは何か?
    • とりうる値は?どういう値だといい?
  • ROCカーブから最適動作点はどうやって選べばよいか?

Drew Houston Commencement Address MIT 2013

DropboxのFounderであるDrew Houston氏のスピーチを聞いた。

"Drew Houston Commencement Address MIT 2013"
Standard YouTube License

tennis ball
すごいやつらは何かにものすごく熱中している。
テニスボールを投げると一生懸命それを追いかける犬のように。
ものすごく働く人は自制心があるからそうするのではなく、問題を解決することが楽しいからそうする。
大切なのは無理をして頑張ることではなくテニスボールを見つけることだ。

having a circle
あなたは、よく時間を過ごす5人の人の平均。
尊敬できる人に囲まれているというのはとても大切。
MIT、シリコンバレー、ハリウッド...。何をやるにしても、ベストな人が集まるのは世界に1箇所しかない。君たちはそこに行くべきだ。他のどこにも留まってはいけない。
ヒーローに出会い、その人から学べるというのは大きな利点がある。ヒーローも君の5人のサークルのうちの一人だ。

30,000
何かを学ぶ最速の方法は実際にやってみること。
ビルゲイツの最初の会社もスティーブジョブズの最初の会社も成功しなかった。でもこの2人が失敗したことにくよくよ思い悩んだなんて想像するのは難しい。
失敗することは大したことじゃない。1回でも成功すればいい。
人生は30,000日しかない。完璧を目指して、勉強したり、準備をしたりしているうちに時間はすぐに過ぎていく。
人生を完璧なものにするのではなく、面白いものにするという発想を持ってみよう。
完璧を目指すのではなく、自分に冒険をする自由を与えてみよう。そして常に上に進み続けよう。

2017年4月1日土曜日

PyStan入門(2)線形モデル

 PyStanで線形モデルのパラメータ推定をしてみた。
データyはNormal(w*x + b, σ2)から生成されるとしてモデルを組んで、パラメータw, b, σの分布を求めてみた。

分かったこと
  • pystan.stan()メソッドの引数で以下を指定できる
    • iter 各チェーンが行うサンプリング回数(warmup含む)
    • warmup 初回の捨てるサンプリング数
    • chain チェーンの数
  • fit.extract()メソッドでサンプリング値を取れる
    • permuted引数=Trueとすると、チェーンをマージしたサンプル列を取得できる
    • permuted引数=Falseとすると、(iter, chain, parameters)の次元のnumpy.ndarrayを取得できる
コード

data {
    int<lower=0> N;
    vector[N] x;
    vector[N] y; 
} parameters {
    real w;
    real b;
    real<lower=0> sigma;
} model {
    y ~ normal(w * x + b, sigma);
}

import numpy as np
import pystan
import matplotlib.pyplot as plt
import pandas as pd


def gen():
  N = 30
  w = 10
  b = 3
  sigma = 0.5

  x = np.random.rand(N)
  err = sigma * np.random.randn(N)
  y = w * x + b + err

  return (x,y)


def train(x, y):
  data = {
    'N': len(x),
    'x': x,
    'y': y
  }
  fit = pystan.stan(file='linear_model.stan', data=data, iter=2000, warmup=1000, chains=4)
  return fit


def predict(param, x):
  return param['w'] * x + param['b']


if __name__ == '__main__':
  x, y = gen()
  fit = train(x, y)
  fit.plot()
  plt.show()

  xt = np.linspace(0, 1, 100)
  params = pd.DataFrame({
    'w': fit.extract(permuted=True).get('w'),
    'b': fit.extract(permuted=True).get('b')
  })

  median_params = params.median()
  yt = predict(median_params, xt)
  yt_lower = [np.percentile(predict(params, x), 2.5) for x in xt]
  yt_upper = [np.percentile(predict(params, x), 97.5) for x in xt]

  plt.plot(x, y, 'bo')
  plt.plot(xt, yt, 'k')
  plt.fill_between(xt, yt_lower, yt_upper, facecolor='lightgrey', edgecolor='none')
  plt.show()

結果
まずトレースプロットと、各パラメータの分布。それっぽい分布になっている。
次に学習結果を使って予測値のプロット。実線はパラメータの中央値を使ってyを予測した値。グレーの部分は得られたパラメータのサンプルを使って予測値を計算したときの2.5%パーセンタイルと97.5パーセンタイルの値以内の領域。

はじパタ第2章 識別規則と学習法の概要

 平井有三先生の「はじめてのパターン認識」を読んでいる。今日は第2章「識別規則と学習法の概要」を読んだので学習メモを残しておく。

理解度チェック
  • アフィン関数とは?
    • アフィン関数を線形関数にするためによく使われる方法は?
  • バイアスとは?
  • バリアンストは?
  • バイアスとバリアンスの関係は?
  • 次の手法について説明せよ。
    • holdout法
    • cross validation法
    • leave one out法
    • bootstrap法
悩んだところ
バイアス項とバリアンス項(分散項)の導出のところで悩んだ。
一般にE[X2] = E[X]E[X]が成り立つと思い込んでしまっていて、バリアンス項って消えるんじゃない?何これ?となっていた..。
XとYが独立な場合はE[XY] = E[X]E[Y]ですが、XとXは独立じゃないので、バリアンスの項は消えませんね。

基本中の基本を忘れてました。
例えばサイコロに{-2, -1, 0, 0, 1, 2}という数字を書いて、2回投げたときにそれぞれ出る値をX、Yとすると、E(XY) = E(X)E(Y) = 0となりますが、E(XX)は0にはなりませんね。

試してみたこと
N個のデータからN回復元抽出(sampling with replacement)すると、あるデータが少なくとも1回は抽出される確率は1-e-1となるってのが面白かったので試してみた。1-e-1は 劣モジュラ関数の最適化や回路の時定数などでも出現する値で興味深い。

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt


def bootstrap_unique_freq(n):
 x = np.random.randint(0, n, n)
 return len(set(x)) / n


if __name__ == '__main__':
 t = 100  # iteration
 ns = []  # sample size
 xs = []  # unique sample freq

 for n in [10, 100, 1000, 10000]:
  ns.extend(n for _ in range(t))
  xs.extend(bootstrap_unique_freq(n) for _ in range(t))

 sns.boxplot(x=ns, y=xs)
 plt.show()

箱ひげ図を書いてみた。
x軸は元のデータ数、y軸はブートストラップサンプルしたときのユニークなデータ数。