記事・メモ一覧

個人的には、言語に入門する際にはオンラインより書籍が適切と考えている。 そして最初の一冊はなるべく薄いものを選ぶべきと考えている。 読んだ後、あえて文章として書き出すと良さそう。行間なりツボを自分の言葉でおまとめ。 あそこにこう書いてあった..というように思い出せると良いかなと。 まずは制御構造。 if ... elif ... else ブロックの開始はコロン。ブロックはインデントで表現。 else ifは繋げて書く。elif,elseはオプション。 x = 20 if x < 0: print("負の値") elif x == 0: print("ゼロ") elif x == 1: print("いち") else: print("その他") for 初期値、終了値、増分を書くCスタイルと違う。シーケンスに対してforeachをかけるスタイルのみ。 ループにかけたからといってシーケンスがコピーされることはなく変更することができる。 当然不安定なのでシーケンスをコピーしてループに使うべき、。 words = [1,2,3,4,5] for x in words: print(x) # for x in words[:] words.insert(0,x) シーケンスを最後まで読み終わった後に実行するブロックを定義できる。 forと同じインデント位置にelseを書く。 words = [1,2,3,4,5] for x in words: print(x) else: print(\"finish\") 1 2 3 4 5 finish range ループはシーケンスを処理するのみなので、シーケンスが無い場合は作る必要がある。 nからn-1ならrange(n)。初期値、終了値、ステップはrange(初期値,終了値,ステップ)。 range(n)はシーケンスを作るのではなく反復可能体を返す。つまり、range()の応答時には シーケンスはメモリ確保されておらず評価時に初めてメモリ確保される。ref.C++ iterator。 for x in range(10): print(x**2) 0 1 4 9 16 25 36 49 64 81 for x in (range(0,10,3)) print(x) 0 3 6 9 シーケンスを最初から最後までなめるのは以下。 ary = [\"hoge1\",\"hoge2\",\"hoge3\"] for x in range(len(ary)): print(ary[x]) hoge1 hoge2 hoge3 break,continue forループ,whileループを抜ける構文。 breakでforループを抜けた場合、for,whileに対応するelse:ブロックは評価されない。 words = [1,2,3,4,5] for x in words: print(x) if x==3: break else: print(\"finish\") 1 2 3 pass, ... pass,または...と書くと何もしない行を書くことができる。 コードに対称性が無くなったときにあえて書いておきたくなりそう。 これは悪くないかも.. words = [1,2,3,4,5] for x in words: print(x) if x==3: ... # don\'t forget elif x==4: pass # don\'t forget else: print(\"finish\") 1 2 3 4 5 finish 関数 関数の書き方。大方の言語のようにドキュメンテーションのやり方がある。 何故か関数の中に書く。docstringという。不思議。 後で自動集計してくれる。他のフォーマットを知っておきたい。 形式言語処理っぽい書き方だと、defの実行により新しいシンボル表が作られる。 関数内のあらゆる代入は新しいシンボル表に書かれる。(ローカルスコープ)。 関数内ではこのシンボル表しか参照することができず、implicitにグローバル変数を 参照できない。ほぅほぅ。 ローカルスコープからは、global識別子によりシンボルを修飾することで初めてグローバルの シンボル表にアクセスできる。 実引数がコピーされてローカルシンボル表に加えられるのではない。 関数の実引数の参照がローカルシンボル表に加えられる。 def fib(n): \"\"\"nまでのフィボナッチ級数\"\"\" a,b = 0,1 while a < n: print(a, end=' ') a,b = b, a+b print() fib(100) 0 1 1 2 3 5 8 13 21 34 55 89 関数名は1つのシンボルとして使える。以下のように f = fib f(100) 0 1 1 2 3 5 8 13 21 34 55 89 関数は\"return 値\"により値を返す。\"return\"またはreturnを書かない場合,Noneを返す。 関数のデフォルト値 大方の言語と同じ書き方でデフォルト値を書ける。 def fib(n=100): \"\"\"nまでのフィボナッチ級数\"\"\" a,b = 0,1 while a < n: print(a, end=' ') a,b = b, a+b print() fib() 0 1 1 2 3 5 8 13 21 34 55 89 関数のデフォルト値は1回しか評価されない。最初に1度だけ評価された後使い回される。 以下は\"[1],[2],[3]\"とならない。 def hoge(x,L=[]): L.append(x) return L [1] [1, 2] [1, 2, 3] こう書いておくと期待通りになる。 def fuga(x,L=None): if L is None: L = [] L.append(x) return L キーワード引数 大方の言語と同じ書き方でキーワード引数を書ける。 引数の位置から解放される。キーワード引数でない引数は前に持ってくる。 def fib(n=100): \"\"\"nまでのフィボナッチ級数\"\"\" a,b = 0,1 while a < n: print(a, end=' ') a,b = b, a+b print() fib(n=200) 0 1 1 2 3 5 8 13 21 34 55 89 144

