記事・メモ一覧

[mathjax] 深層学習入門。Python,Numpyにも少し慣れてきたので、 Numpyだけで伝搬,逆伝搬を計算することで深層学習に慣れていく。 単に自分の理解のためだけの記事なので、誤りがあっても気にしない。 活性化関数 まず活性化関数。activation function。 3つ(x_1,x_2,x_3)の入力があったとする。それぞれに重み(w_1,w_2,w_3)がかかるとする。 ノードは入力に重みをかけた和、つまり(w_1x_1 + w_2x_2 + w_3x_3)を受けるものとする。 ノードは受けた値の大きさに応じて出力を返す機能を持つ。 受けた値と出力の関係（まさに関数）を活性化関数と言って、 ステップ関数,Sigmoid関数,ReLU,Softmax関数などいくつか種類がある。 Step関数,Sigmoid関数,ReLUをNumpyだけで実装してみる。 import numpy as np import matplotlib.pyplot as plt # sigmoid function def sigmoid(x): return 1 / (1+np.exp(-x)) # step function def step(x): return np.array(x > 0, dtype=int) # ReLU def relu(x): return np.maximum(0,x) / 5 x = np.arange(-5.0,5.0,0.1) y = sigmoid(x) x = np.arange(-5.0,5.0,0.1) y1 = sigmoid(x) y2 = step(x) y3 = relu(x) plt.plot(x,y1) plt.plot(x,y2) plt.plot(x,y3) plt.ylim(-0.1,1.1) plt.show() どんなに入力信号が大きくても、出力を0から1の間に押し込める。 入力が大きければ出力が大きいという意図はあるものの、 入力と出力の関係が非線形になっているものが多い。 活性化関数が線形関数だと、ネットワークを重ねていったとしても、 ネットワークの最初の入り口と最後の出口を見たとして、1つの線形関数を通しただけ... ということになる。つまり加算と定数倍は何回実行したとしても、 一つの定数倍、加算の計算にまとめることができる。 これだと、層を重ねる意味がなくなってしまう。 活性化関数が非線形だと、重ねた活性化関数をまとめることはできず、 複雑な入出力を表現できるため、活性化関数として非線形関数を使用する。

あまりこういうことは書かないのだけれども、あんまりだったので記事にしておく。 FacebookやTwitterのシェア数、フォロー数などをバックグラウンドで取得するWordPressプラグイン \"SNS Count Cache\"。 なんか大量にエラーログを吐くので調べてみたら、 WP_DEBUGが立っているとinfoレベルのログをerror_log()で吐く仕様...。 開発環境でエラーログ確認しないのだろうか...。 class SCC_Logger { /** * Class constarctor * Hook onto all of the actions and filters needed by the plugin. */ protected function __construct() { self::log( \'[\' . __METHOD__ . \'] (line=\' . __LINE__ . \')\' ); } /** * Output log message according to WP_DEBUG setting * * @param string $message Message. * @return void */ public static function log( $message ) { if ( WP_DEBUG === true ) { if ( is_array( $message ) || is_object( $message ) ) { error_log( print_r( $message, true ) ); } else { error_log( $message ); } } } } log()というメソッド名でerror_log()を呼ぶというのは想像の斜め上で、 デバッグ初手のアイデアとして浮かばなかった。 /** * Class constarctor * Hook onto all of the actions and filters needed by the plugin. */ private function __construct() { SCC_Logger::log( \'[\' . __METHOD__ . \'] (line=\' . __LINE__ . \')\' ); load_plugin_textdomain( self::DOMAIN, false, basename( dirname( __FILE__ ) ) . \'/languages\' ); add_action( \'init\', array( $this, \'initialize\' ) ); register_activation_hook( __FILE__, array( $this, \'activate_plugin\' ) ); register_deactivation_hook( __FILE__, array( $this, \'deactivate_plugin\' ) ); add_action( \'admin_menu\', array( $this, \'register_admin_menu\' ) ); add_action( \'admin_print_styles\', array( $this, \'register_admin_styles\' ) ); add_action( \'admin_enqueue_scripts\', array( $this, \'register_admin_scripts\' ) ); // add_action( \'admin_notices\', array( $this, \'notice_page\' ) ); add_action( \'wp_ajax_\' . $this->ajax_action, array( $this, \'get_cache_info\' ) ); add_action( \'wp_dashboard_setup\', array( $this, \'add_wp_dashboard_widget\' ) ); add_action( \'deleted_post\' , array( $this, \'clear_cache_deleted_post\' ) ); add_filter( \'plugin_action_links_\' . plugin_basename( __FILE__ ), array( $this, \'add_plugin_action_links\' ), 10, 4 ); }

