中華料理店過程 (CRP: Chinese Restaurant Process) のJulia 1.0実装
観測データ$x_1,\cdots,x_N \sim \mathcal N(x; \mu, \Lambda^{-1})$について， クラスタ$s_1,\cdots,s_N \sim \mathrm{C}\mathrm{R}\mathrm{P} (\alpha)$を生成するプログラム．
観測データが中華料理店の客，クラスタ番号がテーブル番号に対応する．
Juliaは直感的に書けて中々楽しい．
using Distributions using Plots using LinearAlgebra function sample_CRP(X::Matrix{Float64}, alpha::Float64) N = size(X, 2) p = Vector{Float64}([alpha]) S = zeros(Int32, N) K = 0 for n in 1:N k = rand(Categorical(p / (n-1+alpha))) if k - 1 == K # 新しいクラスタが生成された insert!(p, K + 1, 1) K = K + 1 S[n] = K else # 既存のクラスタが生成された p[k] = p[k] + 1 S[n] = k end end # p = p / (N+alpha) return S end mu = zeros(Float64, 2) covmat = Diagonal{Float64}(I, 2) X = rand(MvNormal(mu, covmat), 1000) S = sample_CRP(X, 2.0) K = size(unique(S), 1) # クラスタ数 plt = histogram(S, bins=K) savefig(plt, "figure.png")
参考文献
 | わかりやすいパターン認識 続／石井健一郎／上田修功【合計3000円以上で送料無料】 価格:3,520円 |
 | ノンパラメトリックベイズ 点過程と統計的機械学習の数理 （機械学習プロフェッショナルシリーズ） [ 佐藤 一誠 ] 価格:3,080円 |
Python2のprint文をVimで置換する正規表現
以下のような変換をかけたいときに使う
他人からもらったコードがPython2ベースだったけど print文まわり置換すれば使えそうかなーーーーーみたいな シチュエーション
Before:
print 'Hello World' print 1 + 2 + 3 print message
After:
print('Hello World') print(1 + 2 + 3) print(message)
使ったコマンド
:%s/print\s\(.*\)$/print\(\1\)/g
ガウス過程からサンプリングを行う最小のPythonプログラム
モチベーション
ガウス過程は勉強したことあったけど，理論と実装がイマイチ頭の中で結びつかなかったので色々調べていた． いきなり回帰とか実装するのは重たいので，手始めにガウス過程からの関数のサンプリングを実装してみた．
実装
ガウス過程は連続関数を生成する確率分布だけどプログラムに落とすときは生成される関数を離散な関数として扱うみたい． データ長も有限になるので，ガウス過程と言いつつ実装上は めっちゃ高次元のガウス関数からめっちゃ高次元のベクトルをサンプリングしている形になる． それにカーネル法を組み合わせている感じ．
フーリエ変換を実装するときに離散フーリエ変換使うけど信号の長さが有限になるので 実質的にフーリエ級数展開になってる感覚に近い？
以下は平均関数 $m(x) = 0$，カーネル関数 $k(x, x') = \exp\left(-|x - x'|^2/σ^2\right)$の ガウス過程$\mathcal{GP}\left(m(x), k(x, x')\right)$からサンプル$y(x)$を10個サンプリングするPythonコード．
import numpy as np from matplotlib import pyplot as plt def mean_function(x): return np.zeros_like(x) def covariance_function(x1, x2, s): return np.exp(-(x1 - x2) ** 2 / s ** 2) x = np.linspace(-10, 10, 100) x1, x2 = np.meshgrid(x, x) sigma = 1.0 m = mean_function(x) gram_matrix = covariance_function(x1, x2, sigma) plt.figure(figsize=(12, 8), facecolor='w') # 関数を10個サンプリングしつつプロット for k in range(10): sample = np.random.multivariate_normal(m, gram_matrix) plt.plot(x, sample, label=f'Sample {k}') plt.legend(loc='upper right')
実行結果 ($σ = 1.0$)
実行結果 ($σ = 0.25$)
JupyterLabでMATLABとRを使えるようにする方法 (2018/9/20 Julia v1.0追加)
背景
自分の主力はPythonだけど，大学院に入ってからMATLABとRを使う機会がかなり増えた． 授業中の演習でJupyterLabを動かせるとメモを取りつつコードも書けて便利だったので環境構築の方法をメモ．
環境
JupyterLabのインストール
最近のAnaconda環境だと標準で入っているっぽい？ 入っていないなら下記コマンドでインストールできる．
$ conda install -c conda-forge jupyterlab
OpenCVのインストール
使用頻度が高いのでついでに．下記コマンドでインストールできる．
$ conda install -c conda-forge opencv
JupyterLabからMATLABカーネルを使えるようにする
参考：https://jp.mathworks.com/help/matlab/matlab_external/install-the-matlab-engine-for-python.html
以下のパスは例なので，自分の環境に合わせて変えること．
- Windowsの場合
$ pip install matlab_kernel $ cd C:\Program Files\MATLAB\R2018a\extern\engines\python $ python setup.py install
- Macの場合
$ pip install matlab_kernel $ cd /Applications/MATLAB_R2018a.app/extern/engines/python $ python setup.py install
MATLABをPythonから呼び出す例
Pythonスクリプトと同じディレクトリに以下の内容のファイルmatlab_function.mを用意しておく．
function out=matlab_function(x) out = x .^ 2; end
MATLABエンジンの起動は以下の通り．
import matlab import matlab.engine eng = matlab.engine.start_matlab()
ファイルを置いておくだけでMATLABの関数を呼び出せるようになる．
for k in range(0, 5): print(eng.matlab_function(k))
出力：
0 1 4 9 16
Numpyのndarrayを渡すときは変換が必要なので注意．
import numpy as np a = np.arange(10, dtype=np.float64) a = matlab.double(a.tolist()) a = eng.matlab_function(a) print(a)
出力：
[[0.0,1.0,4.0,9.0,16.0,25.0,36.0,49.0,64.0,81.0]]
JupyterLabからRを使えるようにする
Rをインストールした後，RのコンソールからIRkernelをインストールすれば良い．
参考：
- Rのインストール R: The R Project for Statistical Computing
- RStudioのインストール Home - RStudio
- Rのコンソール上で以下のコマンドを実行し，Rのカーネル (IRkernel) を入れる
> install.packages('devtools') > devtools::install_github('IRkernel/IRkernel') > IRkernel::installspec()
おわり
追記：JupyterLabからJulia v1.0を使えるようにする
Julia v1.0は以前のバージョンと違って，Pkg.add()ではなくpkgモードを使う．
ターミナルでjuliaを実行してJuliaを起動した後，]を打つとpkgモードとなり，
プロンプトが(v1.0) pkg>のようになるので，そこでadd IJuliaを実行すればOK．