[mathjax] マルコフ連鎖モンテカルロ法。2変量正規分布からGibbs Samplingする方法を考えてみました。 式を流してもよくわからないので、行間ゼロで理解できるまで細切れにして書いてみます。 Gibbs Sampling 無茶苦茶一般的に書かれている書籍だと以下みたいになっている。 ステップごとに(theta=(theta_1,theta_2,cdots,theta_n))の各要素を、 その要素以外の要素の条件付き確率でランダムに発生させる。 begin{eqnarray} theta^0 &=& (theta_1^0,theta_2^0,cdots,theta_n^0) \\ theta_i^{t+1} &sim& p(theta_i|theta_1^{t+1},cdots,theta_{i-1}^{t+1},theta_{i+1}^t,cdots,theta_n^t) end{eqnarray} 2次元空間であれば、ある点を決める際に、 片方の変数(theta_1)を固定して条件付き確率(P(theta_2|theta_1))を最大にするパラメタ(hat{theta_2})を決める。 その時点で((theta_1,hat{theta_2}))が決まる。 次は変数(hat{theta_2})を固定して条件付き確率(P(theta_1|hat{theta_2}))を最大にするパラメタ(hat{theta_1})を決める。 2次元正規分布であれば、片方の確率変数を固定したときの条件付き確率が1変数の正規分布になるから、 これがすごくやりやすい。 Gibbs Samplingは、このように条件付き確率を計算できる必要がある。 2変量正規分布の条件付き確率 2変量正規分布の片方の確率変数を固定すると1変数の正規分布が出てきます。 山の輪切りにして出てきた正規分布の母数を調べたいのですが、導出が無茶苦茶大変そうです。 出てきた正規分布の母数について式変形すると出てくるようです。こちらを参考にさせて頂きました。 この関係式を利用してひたすら式変形していきます。 begin{eqnarray} f(x|y) &=& frac{f(x,y)}{f(x)} end{eqnarray} 今、確率変数(y=y\')を固定し確率変数(x)の正規分布を考えます。 2変量正規分布の分散共分散行列が(sum)であるとします。 begin{eqnarray} sum &=& begin{pmatrix} sigma_{xx} & sigma_{xy} \\ sigma_{yx} & sigma_{yy} end{pmatrix} end{eqnarray} 出てくる正規分布の平均は以下となります。 確率変数が(x)なのに固定した(y)が出てくるのがポイント。 begin{eqnarray} mu\' = bar{x} + frac{sigma_{xy}}{sigma_{yy}} (y\'-bar{y}) end{eqnarray} また、出てくる正規分布の分散は以下となります。 begin{eqnarray} sigma\'^2 &=& sigma_{xx} - frac{sigma_{xy}}{sigma_{yy}} sigma_{yx} end{eqnarray} Python実装例 実際にPythonコードにしてみた図です。2変量正規分布の条件付き確率さえわかってしまえば あとはコードに落とすだけです。 import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats n_dim = 2 def gibbs_sampling(mu, sigma, sample_size): samples = [] start = [0, 0] samples.append(start) search_dim = 0 for i in range(sample_size): search_dim = 0 if search_dim == n_dim-1 else search_dim + 1 prev_sample = samples[-1][:] s11 = sigma[search_dim][search_dim] s12 = sigma[search_dim][search_dim - 1] s21 = sigma[search_dim -1 ][search_dim] s22 = sigma[search_dim - 1][search_dim - 1] mu_x = mu[search_dim] mu_y = mu[search_dim-1] _y = prev_sample[search_dim - 1] new_mean = mu_x + s12/float(s22) * (_y - mu_y) new_sigma = s11 - s12/float(s22) * s21 sample_x = np.random.normal(loc=new_mean, scale=np.power(new_sigma, .5), size=1) prev_sample[search_dim] = sample_x[0] samples.append(prev_sample) return np.array(samples) nu = np.ones(2) covariance = np.array([[0.5, 0.5], [0.5, 3]]) samples = gibbs_sampling(nu, covariance, 1000) fig, ax1 = plt.subplots(figsize=(6, 6)) ax1.scatter(sample[:, 0], sample[:, 1], marker=\"o\", facecolor=\"none\", alpha=1., s=30., edgecolor=\"C0\", label=\"Samples\" ) 実行結果

[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} 次回、微小体積とヤコビアンが等しい理由を書いてみます。