引き続きPythonデータ構造。 入れ子のリスト内包 入れ子のリスト内包を説明するために転置行列が使われていてわかりやすい。 このあたりがデータ処理用言語である所以な気がする。 matrix = [ [1,2,3,4], [5,6,7,8], [9,10,11,12] ] # 全行を出力 for i in range(4): for x in matrix: print(x) [1, 2, 3, 4] [5, 6, 7, 8] [9, 10, 11, 12] # 全行のi列目をリスト化 for i in range(4): v = [row[i] for row in matrix] print(v) [1, 5, 9] [2, 6, 10] [3, 7, 11] [4, 8, 12] # 同じことをワンライナーで v = [[row[i] for row in matrix] for i in range(4)] print(v) [1, 5, 9] [2, 6, 10] [3, 7, 11] [4, 8, 12] リストの要素削除 del文で要素を削除する。 a = [1,2,3,4,5] del a[3] print(a) [1, 2, 3, 5] del a[:] print(a) [] タプル 似て非なるデータ構造リストとタプル。リストは変更可能。タプルは変更不能。 タプルの要素として変更可能な変数を含めることはできる。 作り方が特殊。要素をカンマ区切りで並べるとタプルができる。 要素が1個の場合は、他の変数と区別が付かないので末尾にカンマを配置する。 要素が0個の場合は、()とする。 empty = () print(empty) () tupple = 100,200,300 print(tupple) (100, 200, 300) onetupple = 100, print(onetupple) (100,) 集合 データ分析用言語なのでこういうのは豊富。順序無し、重複無しのデータを格納する。 重複データを加えるときに重複が除去される。 hoge = {\'hoge1\',\'hoge2\',\'hoge3\',\'hoge4\'} print(hoge) set([\'hoge4\', \'hoge1\', \'hoge2\', \'hoge3\']) fuga = {\'fuga1\',\'fuga1\',\'fuga1\'} print(fuga) set([\'fuga1\']) リスト内包と同じ書き方で集合内包を記述できる。 v = {x for x in \'abcdefg\' if x not in \'abc\'} set([\'e\', \'d\', \'g\', \'f\']) ディクショナリ 連想配列。key-valueを格納する。 初期化の方法が複数あるので流してみる。 dic = dict([(\'hoge\',100),(\'fuga\',500),(\'hogehog\',1000),(\'fugafuga\',2000)]) print(dic) {\'fugafuga\': 2000, \'fuga\': 500, \'hogehog\': 1000, \'hoge\': 100} dic2 = dict(hoge=100, fuga=500, hogehgoe=1000, fugafuga=2000) print(dic2) {\'fugafuga\': 2000, \'fuga\': 500, \'hogehgoe\': 1000, \'hoge\': 100} ループ内でディクショナリのKeyValueを一度に取る方法を試してみる。 for k,v in dic.items(): print(k,v) (\'fuga\', 500) (\'hogehog\', 1000) (\'hoge\', 100) ちなみに、普通のリストは先頭から0,1,2,.... のように順序数が振られているはずで、 それを取得するにはenumerate()関数を使う。ループ外で変数を確保してループ内でインクリメントしたりしない。 l = [\'hoge\',\'fuga\',\'hogehoge\',\'fugafuga\'] (0, \'hoge\') (1, \'fuga\') (2, \'hogehoge\') (3, \'fugafuga\') 2つ以上のシーケンスにループをかけるとき普通はn重ループを書くと思うけれども、 Pythonはzip()関数とアンパッキングを使って1行で書くことができる。 これは良くみるし、Pythonのキモである気がする。 qq = [\'aaa1\', \'bbb1\',\'ccc1\',\'ddd1\'] rr = [100, 200, 300, 400] for q, r in zip(qq,rr): print(\'The value of qq is {0}. The value of rr is {1}\'.format(q,r)) The value of qq is aaa1. The value of rr is 100 The value of qq is bbb1. The value of rr is 200 The value of qq is ccc1. The value of rr is 300 The value of qq is ddd1. The value of rr is 400 逆順,ソート後ループ 正順のシーケンスにreversed()関数をかますと非破壊で逆順ソートできる。 それをループに使うと逆順ループができる。 sorted()関数をかますと非破壊でシーケンスをソートできる。 それをループに使うとソート後ループができる。 n = [1,2,3,4,5,6] for i in reversed(n): print(i) 6 5 4 3 2 1 r = [\'alpha\',\'beta\',\'gamma\',\'omega\',\'epsilon\'] for i in sorted(r): print(i) alpha beta epsilon gamma omega 演算子の優先順位 比較演算子は全て同じ優先順位。以下は左から順に比較が行われる。 まずa<bの比較が行われ、次にb==cの比較が行われる。a<bかつb==cであれば以下は真。 a < b == c 不思議な感じだけれども、これにより以下の比較がワンライナーでできる。 a > 1 and a <3 1 < a < 3 ブール演算子は比較演算子よりも優先順位が低い。 and,orの優先順位は同じ。notはand,orより高い。以下は同じ。 A and not B or C (A and (not B)) or C ブール演算は左から順番に評価され、途中で結論が確定したときに評価は中断される。 Cと異なり、式の途中で代入することはできない。つまり以下みたいなことはできない。 if a=10<20: シーケンスの比較 同じシーケンス型の変数を比較できる。辞書順で要素を比較していき全てが同じであれば真を返す。 両者が異なっていれば、異なっていた箇所の要素の辞書順の大小を返す。 いずれかが片方のサブセットである場合、短い方が小、長い方が大となる。

