mots quotidiens.
Daichi Mochihashi (持橋大地) daichi <at> ism.ac.jp by hns, version 2.10-pl1.

先月 2018年07月 来月
1 2 3 4 5 6 7
8 9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30 31

2018年05月16日(水) [n年日記]

#1 DataGlyph

増井さんの インターフェイスの街角 を読んでずっと以前から気になっていた, QRコードに本来代わるような Xerox PARCで開発されたデータの画像化法 "DataGlyph" を, 夜中にふと思い立って C言語でエンコーダとデコーダを簡単に実装してみました。
% glyph ファイル でファイルをエンコードし, % unglyph ファイル でエンコードされた文字列(改行は任意)をもとの内容に戻します。 EUC専用なので, 端末をEUCモードにしてから実行して下さい。
glyph.c  unglyph.c

やっていることは"ロリロリ変換"と基本的に同じで, データの1/0をそのまま表示 しているだけです。
DataGlyphは, データを何気なく印刷物に潜ませる洗練された方法だなぁ, ということに 感銘を受けていて, Xeroxのページ にも説明のPDFがあります。 より詳しくは こちらのページ にあり, とにかく1/0を機械が判別さえできればよいので, 線の太さを変えたり, 色を変えたりすることで, ごく普通の写真や印刷物にも全くわからないようにデータを潜ませる ことができる技術で, 本当に凄いと思います。

ただこれが広まらず, QRコードのようなあまり美しくない方法が広まっているのは, QRコードの方が目のようなガイドがあって, ロバストだからでしょうか。 ただ, DataGlyphでも任意に読み取ったうち, どこから読めばデータなのかを自動的に 判定させることは技術的にはできる気がしますので, Xeroxの広報力の問題だったの かもしれません。残念です。
とはいえ, いつかこの技術が使われないかなあ, と僕は思っています。 こういうコードを書いていると, それでは生のビットではなくてzlibで圧縮 したり, エラー訂正を含めたくなるので, 少なくとも学習用としては面白い課題ですね。


2017年09月29日(金) [n年日記]

#1 呉越同舟

前々から, マンションにはかなりの人が住んでいるはずなのに, エレベーターで 鉢合わせすることが期待よりかなり少ないのが不思議だった。
今日地下鉄に乗る時に少し考えてみて, これは確率的には確かにそうなのかも しれないと納得できた。 以下の話はネットワークのパケットの衝突の話と同じなので, 恐らく専門的な理論 があると思いますが, 簡単な議論として。

抽象的に考えると, 問題は区間 [0,1] に長さ p (<< 1) の短冊をN個落とした ときに, ある短冊がどれとも重ならない確率と等しい。 1日の時間が区間 [0,1] を, 短冊がそれぞれの住人がエレベーターを使っている時間を表している。

| □    □ □             □      |
これは, N個の短冊がすでにランダムに [0,1] 上にばらまかれている時, [0,1] 上に一様分布でピンを落とした時にどの短冊にもヒットしない確率と等価。

よって, エレベーターを使っている時間を仮に1分とすると, 1日は60分x24時間 = 1440分なので, 上の確率pは1/1440になる。実際には夜中に使う人は少ないので, 若干大き目に見積もるとp=1/1000=0.001くらい。
すると, 「N個の短冊のどれにもヒットしない」という確率は, (1-p) のN乗だから, 1-p=0.999とすれば, 計算すると

となる。今のマンションは住人は100人くらいだと思われるので, かち合う確率が 1/10程度というのは経験とも近い。 なお人口1000人のタワーマンションだと, この確率は0.368になり, 7/10は住人と遭遇するようです。
朝夕の時間に特に出入りが激しくなるので, 上の確率を2倍にして p = 1/500 としてみても, (1-p)^100 = 0.818 になり, 5回に4回は住人と衝突しない。
というわけで, また長年の疑問の一つが解決したのでした。


2017年06月16日(金) [n年日記]

#1 Mathematica autoload

Mathematicaには, MATLABのように特定のディレクトリにユーザー関数を置いて おくと自動的に使える, という機能がないようです(これはRも同じ)。
少し試したところ, 下のようなファイルを書いて init.m でロードすると, 指定したディレクトリにある全てのファイルを読み込むことができるようです。 << はMathematica的にはGetの略記らしい。
(* Library.ma *)
(* $Id: library.ma,v 1.2 2017/06/16 08:56:05 daichi Exp $ *)

dir = Directory[]
SetDirectory["~/share/math"]

Map[ Get,
     Select[ FileNames[],
             Function[x, !DirectoryQ[x] && (x != "library.ma")
                                        && (x != "init.m")]]]

SetDirectory[dir]
init.m は僕の場合(Mac)では ~/Library/Mathematica/Kernel/init.m にありました。 このディレクトリは長いので, 個人的に ~/share/math/ にシンボリックリンクを 張って使っています。
これで, ユーザー定義関数を書いて無事に機能を増やせるようになりました。 かなり基本的な話だと思うのですが, ほとんど書かれていないのが不思議です。

