記事・メモ一覧

[mathjax] 最尤推定とMAP推定とベイズの定理は繋がっていたので、 記憶が定かなうちに思いの丈を書き出してみるテスト。俯瞰してみると面白い。 あるデータ達(x)が観測されていて、それらは未知のパラメータを持つ確率分布から発生している。 観測されたデータ達(x)を使って、それらのデータを発生させたモデルのパラメータを推定したい。 確率密度関数の中に2つの変数があって。 片方を定数、片方を確率変数として扱うことで2通りの見方ができる。 例えば(n)回のコイントスで(k)回表が出る確率が(theta)だとしてベルヌイ分布の確率密度関数は (k)と(theta)のどちらが確率変数だとしても意味がある。 begin{eqnarray} f(k;theta)=theta^k (1-theta)^{n-k} end{eqnarray} 表が出る確率(theta)が定数だと思って、確率変数(x)の確率密度関数と思う。 単なる確率変数(x)の確率密度関数の中に(theta)という定数がある。尤度。 尤度は確率変数(x)の確率密度関数!。 begin{eqnarray} p(X=x|theta) end{eqnarray} 尤度(p(x|theta))を最大にする(theta)を推定するのが最尤推定。 begin{eqnarray} newcommand{argmax}{mathop{rm arg~max}limits} hat{theta} = argmax_{theta} p(X=x|theta) end{eqnarray} 事後確率と事前確率には関係があって以下のようになる。ベイズの定理。 begin{eqnarray} p(theta|x)=frac{p(x|theta) p(theta)}{p(x)} end{eqnarray} ちなみに、(p(x))は以下のようにしておくとわかりやすい。 同時確率と周辺確率の関係。表を書いて縦、横がクロスするところが同時確率だとして、 縦、横いずれかの方向に同時確率を足し合わせる操作にあたるらしい。 なにか確率変数が独立でなければならない、というのは気にしない。 begin{eqnarray} p(x) = int p(x,theta) dtheta end{eqnarray} なので以下みたいに書き直せる。最後の比例のところは...。 左辺は事後確率分布、右辺は尤度と事前確率分布の積!!。 begin{eqnarray} p(theta|x) &=& frac{p(x|theta) p(theta)}{p(x)} \\ &=& frac{p(x|theta) p(theta)}{int p(x,theta) dtheta} \\ &propto& p(x|theta) p(theta) end{eqnarray} (p(theta))は確率変数(theta)の確率分布。尤度(p(x|theta))は(theta)に対して定数。(x)に対して変数。 ということで、右辺は確率分布(p(theta))を尤度(p(x|theta))を使って変形した確率分布。 で、左辺の(p(theta|x))は右辺の(p(x|theta) p(theta))を定数倍した確率分布。 データを観測していない状態で立てた(p(theta))があって、 観測したデータを使って求めた尤度(p(x|theta))が得られたことで、 左辺の(p(theta|x))が得られた、という状況。 (p(theta|x))は確率変数(theta)の確率分布なので、 最尤推定とベイズの定理を俯瞰してみると、最尤推定が点推定である一方で、 ベイズの定理では確率分布が得られるという具合で異なる。 (観測値が極端なデータだったとき、最尤推定は極端な推定結果が得られるだけだけれども、 ベイズの定理で得られる事後確率分布は確率分布なので様子がわかる..??) 事後確率分布を最大化する(theta_{MAP})を求めるのがMAP推定。(点推定) begin{eqnarray} hat{theta_{MAP}} = argmax_{theta} p(theta|x) end{eqnarray} 尤度をかけて得られた事後分布と同じ形になる便利な分布があって、 観測データ達の分布と対応して決まっている(共役事前分布)。 ベルヌイ分布の共役事前分布はベータ分布。