リスト 覚えるものでもないのだけれども一度は通しておきたいリストのメソッド。 大方の言語と違って、基本的に破壊的メソッドだらけ。 append()で末尾に追加。 v = [] v.append(10) v[len(v):] = [30] print(v) [10 20] リストを伸長。 r = [20] v.extend(r) print(v) [10 20 30] リストに挿入。 v.insert(0,100) print(v) [100, 10, 30, 20] 要素の削除。 v.remove(10) print(v) [100, 30, 20] 要素のインデックス。indexOf(...)みたいな。 要素が無いと-1を返しそうだけれどもエラー。 idx = v.index(20) print(idx) idx2 = v.index(1000) print(idx2) ValueError: 1000 is not in list 含まれる要素xの個数を返す不思議なメソッド。 cnt = v.count(20) print(cnt) 要素を逆順にする。コピーしたものを逆順にするのではなくそのものを逆順にする。 v.reverse() [20, 30, 100] コピーする。 vv = v.copy() vv[0] = \'hoge\' print(v) print(vv) [20, 30, 100] [\'hoge\', 30, 100] リスト内包 例えば0から9までの値の2乗のリストを作成するコードは以下。 squares = [] for x in range(10): squares.append(x**2) [0, 1, 4, 9, 16, 25, 36, 49, 64, 81] Pythonっぽくワンライナーで書くとこうなる。 map関数は高階関数。第1引数として関数、第2引数としてシーケンスを取る。 lambdaは関数のSyntaxSugarなので第1引数はlambdaも渡せる。 渡したシーケンスの全要素についてlambda式を実行した結果を返す。 range(10)は0から9まで。lambda x:x**2は引数を2乗する。 なのでmap(lambda x: x**2, range(10))は、0 1,4,9,...を返す。 それをlist()でリスト化する。 squares = list(map(lambda x: x**2, range(10))) print(squares) [0, 1, 4, 9, 16, 25, 36, 49, 64, 81] Pythonでは、これが言語仕様で出来るようになっている。リスト内包。 []の中の最初は式。式に続けてfor式,if式を書ける。 squares2 = [x**2 for x in range(10)] print(squares2) [0, 1, 4, 9, 16, 25, 36, 49, 64, 81] より複雑な式,for式をリスト内包にしてみる。以下は等価。 cs = [] for x in [1,2,3]: for y in [3,1,4]: if x != y: cs.append((x,y)) print(cs) [(1, 3), (1, 4), (2, 3), (2, 1), (2, 4), (3, 1), (3, 4)] cs2 = [(x,y) for x in [1,2,3] for y in [3,1,4] if x != y] print(cs2) [(1, 3), (1, 4), (2, 3), (2, 1), (2, 4), (3, 1), (3, 4)] 色々な書き方が出来る。戻り値をタプルにしたい場合は()で囲う必要がある。 vec = [-4, -2, 0 ,2, 4] cs3 = [x*2 for x in vec] print(cs3) [-8, -4, 0, 4, 8] cs4 = [x for x in vec if x >= 0] print(cs4) [0, 2, 4] cs5 = [(x, x**2) for x in range(10)] print(cs5) [(0, 0), (1, 1), (2, 4), (3, 9), (4, 16), (5, 25), (6, 36), (7, 49), (8, 64), (9, 81)] 以下面白い。リスト内包には入れ子の関数も含められる。 from math import pi cs6 = [str(round(pi, i)) for i in range(1,6)] print(cs6)

