Search on the blog

2017年5月7日日曜日

Kさんの教え

 半年間一緒に働いたエンジニアのKさんの教え。

おれたちはソフトウェアを作っている
ハードウェアを作る場合は、綿密な計画、設計無しにものづくりを始めてはいけない。作り直しのコストが非常に大きいからね。でも、おれたちはソフトウェアを作っている。綿密な計画、設計をしているうちに作っちゃった方が早いよね。

必要以上に複雑にするな
ソフトウェアは使ってもらいながら進化していくものなので、最初から複雑にしすぎてはいけない。絶対必要なものに集中してとにかくシンプルなものからスタートせよ。

変更は小さく頻繁に
ソフトウェアの変更は細かく、高頻度に行うべき。一度に多くのものを変えようとするのはよくない。簡単な機能なのにすぐに変更、リリースができないなんて論外。

学びこそ至福の喜び
常に学び続けろ。新しいことを学ぶのがエンジニアの大きな喜びの一つだ。学ぶことをやめてしまったエンジニアに価値はない。

おまえはどうしたいのか
あの人がこう言った、あのチームの人たちはこう言っている、じゃない。
お前はどうしたいんだ?それが一番大事なことだ。

一心不乱に夢中になること
何かに夢中になっているヤツのまわりにいると、こっちもなんだか楽しい気分になってくるぜ。おれの会社の創業者もそうさせてくれるやつだった。

セキュリティと生産性のトレードオフ
セキュリティを厳しくしすぎると生産性が著しく低下する。両者はトレードオフの関係にあるけど、日本では生産性の問題が軽視される傾向にあるみたいだ。

(番外編)マラソン
マラソンはいろんなことを教えてくれる。精神状態をいかにコントロールするか、ペース配分をどうするかなど。

(番外編)やばたん
YABATAN知らないのか?日本語だぜ。若い女の子のことのこと分かってないんじゃないか?

2017年5月4日木曜日

便利な英語表現: To our best knowledge

 論文を読んでいるとよく出てくる表現。「我々が知るかぎり〜だ」と言うときに使える。

To our best knowledge, this is the first time that ....
我々が知るかぎり、....したのはこれが初めてだ。

To our best knowledge, our model is better than all previously published methods.
我々が知るかぎり、我々の手法は過去に発表された他の手法より優れている。

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

はじパタ第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カーブから最適動作点はどうやって選べばよいか?

2017年4月8日土曜日

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月2日日曜日

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パーセンタイルの値以内の領域。

2017年4月1日土曜日

はじパタ第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軸はブートストラップサンプルしたときのユニークなデータ数。

PyStan入門(1)とりあえず動かす

 緑本を読んで、MCMCで統計モデルのパラメータ推定をすることに興味がわいてきたので、PyStanをさわっていこうと思う。
 まずは、とりあえずインストールして動かしてみました。

PyStanのインストール
pipで入ります。

$ pip install pystan

モデルの記述
stanのモデルは.stanという拡張子のファイルに書くのが一般的なようです。他にもPythonのコードにベタ書きするという方法もあるようです。
ここではstanファイルにモデルを記述することにし、以下の内容を8schools.stanというファイルに記述します。

data {
  int<lower=0> J; // number of schools 
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates 
}
parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}

上のモデルを以下のPythonコードから参照します。以下のコードをsample.pyという名前のファイルに記述します。
import pystan
import matplotlib.pyplot as plt

schools_dat = {
 'J': 8,
 'y': [28,  8, -3,  7, -1,  1, 18, 12],
 'sigma': [15, 10, 16, 11,  9, 11, 10, 18]
}

fit = pystan.stan(file='8schools.stan', data=schools_dat, iter=1000, chains=4)

print(fit)
fit.plot()
plt.show()

実行結果
実行してみます。