ところで, Python でいう zip の機能, すなわち ((a,b),(x,y)) というリストから ((a,x),(b,y)) というリストを作る機能がMathematicaでも欲しいと思って 書こうかと思っていたところ, 先に調べたら

In[1]:= x = {{"a","b"},{"x","y"}}
Out[1]= {{a, b}, {x, y}}
In[2]:= Transpose[x]
Out[2]= {{a, x}, {b, y}}
でできることを知って, まさに目から鱗でした。


2017年05月01日(月) [n年日記]

#1 Bayesian inference of Logistic regression

下で書いたような離散確率の時系列の他に, 特に社会科学などで, 観測値が ロジスティック正規分布に従っている場合が多くあると思います。 つまり具体的には, 観測値 y ∈ {1,0} (1:生起, 0:非生起)だとして,
y 〜 Bernoulli(σ(x)) = Bernoulli (1 / (1 + exp(-x)))
x 〜 N(0,σ2)
になっているようなモデル。回帰モデルでは x がさらに wTxと回帰になっている場合を考えますが, 議論は基本的に同じです。

これは多項分布の場合はいわゆる対数線形モデルで, 自然言語処理では通常 gradientを計算してL-BFGSやSGDなどの最適化で解くことが多いと思います。 ただし, 最適化の前提となる共変量xが既知ではなく, 学習途中に決まる 潜在変数だったりすると, 最適化してしまうと最初に変な局所解にトラップされて しまい, 学習がうまく動かなくなります。

このため, 例えばBleiの Dynamic Topic Modelではロジスティック正規分布を 下から変分近似して解く方法が示されています。 僕もなぜか長い間, ロジスティック正規分布はこうやって近似するしか ベイズ推定する方法はないのかと思っていたのですが, 最近研究で必要があって 調べたところ, 補助変数を使ってまったく問題なく正しくベイズ推定できることを 知りました。以下, すでに知っているという方は飛ばして下さい。
論文は Groenewald&Mokgatlhe (2005) "Bayesian computation for logistic regression"です。 実は Paul Damien, Jon Wakefield, Stephen Walker(1999), "Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables", JRSS B. の方が著者は有名で, 本質的に同じ内容が 一部書かれていますが, この論文は最初見たときさっぱり分からず, 最初の論文の方が回帰の場合, 多項分布の場合や順序変量の場合など一般的な場合を丁寧に扱っていてわかりやすい ので, そちらの方がお薦めです。

さて, 基本的な戦略はプロビット回帰の場合の有名な Albert&Chib (1993)と同じで, 「観測値を得るために本来あるはずの変数」を潜在変数にして, 補助変数として サンプリングすること。
以下は x を直接推定する場合を考えていますが, xが wTxのように回帰になっている場合も同じで, 論文にはむしろこちらが書かれています。

具体的には, ロジスティック回帰で p=σ(x)=1/(1+exp(-x)) のとき, 観測値 y が1だということは, [0,1] の一様乱数 u があって u < p だったということなので,

p(y|x) = (σ(x))y (1-σ(x))1-y
の補助変数 u との同時確率は,
p(y,u|x)
= I[u≦σ(x)]y I[u>σ(x)]1-y …(1)
= I[x≧log (u/(1-u))]y I[x<log (u/(1-u))]1-y …(2)
になります。
よって, (1)式から u に依存する項を取り出せば, p(u|y,x)∝p(u)p(y,u|x) なので になります。 同様に x の場合は, 上の(2)式から になります。 これから, u の分布は σ(x) を端点の片方とした一様分布, x の事後分布は事前分布p(x)が正規分布N(0,1)だったとしたとき, 正規分布をyに応じて右側あるいは左側を切り落とした truncated Gaussian になる, ということがわかります。
一般には y は1つではなく, Y = [1 1 0 1] のように複数の観測値があるのが普通 なので, この場合は上式の掛け算を考えると, truncated Gaussian は
p(x|Y,U) ∝ p(x)I(a<x<b),
a = maxi log ui/(1-ui)  for yi=1
b = mini log ui/(1-ui)  for yi=0
になります。なお, 論文ではこれが掛け算(Π)ではなく, 突然Σで書かれていて, 恐らく間違っているので注意が必要です。Σにしてしまうと, 上式のmaxやminが 導けないはずです。

というわけで, MATLABでプログラムを書いて確かめてみました。 スクリプトはこちらです。 blogis.m
なお, truncated Gaussianからの生成は, 単純に正規分布から生成してから, 範囲に入るまでrejectすればokです。原理的には, 幅が小さいときは効率が悪くなる ので, Chopin (2012)ではテーブルを使って高速化する方法が示されているようです。

y=[1]y=[0 0 1 0]y=[1 1 1 1 1 1 1 1 1 1 1 1 1]

下のyが観測値です。 これからわかるように, 綺麗にガウス分布状となる事後分布が推定できています。 y=111111111111のような分布は横が歪んでいるので, 単純な正規分布ではないことが わかります。
上で書いたように, 回帰係数 w を推定する場合も方法は同じなので(論文参照), これでロジスティック回帰も問題なくベイズ推定することができるので安心しました。


2017年04月27日(木) [n年日記]

#1 Discrete probability time series