任意引数 大方の言語と同じように任意引数を定義できる。任意引数の位置に渡した引数はタプルに変換される。 複数の変数を1つのタプルに変換する操作を、詰める(pack)と言ったりする。 もちろん任意引数の後はキーワード引数しか置けない。 def concat(*args,sep=\"/\"): return sep.join(args) val=concat(\'hoge\',\'fuga\',\'hogehoge\',\'fugafuga\') print(val) hoge/fuga/hogehoge/fugafuga unpack 逆に1つのタプルなりリストを複数の変数にバラすことをunpackと言ったりする。 hoge = [1,2,3] x,y,z = hoge print(x) print(y) print(z) 1 2 3 lambda 大方の言語と同じように無名関数を書ける。lambdaは1行しか書けない。 中では関数になっている。関数が関数外の変数を参照できるのと同様に lambdaからlambda外の変数を参照できる。 def make_incrementor(n): return lambda x:x+n f = make_incrementor(42) v =f(0) print(v) 42 r =f(10) print(r) 52 関数アノテーション 関数の引数、戻り値に付加情報を付けることができる。 PHPのTypeHintingみたいに、引数の型,戻り値の型を書ける。 型として実在しないクラスを指定するとエラーになる。 def f(hoge: str, fuga: str = \'fuga\') -> str: print(\"Annotations:\", f.__annotations__) print(\"Arguments:\", hoge, fuga) return hoge + \' and \' + fuga f(\'1\',\'2\') Annotations: {\'hoge\': , \'fuga\': , \'return\': } Arguments: 1 2 コーディングスタイル Pythonのコーディングスタイル。なんだかその理由がどうでも良い感じ。 やるのであればRubyのように外れたらエラーにするべきだし。 インデントはスペース4個。タブは使わない。 79文字以下で折り返す。 関数、クラス、関数内の大きめのブロックを分離するために空白行を使用 コメントは文末ではなく行として独立すべき docstringを使う 演算子の周囲やカンマの後ろはスペースをいれる。 括弧のすぐ内側にはスペースをいれない。 関数名、メソッドはlower_case_with_underscores UTF-8を使う。

個人的には、言語に入門する際にはオンラインより書籍が適切と考えている。 そして最初の一冊はなるべく薄いものを選ぶべきと考えている。 読んだ後、あえて文章として書き出すと良さそう。行間なりツボを自分の言葉でおまとめ。 あそこにこう書いてあった..というように思い出せると良いかなと。 まずは制御構造。 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);