$ python sample.py
DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal): assignment operator <- deprecated in the Stan language; use = instead.

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_286b3180dfa752c4cfedaf0241add0e4 NOW.
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 1)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 3)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 2)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 0)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 1)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 2)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 3)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 0)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 1)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 0)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 0)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 1)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 2)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 1)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 0)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 3)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 3)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 2)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 0)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 1)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 0)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 1)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 3)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 0)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 2)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 1)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 3)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 1)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 2)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 2)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 0)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 3)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 1)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 3)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 2)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 0)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 3)
# 
#  Elapsed Time: 0.044725 seconds (Warm-up)
#                0.037462 seconds (Sampling)
#                0.082187 seconds (Total)
# 
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 1)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 0)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 2)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 1)
# 
#  Elapsed Time: 0.05144 seconds (Warm-up)
#                0.040057 seconds (Sampling)
#                0.091497 seconds (Total)
# 
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 0)
# 
#  Elapsed Time: 0.04517 seconds (Warm-up)
#                0.046124 seconds (Sampling)
#                0.091294 seconds (Total)
# 
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 2)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 2)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 2)
# 
#  Elapsed Time: 0.053818 seconds (Warm-up)
#                0.040396 seconds (Sampling)
#                0.094214 seconds (Total)
# 
Inference for Stan model: anon_model_286b3180dfa752c4cfedaf0241add0e4.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu         7.81    0.17   5.09  -2.41   4.61   7.77  11.13  17.63    907    1.0
tau        6.63    0.24    5.6   0.29   2.45   5.33   9.17  21.19    531    1.0
eta[0]     0.42    0.02   0.94  -1.53  -0.17   0.44   1.05   2.17   1701    1.0
eta[1]     0.01    0.02   0.87  -1.67  -0.55 3.2e-3    0.6   1.74   1751    1.0
eta[2]    -0.23    0.02   0.89  -1.95  -0.81  -0.22   0.38   1.51   1866    1.0
eta[3]    -0.01    0.02    0.9  -1.82  -0.61-8.1e-3   0.58   1.79   1669    1.0
eta[4]    -0.33    0.02   0.87   -2.1  -0.87  -0.34   0.21   1.44   1603    1.0
eta[5]    -0.21    0.02    0.9   -1.9  -0.82  -0.23   0.39   1.55   1599    1.0
eta[6]     0.34    0.02    0.9  -1.42  -0.24   0.36   0.95   2.07   1537    1.0
eta[7]     0.06    0.02   0.94  -1.84  -0.56   0.07   0.66   1.87   1604    1.0
theta[0]  11.61    0.24   8.19  -1.43   6.15  10.55  15.65  31.15   1211    1.0
theta[1]   7.84    0.14   6.35  -5.56   4.05   7.76   11.7  20.83   1988    1.0
theta[2]   5.97     0.2    7.7  -11.7   1.85   6.38  10.91  20.08   1505    1.0
theta[3]   7.61    0.15   6.45  -5.51   3.46   7.67  11.86  20.35   1897    1.0
theta[4]   5.13    0.14   6.45  -8.65   0.99   5.58   9.56   17.0   2000    1.0
theta[5]   6.21    0.17   6.51  -7.34   2.29   6.46  10.55  18.69   1522    1.0
theta[6]  10.64    0.17   6.88  -1.89   6.09  10.06  14.68  26.06   1623    1.0
theta[7]   8.33     0.2   7.94  -7.88   3.77   8.07  12.54  25.84   1522    1.0
lp__      -4.82    0.11   2.66 -10.43  -6.54  -4.59  -2.84  -0.42    598    1.0

Samples were drawn using NUTS at Sat Apr  1 00:55:24 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

C++のコンパイラが走って、MCMCのサンプリングが走って、パラメータの分布が出たようです。続いて各パラメータの分布図とiterationごとの変化がプロットされます。
細かいことは置いといて、何やら楽しそうなことが出来るということは分かりました。
次回は自分で簡単な乱数生成モデルを作って、それをモデリングしてパラメータ推定ができるかどうかを試してみたいと思います。

参考
Getting started — PyStan 2.14.0.0 documentation

はじパタ第1章 はじめに

 平井有三先生の「はじめてのパターン認識」を読み始めた。パッと目次を見た感じだとだいたい知ってることが多そうだったが、基礎は大事なので丁寧に読み進めたい。
 今日は第1章「はじめに」を読んだ。100次元の超立方体は毬栗のような形をしているという話がおもしろかった。
 以下に理解度チェックリストを作っておく。
  • よいパターン認識とは何かいかのキーワードを使って述べよ。
    • 特徴抽出
    • 識別規則の学習
    • 般化能力
  • 以下の尺度の定義、具体例を述べよ。
    • 名義尺度
    • 順序尺度
    • 間隔尺度
    • 比例尺度
  • ダミー変数とはなにか?
  • 次元の呪いとは?
  • d次元超立方体について以下を計算せよ。
    • 頂点の数はいくつ?
    • 辺の数は?
    • 超立方体を構成するm次元超平面(1 < m < d)の個数はいくつ?