以前、\"正規分布に従う確率変数の二乗和はカイ二乗分布に従うことの証明\"という記事を書いた。 記事タイトルの通り、正規分布に従う確率変数の二乗和はカイ二乗分布に従う。 [clink url=\"https://ikuty.com/2018/08/01/chi-square-distribution\"] 実際にデータを生成して確かめてみる。 まずは、scipi.stats.chi2.pdfを使って各自由度と対応する確率密度関数を書いてみる。 \"pdf\" ってPortableDocumentFormatではなく、ProbabilityDensityFunction(確率密度関数)。 import numpy as np import matplotlib.pyplot as plt from scipy import stats # 0から8まで1000個のデータを等間隔で生成する x = np.linspace(0, 8, 1000) fig, ax = plt.subplots(1,1) linestyles = [\':\', \'--\', \'-.\', \'-\'] deg_of_freedom = [1, 2, 3, 4] for k in deg_of_freedom: linestyle = linestyles[k-1] ax.plot(x, stats.chi2.pdf(x, k), linestyle=linestyle, label=r\'$k=%i$\' % k) plt.xlim(0, 8) plt.ylim(0, 1.0) plt.legend() plt.show() 定義に従って各自由度ごとに2乗和を足し合わせてヒストグラムを作ってみる。 (ループのリストが自由度) cum = 0 for i in [1,2,3,4]: # 標準正規分布に従う乱数を1000個生成 x = np.random.normal(0, 1, 1000) # 2乗 x2 = x**2 cum += x2 plt.figure(figsize=(7,5)) plt.title(\"chi2 distribution.[k=1]\") plt.hist(cum, 80) 全部k=1ってなってしまった...。左上、右上、左下、右下の順に1,2,3,4。 頻度がこうなので、各階級の頻度の相対値を考えると上記の確率密度関数の形になりそう。 k=1,2が大きくことなるのがわかるし、3,4と進むにつれて変化が少なくなる。

PHP7.4の新仕様 FFI,Foreign Function Interface について。Python的な...。 PHP7の時点で既にバイトコードになっていてPHP8からJITによるネイティブ化だったような。 PHP5.xから7.0で劇的に速くなった記憶があり、速度を追求するような機運があるのは確か。 投票システムの不備で\"賛成24、反対15,過半数越えだから採択\"みたいになり、 その後過半数から2/3に更新された、という経緯から、望まれない仕様なのは確かそう。 WordPressが速くなるのは嬉しいけども不毛すぎ..。 脱PHPが必要。 FFI is one of the features that made Python and LuaJIT very useful for fast prototyping. It allows calling C functions and using C data types from pure scripting language and therefore develop “system code” more productively. For PHP, FFI opens a way to write PHP extensions and bindings to C libraries in pure PHP. <?php / create FFI object, loading libc and exporting function printf() $ffi = FFI::cdef( \"int printf(const char *format, ...);\", // this is regular C declaration \"libc.so.6\"); // call C printf() $ffi->printf(\"Hello %s!n\", \"world\"); // create gettimeofday() binding $ffi = FFI::cdef(\" typedef unsigned int time_t; typedef unsigned int suseconds_t; struct timeval { time_t tv_sec; suseconds_t tv_usec; }; struct timezone { int tz_minuteswest; int tz_dsttime; }; int gettimeofday(struct timeval *tv, struct timezone *tz); \", \"libc.so.6\"); // create C data structures $tv = $ffi->new(\"struct timeval\"); $tz = $ffi->new(\"struct timezone\"); // calls C gettimeofday() var_dump($ffi->gettimeofday(FFI::addr($tv), FFI::addr($tz))); // access field of C data structure var_dump($tv->tv_sec); // print the whole C data structure var_dump($tz);

vagrantを知っていてこれからdockerに入門する方向のチートシートを書く。 dockerコンテナの方が小さいので基本的に粒度が合わないが、 対称性を重視して似た機能を並べてみた。 コマンド一覧(vagrant list-commands)からメジャーなのをピックアップした。 Dockerコンテナのネットワーク設定等、概念違いで比較しづらいものは除いた。 また、各コマンドにはパラメタが必要だが煩雑になるため省略した。 vagrant box / docker image操作 vagrant docker boxの一覧vagrant up imageの一覧docker images boxの追加vagrant box add */*imageの追加docker pull imageの詳細docker inspect * boxの削除vagrant box remove *imageの削除docker rmi * boxの更新vagrant box update */*imageの更新docker pull * 全boxの更新vagrant box update imageからbox初期化vagrant init * imageのタグ設定docker tag * vagrant 仮想マシン操作/コンテナ生成/起動/停止 vagrant docker 仮想マシン起動vagrant upコンテナ生成/起動docker run コンテナバックグラウンド実行docker run -d コンテナ起動docker start 仮想マシン終了vagrant haltコンテナ停止docker stop 仮想マシン再起動vagrant reloadコンテナ再起動docker restart 仮想マシン一時停止vagrant pauseコンテナ中断docker pause 仮想マシン再開vagrant resumeコンテナ再開docker unpause 仮想マシン削除vagrant destroyコンテナ削除docker rm 仮想マシンステータス表示vagrant status 全マシン一覧vagrant global-status稼働コンテナ一覧docker ps -a 仮想マシンログインvagrant ssh スナップショット(vagrantのみ) vagrant スナップショット作成vagrant snapshot save * スナップショット復元vagrant snapshot restore * スナップショット削除vagrant snapshot delete * スナップショット一覧vagrant snapshot list dockerコンテナはホストと共有する範囲が大きく、全環境で完全に同じように動かすことは 難しいかもしれない。 vagrantとdockerは排他的な存在ではなく、vagrant上にdockerコンテナっていう構成もありえる。 Vagrantfile, Dockerfile関連は次回...