トピックモデルや協調フィルタリングなどで, 「離散分布の時系列」を考えたく なることは時々あるのではないかと思います。 一番典型的な例はトピックモデルで, 時刻tのトピック分布 θt が時刻t-1のトピック分布 θt-1 に依存して決まっている, というマルコフモデル。
離散分布の一番基本的な分布はディリクレ分布なので, 最も素直には, これにはαをスカラーとして
θt 〜 Dir(αθt-1)
とすればよいように思えます。 実際, NTT岩田君の2009年のIJCAIの論文 "Topic Tracking Model for Analyzing Consumer Behavior" でも基本的にこの方法が使われています。(ただしαが時変なので, 下の議論とは 正確には異なっています。後述)

これは一見よさそうに思えるのですが, 少し考えていて, 実はこれは時系列モデル としては適切でないことに気付きました。
いま, 問題を最も簡単にして, 二項分布の確率 ptの時系列を考えてみます。これは上の二次元版なので,

pt 〜 Be(αpt-1,α(1-pt-1))
になります(Be()はベータ分布)。 ベータ分布 Be(a,b) からの乱数は, スケール係数1 *1 のガンマ分布からの乱数 γ1〜Ga(a,1),γ2〜Ga(b,1)を生成して p = γ1/(γ12)とすればよいので, 上記の時系列モデルでは Ga(αpt-1,1)のような乱数を引くことになります。

ここでGa(a,1)はaが小さい値のとき, 0方向にきわめて偏った分布なので, pがすでに0に近い値のとき, 次の時刻の確率の値が今より大きくなる確率は pが小さいほど, どんどん小さくなっていきます。つまり, pが小さくなる方向に バイアスがかかることになります。

ということで, 実際にMATLABでプログラムを書いて試してみました。
初期値がp=0.5のとき, 上記の時系列モデルで, α=1,10,100のときの結果が以下です。

α=1のとき

α=10のとき
α=100のとき

α=1やα=10では, 予想通り確率が最終的に0または1に吸い込まれ, 戻って これなくなることがわかります。 α=100では吸い込まれることは少ないですが, 逆にスケールが非常に大きい(=中間の 確率密度がきわめて高い)確率分布を仮定していることになるので, 確率が0や1に寄る ことができず, ほとんど0.5前後の値だけをふらふらすることになります。

というわけで, 多項分布の時系列には多項分布をそのままディリクレ分布で考える のは不適切で, 他の分布, たとえばガウス分布をロジスティック変換したもの (ロジスティック正規分布)のようなものを考えた方がいいように思えます。 ちなみに上の岩田君の論文では, αは固定ではなく, 時刻t毎に推定するとしている ので, この問題は回避されていると思います。
実際にロジスティック分布, つまりガウス分布の通常のランダムウォークを ロジスティック変換した時系列は以下のようになります。

α=0.5のとき α=1のとき α=10のとき
ここでαはランダムウォークの幅です。(x(t)=x(t-1)+αN(0,1)) α=0.5のときは一度小さい値になった後も移動がゆっくりですが, 末端で確率が回復 しています。またα=1,10のときは確率が1または0に振れても, またそこからどんどん 戻ってきていることがわかります。
他の可能性として, そもそもディリクレ分布(ディリクレ過程)は上の作り方から わかるように, 正規化ガンマ分布(ガンマ過程)なので, ガンマ分布のスケールパラメータを er倍にして動かすという方法があると思います (rは正規乱数)。 単純にガウス過程をロジスティック変換するという方向(対数ガウスCox過程)もあります。 ただいずれにしても結局正規分布やそれに似たものが入ってくるので, 時系列モデルを考える際には連続分布としてのガウス分布は避けられないのでは ないか, という気がします。 *2 そういう意味で, トピック分布の時系列にロジスティック正規分布を考える Dynamic Topic Model は(推定が変分ベイズによる近似で非常に難しいという点を除けば), 妥当なことをしているのではないかと改めて思いました。


*1: 1でなくてもよく, 共通であればよいですが, ここでは簡単のため1とします。
*2: ガンマ分布を使ったランダムウォークもMATLABのスクリプトを書いてみましたが, 容易に0に縮退するか∞に発散するかのどちらかになりました。

2017年04月07日(金) [n年日記]

#1 iTHMM

僕が去年出した無限木構造HMMや教師なし形態素解析が こちらのページ で実装されているのを知りました。 iTHMMは自然言語処理の人でも詳しい内容を理解できている人はほとんどいないと 思われるので, 驚き。

iTHMMについて, こちらで 図入りでものすごく詳しく解説してくれています。嬉しいですね。 実装には前提知識がないと半年以上, 上の記事を書くにも一週間以上かかっている とのことで, そうだろうなと思います。
途中で出されている2つの疑問について,

ちょうどiTHMMをガウス分布に対応させるコードを書いていたところだったので (TACL等でこれまで時間がなかった), 非常に驚きました。
他にも多数の実装が公開されているようで, 素晴らしいですね。


6 days displayed.
タイトル一覧
カテゴリ分類
 なかのひと
Powered by hns-2.10-pl1, HyperNikkiSystem Project