2017年3月30日木曜日

緑本 第11章 空間構造のある階層ベイズモデル

 久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。ようやく最終章までたどり着いた!第11章「空間構造のある階層ベイズモデル」を読んだので、理解度チェックリストを作っておく。
  • 空間相関とは?
    • どのようなときに考慮しないといけない?
    • どのように表現できる?
  • 条件つき自己回帰モデル(conditional auto regressive, CAR)とは?
  • 確率場とは?
    • マルコフ確率場とは?
    • ガウス確率場とは?
  • 観測データに欠測がある場合、空間相関を考慮する・しないで予測結果はどのように変わる?

2017年3月27日月曜日

緑本 第10章 階層ベイズモデル -GLMMのベイズモデル化-

 久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第10章「階層ベイズモデル -GLMMのベイズモデル化-」を読んだ。

 そもそもこの本を読み始めたのは、MCMCをどのようにベイズに適用するのかを理解するのが目的だったのだが、この章に来てようやくイメージができた気がする(9章からぼんやりわかり始めていたものが10章でしっかりと分かった)。この本を読む前までは、MCMCというのを最適化の手法のように勘違い(おそらくメトロポリス法が最適化手法の一つであるsimulated annealingに似ていたのでそう思ってしまったのかも)していた気がする。つまり最尤推定でやっているようにMCMCで対数尤度関数を最適化するパラメータを探すのだと勘違いしていて、それってどうやるのか全然わからない、っていうかわざわざそんなことする意味あるの?と悩んでいたのである。
 しかし、MCMCを「ベイズ」に適用するので、やろうとしていることはパラメータの点推定ではなく、パラメータの分布推定である。パラメータの分布が分かれば、特定パラメータの周辺分布や、応答変数の信頼区間や、ある説明変数が応答変数に影響するかどうかなどの統計検定ができる。こういうことをやるためにMCMCを使うというのが分かると、あぁなるほどな〜と使い方(ベイズのモデルをMCMCで計算するという意味)をイメージできるようになった。

 まえおきが長くなってしまったが、いつもの理解度をチェックするための確認リストを作っておく。
  • 階層ベイズモデルとは?
    • GLMMだと何が難しい?
  • 階層事前分布とは?
  • 超事前分布とは?
  • ベイズ統計モデルで主観的な事前分布の使用を避けたい場合は、どのような分布を使えばよいか?
  • 無情報事前分布と階層事前分布のどちらを使うべきかはどのように判断すればよいか?
    • 大域的なパラメータと局所的なパラメータ

2017年3月26日日曜日

緑本 第9章 GLMのベイズモデル化と事後分布の推定

 久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。第9章「GLMのベイズモデル化と事後分布の推定」を読んだので理解度チェックリストを作っておく。
  • 無情報事前分布とは?
    • 一様分布以外にはどのようなものが使える?
  • データの中央化とは?
  • データの標準化とは?
  • MCMCでサンプリングする際に設定すべき値は?
    • どのようにチューニングするのがよいか?
  • 周辺事後分布とは?
  • ギブスサンプリングのアルゴリズムについて説明してください
    • 全条件つき分布とは?
    • メトロポリス法と比較したときの利点は?
  • 共役事前分布とは?

2017年3月25日土曜日

緑本 第8章 マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデル

 久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。8章「マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデル」を読んだので、自分の理解度確認のためのチェックリストを作っておく。
  • 統計モデルのパラメータqの最尤推定値を求めたいとします。qを解析的に求められないときにはどのような方法を使えばよいか?原始的な方法を説明してください。
  • メトロポリス法のアルゴリズムを説明してください。
  • MCMCとは?
    • どの部分がマルコフ連鎖なのか?
    • どの部分がモンテカルロ法なのか?
    • MCMCの目的は?
      • 何を出力する?
      • どういうときに重宝する?
  • MCMCの定常分布は何を表す?
  • 頻度主義とベイズ統計学の違いについて説明せよ。
    • モデリングにデータを当てはめた結果えられるものはそれぞれ何か?