昔使っていた.tmux.confが失われてしまったため、新しいものを作っていく。 不要な機能が満載になっていたので、これを機に必要なものだけ設定する。 複数の鯖にsshで繋ぐと、ssh接続のペインだらけになり、間違えて操作してしまうことが多い。 その対策として、ssh接続ごとにペインの背景の色を変えていた。 結局、間違えるときは色を変えたところで間違えるため、デザインがケバくなるだけで無意味。 ペインに接続先ホスト名を表示するぐらいがちょうど良い感じだと思う。 こちらの完成度が超絶高いので使うことにしてみた。 ホスト名が強調されたりするともっと良いのだけれども。 SSH接続先のホスト名をtmuxのペインに表示する(大事件) #!/bin/zsh #/usr/local/bin/tmux-pane-border if [[ $1 = \"ssh\" ]]; then pane_pid=$2 info=$({ pgrep -flaP $pane_pid ; ps -o command -p $pane_pid; } | xargs -I{} echo {} | awk \'/ssh/\' | sed -E \'s/^[0-9]*[[:blank:]]*ssh //\') port=$(echo $info | grep -Eo \'-p ([0-9]+)\'|sed \'s/-p //\') if [ -z $port ]; then local port=22 fi info=$(echo $info | sed \'s/-p \'\"$port\"\'//g\') user=$(echo $info | awk \'{print $NF}\' | cut -f1 -d@) host=$(echo $info | awk \'{print $NF}\' | cut -f2 -d@) if [ $user = $host ]; then user=$(whoami) list=$(awk \' $1 == \"Host\" { gsub(\"\\\\.\", \"\\\\.\", $2); gsub(\"\\\\*\", \".*\", $2); host = $2; next; } $1 == \"User\" { $1 = \"\"; sub( /^[[:space:]]*/, \"\" ); printf \"%s|%sn\", host, $0; }\' ~/.ssh/config ) echo $list | while read line; do host_user=${line#*|} if [[ \"$host\" =~ $line ]]; then user=$host_user break fi done fi ssh_hostname=\" ssh:$user@$host \" fi echo $ssh_hostname # .tmux.conf set-option -g pane-border-status bottom set-option -g pane-border-format \"#P #(tmux-pane-border #{pane_current_command} #{pane_pid})\" 上記のサイトの通り、各ペインのセッションがsshであればペインのボトムにhostnameを表示している。 sshコマンド自体に手を入れたくないので、process listをポーリングするやり方に賛同します。 次回はシステムのクリップボードコピー。

[mathjax] 球面集中現象理解のための数学シリーズ第2弾。 前の記事でデカルト座標->極座標の変換から体積要素の積分により3次元球体の体積を導出してみました。 極座標の3変数(r,phi_1,phi_2)について定積分を計算していくと(frac{4pi r^3}{3})が出てきます! [clink url=\"https://ikuty.com/2019/08/02/concentration_on_the_sphere_1/\"] より高次元の超球体の体積がわかると、超高次元のときに体積が超球体の表面に集中する様が 論理の飛躍なく理解できるらしいので、今回は高次元の超球体の体積を導出するやつを流してみます。 (n)次元空間において以下を満たす点の集合が半径(r)の(n)次元球体です。 (x,y,z...)とやっていくとアルファベットが尽きるので(x_1,x_2,cdots,x_n)とします。 (3次元だとすると、半径が(sqrt{x_1^2+x_2^2+x_3^2}le r)の全ての点。) begin{eqnarray} x_1^2+x_2^2+cdots+x_n^2 le r^2 end{eqnarray} この集合を全部足したものが体積になります。 例えば3次元のとき、前回の記事の通り以下のようになります。ここではデカルト座標(x,y,z)。 極座標変換することで計算ができます。 begin{eqnarray} int_{x^2+y^2+z^2 lt r^2} U(x,y,z) dV &=& int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{r} U(r,phi_1,phi_2)r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 \\ &=& frac{4}{3} pi r^3 end{eqnarray} ちなみに2次元球体の体積(つまり円の面積)は当たり前ですが(pi r^2)です。 1次元球体の体積(つまり直線の長さ)は(r)です。 なんとなく(n)次元だと(r)の次数が(n)になりそうですが、 実際そうで(時間がないので公式として使ってしまいます)、 (n)次元の超球体の体積は(r^n)に比例します。 (Vr^n)のようにしておきます。 (n)次元の超球体を輪切りにすることを考えます。 球をどうやって輪切りにしても断面図は真円ですので...。 3次元の球体を輪切りにすると2次元の平面が現れる、というイメージです。 (n)次元の超球面を輪切りにし(n-1)次元の平面を作成し、 (n-1)次の平面と微小距離(Delta r)をかけて(n)次元の直方体を作ります。 微小距離(Delta r)を限りなくゼロに近づけることで体積要素とします。 体積要素は詰まるところ輪切りですが、輪切りなのに(n)次元なのです。 体積要素を全範囲で積分することで体積を求めます。 さて、(x_1^2 + x_2^2 + cdots + x_n^2 = 1)が半径1の単位球です。 これを輪切りにします。(n)次の変数を定数にします。(x_n=t)。 つまり、輪切りは(x_1^2 + x_2^2 + cdot + x_{n-1}^2 = 1-t)。 輪切りの半径は(sqrt{1-t})。(t=0)であれば丁度球体を真っ二つにする感じです。 (t=1)であれば球体のキワの限界で面積がゼロ。 (t)で輪切りにしたあと、(t)を微小距離(t+Delta t)だけ増やします。 その体積は(V_{n-1}(sqrt{(1-t)^{n-1}}) Delta t)。 (Delta t)を限りなくゼロに近づけると真の輪切りに近づき、 それを(r)の全範囲で定積分します。 begin{eqnarray} V_n=int_{-1}^{1} V_{n-1}sqrt{(1-t)^{n-1}} dt end{eqnarray} 球体は真ん中で対照なので、 begin{eqnarray} V_n =2 int_{0}^{1} V_{n-1}sqrt{(1-t)^{n-1}} dt end{eqnarray} 計算してから積分するのと、積分してから計算するのが同じなので、 begin{eqnarray} V_n=2 V_{n-1} int_{0}^{1} sqrt{(1-t)^{n-1}} dt end{eqnarray} 超絶懐かしい置換積分を使うと(int f(x) dx = int f(g(t)) frac{dx}{dt} dt)が出来る。 (t=sintheta)とすると、 begin{eqnarray} V_n= 2 V_{n-1}(int_{0}^{frac{pi}{2}} cos^{n-1}theta cos theta d theta ) end{eqnarray} これどうやって積分するんだよ...、と呆然としてしまうけれども、 公式があるようです。時間がないので公式を使います! 階乗が2つ並んでいるのは1つ飛ばしで階乗をする2重階乗というらしい。 begin{eqnarray} int_{0}^{frac{pi}{2}} sin^n x dx &=& int_{0}^{frac{pi}{2}} sin^{n-1}x sin x dx \\ &=& frac{pi}{2} frac{(n-1)!!}{n!!} end{eqnarray} ちなみに(cos^n x)の積分は(sin^n x)の積分から導かれて以下のようになる。 begin{eqnarray} int_{0}^{frac{pi}{2}} sin^n x dx = int_{0}^{frac{pi}{2}} cos^n x dx = begin{cases} frac{(n-1)!!}{n!!} & nが奇数 \\ frac{pi}{2} frac{(n-1)!!}{n!!} & nが偶数 end{cases} end{eqnarray} 体積の話に戻ると、 begin{eqnarray} V_n = begin{cases} 2V_{n-1}frac{(n-1) (n-3) (n-5) cdots 2}{n (n-2) cdots 3} & nが奇数 \\ pi V_{n-1} frac{(n-1)(n-3)(n-5)cdots 3}{n (n-2) (n-4) cdots 2} & nが偶数 end{cases} end{eqnarray} 偶数のときと奇数のときで別になってしまっているのを1つにしたい。 偶数のとき(n=2k)、奇数のとき(n=2k-1)とおいて、それが等しいとする。 begin{eqnarray} V_{2k-1} &=& 2V_{2k-2} frac{(2k-2)(2k-4)(2k-6)cdots 2}{(2k-1) (2k-3) cdots 3 } \\ V_{2k} &=& pi V_{2k-1} frac{(2k-1)(2k-3)(2k-5)cdots 2}{2k (2k-2)(2k-4) cdots 2} end{eqnarray} 2つの式で(V_{2k-1})が現れるので、それで等式を立てる。 begin{eqnarray} 2V_{2k-2} frac{(2k-2)(2k-4)(2k-6)cdots 2}{(2k-1) (2k-3) cdots 3 } = frac{V_{2k}}{pi} frac{2k (2k-2)(2k-4) cdots 2}{(2k-1)(2k-3)(2k-5)cdots 2} end{eqnarray} ガンガン約分する。 begin{eqnarray} V_{2k} = 2pi V_{2k-2} frac{1}{2k} = frac{pi V_{2k-2}}{k} end{eqnarray} (V_{2k})と(V_{2k-2})の漸化式になっていて、(V_{2k-2})が(V_2)になるまで再帰的に計算する。 begin{eqnarray} V_{2k} = frac{pi^{k-1}}{k!} V_2 = frac{pi^k }{k!} end{eqnarray} 階乗の一般化(ガンマ関数,(n! = Gamma{(n+1)}))を使って書くと、 begin{eqnarray} V_n = frac{pi^{frac{n}{2}}}{Gamma{(frac{pi}{2}+1)}} end{eqnarray} (n)次の超球面の体積が出ました...。 次回、これを使って超高次元でメロンパンの皮が厚いのを示します..。

LaravelをAPIサーバにして同期的にsklearnのPCAを実行するアプリを作ってみました。 jQyery/bootstrap/chart.jsがフロント、APIサーバはLaravel+MySQL。 Laravel製APIがGET/POSTに対してPythonコードを実行します(Shellで...)。 exec()でPythonを起動するため無茶苦茶重いし、 ろくにエラーハンドリングできません。 結論から書けば同期的なアプリをこの構造で作るのは無理があります。 バックエンドが無茶苦茶重くてどうせバッチ実行になるのであれば、 上記の問題は結構問題なくなって、これでも良いかなと思い始めます。 MS系のInteroperabilityで、多言語が動的に結合するやつがありますが、 あんな感じでLL言語をglueできれば楽なのになと思います。 PSRの多言語拡張みたいなやつで、PHPからPythonのクラスを使うとか...

[mathjax] 球面集中現象を理解するために記憶にないことが多すぎるので 理解に必要そうなことを少しずつ復習していきます。 まず極座標の直行座標変換。極座標を使って球の体積を求めてみます。 次の記事で球の体積を求める時に使う微小体積がヤコビアンになることを確認します。 極座標の直行座標変換 2次元のとき、極座標((r),(phi)) を直行座標で表すと、 begin{eqnarray} x &=& r cos(phi) \\ y &=& r sin(phi) end{eqnarray} 3次元のとき、極座標((r),(phi_1),(phi_2))を直行座標で表すと、 begin{eqnarray} x &=& r cos(phi_1) \\ y &=& r sin(phi_1) cos(phi_2) \\ z &=& r sin(phi_1) sin(phi_2) end{eqnarray} ここまではイメージできるのだけれども、(N)に飛ばしたときにどうなるのか。 N次元のとき、極座標((r),(phi_1),(phi_2),(cdots),(phi_{N-1}))を直行座標で表すと 以下のようになるらしい。 begin{eqnarray} x_1 &=& r cos(phi_1) \\ x_2 &=& r sin(phi_1) cos(phi_2) \\ x_3 &=& r sin(phi_1) sin(phi_2) cos(phi_3) \\ vdots \\ x_{N-1} &=& r sin (phi_1) cdots sin(phi_{N-2}) cos(phi_{N-1}) \\ x_{N} &=& r sin (phi_1) cdots sin(phi_{N-2}) sin(phi_{N-1}) end{eqnarray} この後使わないけれど、逆変換は以下の通りになるらしい。 begin{eqnarray} r &=& sqrt{x_1^2 + cdots + x_{N-1}^2 + x_N^2} \\ phi_1 &=& arccosfrac{x_1}{sqrt{x_1^2 + cdots + x_{N-1}^2 + x_N^2}} \\ phi_2 &=& arccosfrac{x_2}{sqrt{x_2^2 + cdots + x_{N-1}^2 + x_N^2}} \\ vdots \\ phi_{N-2} &=& arccos frac{x_{N-2}}{sqrt{x_{N-2}^2 + x_{N-1}^2 + x_N^2}} \\ phi_{N-1} &=& begin{cases} arccos frac{x_{N-1}}{sqrt{x_{N-1}^2+X_N^2}} & X_N ge 0 \\ -arccos frac{x_{N-1}}{sqrt{x_{N-1}^2+X_N^2}} & X_N lt 0 end{cases} end{eqnarray} 3次元極座標系における体積計算の例 極座標((r),(phi_1),(phi_2))において(r)を微小量(Delta r)だけ増やす。 (2pi r)は半径を(r)とする円の直径だということを考慮すると、 (phi_1)から微小角(Delta phi_1)増やしたときの(r Delta phi_1)は円弧の一部。 (Delta phi_1)を限りなく小さくすると(r Delta phi_1)は直線と考えられる。 (2pi rsin(phi_1) )が半径を(rsin(phi_1))とする円の直径と考えると、 (phi_2)から微小角(Delta phi_2)増やしたときの(r sin(phi_1) Delta phi_2)は円弧の一部。 (Delta phi_2)を限りなく小さくすると(r sin(phi_1) Delta phi_2)は直線と考えられる。 (Delta r, Delta phi_1, Delta phi_2)が限りなくゼロに近くときに 上記のそれぞれの直線を辺とする直方体ができる。その体積を(dV)とすると、 begin{eqnarray} dV &=& Delta r cdot r Delta phi_1 cdot r sin(phi_1) Delta phi_2 \\ &=& r^2 sin(phi_1) Delta r Delta phi_1 Delta phi_2 end{eqnarray} 直交座標(x,y,z)の位置を(U(x,y,z))とする。 直交座標系における球の体積はざっくり以下で表せる。 begin{eqnarray} int U(x,y,z) dV end{eqnarray} 直行座標の極座標変換は前述の通り以下。この座標を(U(r,phi_1,phi_2))とする。 begin{eqnarray} x &=& r cos(phi_1) \\ y &=& r sin(phi_1) cos(phi_2) \\ z &=& r sin(phi_1) sin(phi_2) end{eqnarray} 極座標系において以下だから、 begin{eqnarray} 0 le r le 1 \\ 0 le phi_1 le pi\\ 0 le phi_2 2pi end{eqnarray} それぞれについて積分すると体積になる。 begin{eqnarray} int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{1} U(r,phi_1,phi_2)r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 end{eqnarray} つまり以下。 begin{eqnarray} int U(x,y,z) dV = int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{1} U(r,phi_1,phi_2)r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 end{eqnarray} ここで( U(r,phi_1,phi_2) =1 )として積分範囲を(r)とすると球の体積の公式になる。 begin{eqnarray} int U(x,y,z) dV &=& int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{r} r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 \\ &=& int_{0}^{2pi} int_{0}^{pi} Bigl( int_{0}^{r} frac{sin(phi_1)}{3} bigl[ r^3 bigr]_{0}^{r} Bigr) dphi_1 dphi_2 \\ &=& int_{0}^{2pi} int_{0}^{pi} frac{r^3 sin(phi_1)}{3} dphi_1 dphi_2 \\ &=& int_{0}^{2pi} -frac{r^3}{3} bigl[ cos(phi_1) bigr]_{0}^{pi} dphi_2 \\ &=& int_{0}^{2pi} frac{2r^3}{3} dphi_2 \\ &=& frac{2r^3}{3} bigl[ phi_2bigr]_{0}^{2pi} \\ &=& frac{4pi r^3}{3} end{eqnarray} 次回、微小体積とヤコビアンが等しい理由を書いてみます。

[mathjax] sihoutteって何て読むのか...と思うけども\"シルエット\"だそう。フランス語源。 ポートレート写真を背景白、顔を黒に減色した例の\"シルエット\"。\"輪郭\"みたいな。 データ達を複数のクラスタに分割したとして、各々のデータがそのクラスタに存在することの 収まりの良さを表すことができる。 本来(相対的に)他のクラスタに属しているべきデータ達の割合が分かったり、 クラスタ境界で(相対的に)どちらのクラスタに属しても良さそうなデータ達の割合が分かったりする。 sklearnは教師なしクラスタリングまでで、それをグラフ化しないといけないのだけど、 matplotlibに対応するのがなく、\"線を書く\"みたいなコマンドを並べて作っていかないといけない様子。 yellowbrickというパッケージを使うとそれもやってくれる。 Yellowbrick is a suite of visual diagnostic tools called “Visualizers” that extend the Scikit-Learn API to allow human steering of the model selection process. In a nutshell, Yellowbrick combines scikit-learn with matplotlib in the best tradition of the scikit-learn documentation, but to produce visualizations for your models! sklearnとmatplotloibだけで出力するバージョンと、yellowbrickで出力するバージョンの 両方を試してみた(前者はクラスタのラベルを出力できずsihoutteグラフとの対応関係が理解できずに 中途半端で終了)。 sihoutte係数 全データ達は(x_i)。データ(x_i)がクラスタ(A)に属し、最近傍にクラスタ(B)があるという状況。 (A)のクラスタ中心は(x_A)、(B)のクラスタ中心は(x_B)。 \"最近傍のクラスタ中心との距離の平均\"から\"自分のクラスタ中心との距離の平均\"を引いた値。 わざわざ1行に全ての変数が出てくるように書いてみる。 begin{eqnarray} s_i = frac{sum_{i=1}^{n}|x_i-x_B|/n - sum_{i=1}^{n}|x_i-x_A|/n }{max bigl{sum_{i=1}^{n}|x_i-x_A|/n,sum_{i=1}^{n}|x_i-x_B|/n bigr}} end{eqnarray} データが(n)個あるので、sihoutte係数も(n)個できる。 自分が属しているクラスタ(A)よりも、隣のクラスタ(B)の方が居心地が良いと 分子はマイナスになる。 今属しているクラスタでも隣のクラスタでも、どちらもさほど居心地が変わらないとゼロ近辺になる。 大きければ迷いなく今のクラスタで良いことを表せる。 sklearnとmatplotlibだけでsihoutte係数 乱数で作った偽データに6-meansをかけた図。 どうやってもクラスタ中心のラベルが取り出せない!!..(スルー..) import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [[0,0], [20,20], [50,50], [40,30], [40,10], [20,40]] sigma = [ [[30,20],[20,50]], [[20,30],[10,20]], [[60,40],[20,20]], [[60,20],[20,60]] ,[[30,10],[10,30]],[[50,20],[20,50]] ] points = 100 clusters = [] for index in range(len(mu)): cluster = np.random.multivariate_normal(mu[index], sigma[index], points) dig = np.full((points,1),index+1, dtype=int) cluster = np.hstack((cluster,dig)) clusters = np.r_[clusters,cluster] if len(clusters) > 0 else cluster plt.scatter(x=clusters[:,0], y=clusters[:,1],c=clusters[:,2]) from sklearn.cluster import KMeans model = KMeans(n_clusters=6, init=\'random\',max_iter=10) y_km = model.fit_predict(clusters[:,:2]) labels = model.labels_ centers = model.cluster_centers_ centers fig = plt.figure() ax1 = fig.add_subplot(1,1,1) ax1.scatter(x=clusters[:,0], y=clusters[:,1],c=labels) ax2 = fig.add_subplot(1,1,1) ax2.scatter(x=centers[:,0], y=centers[:,1], alpha=0.5,s=600,c=\"pink\",linewidth=2,edgecolors=\"red\") sklearnとmatplotlibだけでsihoutte係数のグラフを出してみる。 \"線を書く\"みたいなコマンドを並べて規格通りの図を作るのか...。 from sklearn.metrics import silhouette_samples # cluster数 num_clusters=6 #全データのsilhouette係数を取得 silhouette_vals = silhouette_samples(clusters[:,:2],y_km,metric=\'euclidean\') cluster_labels = np.unique(y_km) min,max = 0,0 yticks = [] for i,c in enumerate(cluster_labels): c_silhouette_vals = silhouette_vals[y_km == c] c_silhouette_vals.sort() max += len(c_silhouette_vals) plt.barh( range(min,max), c_silhouette_vals, height=1.0, edgecolor=\'none\' ) yticks.append((min+max)/2.) min += len(c_silhouette_vals) avg = np.mean(silhouette_vals) plt.axvline(avg, color=\'red\', linestyle=\"--\") plt.yticks(yticks, cluster_labels + 1) plt.ylabel(\'Cluster\') plt.xlabel(\'Silhouette coefficient\') plt.show() 全体的に茶色のクラスタのsilhouette係数が低め。 残念ながら、どのクラスタと対応するのか出力できず何の考察も出来ず...。 yellowbrickでsilhouette係数を出力 無茶苦茶簡単に出せる。 from yellowbrick.cluster import SilhouetteVisualizer sv = SilhouetteVisualizer(model) sv.fit(clusters[:,:2]) sv.poof() 出てきた図。自作したものと全然合っていないように見える.. 自作したものとラベルが合っていないだけだと信じたい。 似た形の塊があるので.. 以上、完全な失敗だけれども一旦終了...

k-means法を実行する際に妥当なkを決めたいという欲求があります。 クラスタ集合の凝集度を定量化することでkと凝集度の関係を得られます。 複数のkについてkと凝集度の関係を取得し、 そこから妥当なkを発見する方法がエルボー法です。 この記事においては以下を書いてみます。 クラスタ集合の凝集度を定量化する方法 kと凝集度の関係を表すグラフの書き方 クラスタ集合の凝集度を定量化する方法 これ以上無いくらいシンプルなやり方です。 データポイントとcentroidの距離の2乗和(SSE)を凝集度として使います。 それを全クラスタに関して計算し足し合わせます。 足し合わせた値を歪み(distortion)と言っているところが多いです。 クラスタ数kが増えるにつれてこの値は小さくなるだろうと予測できます。 横軸にk、縦軸にSSEの和をとるグラフを書いてみると、 kの増加に伴ってSSEの和が小さくなること、減少幅は次第に小さくなりそうです。 (すべてのケースでそうなるのかは不明) 減少幅が緩やかになる地点のkが費用対効果が高いkですし、 ほぼ減少が止まる(サチる)地点のkを採用することは、 それ以上kを増加させても意味がないという事実を根拠として使えそうです。 擬似データへの適用 前回から使っている多次元正規分布に従う擬似クラスタデータに対してk-meansを適用します。 k=1から9まで適用し、k-distortionグラフを書いてみます。(だいぶPythonに慣れてきた..。) import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [[0,0], [20,20], [50,50], [40,30], [40,10], [20,40]] sigma = [ [[30,20],[20,50]], [[20,30],[10,20]], [[60,40],[20,20]], [[60,20],[20,60]] ,[[30,10],[10,30]],[[50,20],[20,50]] ] points = 100 clusters = [] for index in range(len(mu)): cluster = np.random.multivariate_normal(mu[index], sigma[index], points) dig = np.full((points,1),index+1, dtype=int) cluster = np.hstack((cluster,dig)) clusters = np.r_[clusters,cluster] if len(clusters) > 0 else cluster plt.scatter(x=clusters[:,0], y=clusters[:,1],c=clusters[:,2]) from sklearn.cluster import KMeans kmeans_model = KMeans(n_clusters=6, init=\'random\',max_iter=10).fit(clusters[:,:2]) labels = kmeans_model.labels_ distortions = [] numClusters = 10 for i in range(1,numClusters): kmeans_model = KMeans(n_clusters=i, init=\'random\',max_iter=10) kmeans_model.fit(clusters[:,:2]) distortions.append(kmeans_model.inertia_) plt.plot(range(1,numClusters),distortions,marker=\'o\') plt.xlabel(\'Number of clusters\') plt.ylabel(\'Distortion\') plt.show() 複数の擬似点を中心に2次元正規分布に従う散らばりをもったデータ達です。 k-distortionグラフです。もの凄い期待通りなグラフができました。 k=5か6あたりでdistortionがサチっています。 6-meansを適用し、centroidを重ねてみました。 散布図を見ての想像になりますが、 確かに4-meansと5-meansでは、5-meansの方が凝集していそうです。 5-meansと6-meansだと、4と5程凝集度合いに変化がなさそうです。