MNISTの手書き文字認識用データ取得クラスの作成
MNISTから手書き文字認識用のデータセットをロードするクラスを作ってみた。 ロードしたデータセットをpickleでシリアライズ、デシリアライズする機能付き。 後々改造する予定でここに貼ったのはメンテしない予定。 たったこれだけ書くのに40分も要してしまった...。 Pythonの冗長感が半端ないけども慣れるしかない。 マジックコードだらけだけども、MNISTの手書き文字認識用データ取得専用だから仕方ない。 get_image()関数により訓練データ、訓練ラベル、テストデータ、テストラベルを取れる。 データは1次元になって入る。つまり1ファイルごとに28*28のサイズがある。 import os import urllib.request import numpy as np import gzip import pickle class MnistLoader: def __init__(self, mnistdir): self.url_base = \'http://yann.lecun.com/exdb/mnist/\' self.dataset_dir = mnistdir self.save_path = self.dataset_dir + \'mnist.pkl\' self.key_file = { \'train_img\':\'train-images-idx3-ubyte.gz\', \'train_label\':\'train-labels-idx1-ubyte.gz\', \'test_img\':\'t10k-images-idx3-ubyte.gz\', \'test_label\':\'t10k-labels-idx1-ubyte.gz\' } self.dataset = {} if os.path.exists(self.save_path): with open(self.save_path,\'rb\') as f: self.dataset = pickle.load(f) else: self.__load_mnist() self.dataset[\'train_img\'] = self.__load_img(self.key_file[\'train_img\']) self.dataset[\'train_label\'] = self.__load_label(self.key_file[\'train_label\']) self.dataset[\'test_img\'] = self.__load_img(self.key_file[\'test_img\']) self.dataset[\'test_label\'] = self.__load_label(self.key_file[\'test_label\']) with open(self.save_path, \'wb\') as f: pickle.dump(dataset, f, -1) def __load_mnist(self): \'\'\' load mnist data and store to file \'\'\' for v in self.key_file.values(): file_path = self.dataset_dir + v urllib.request.urlretrieve(url_base + v, self.dataset_dir + v) def __load_img(self, file_name): file_path = self.dataset_dir + file_name with gzip.open(file_path, \'rb\') as f: data = np.frombuffer(f.read(), np.uint8, offset=16) data = data.reshape(-1, 784) return data def __load_label(self,file_name): file_path = dataset_dir + \'/\' + file_name with gzip.open(file_path, \'rb\') as f: labels = np.frombuffer(f.read(), np.uint8, offset=8) return labels def get_image(self): return (self.dataset[\'train_img\'] , self.dataset[\'train_label\']), (self.dataset[\'test_img\'], self.dataset[\'test_label\']) mnist_loader = MnistLoader(mnistdir=\'/Users/ikuty/Documents/mnist/\') PILを使って訓練データの1枚目を表示するテスト。 import sys, os sys.path.append(os.pardir) import numpy as np from PIL import Image def img_show(img): pil_img = Image.fromarray(np.uint8(img)) pil_img.show() (x_train, t_train), (x_test, t_test) = mnist_loader.get_image() img = x_train[0] label = t_train[0] bimg = img.reshape(28, 28) img_show(bimg) こんなのが出る。(28x28しかなくて小さすぎなので256x256に引き伸ばして表示。)
SoftMax関数
[mathjax] 出力層の総和を1にするように調整できれば、出力層を確率としてとらえることができるようになる。 入力層に画像を放り込み、出力層でラベルに属する確率を出せば、画像にラベルをつける分類器になる。 入力が(n)ノードあるとして、全てのノードを入力を受けて出力する。 出力側のノードが(m)ノードあるとして、和が(1)になるようにどう操作するか。 分母は入力側の全ノードの和、分子は対応する1ノードだけ。そうやって1にする。 指数関数に入力を食わせて和を取る。指数関数は単調増加関数だから、 入力の大小と出力の大小は一致する。定義上は以下の通り。 begin{eqnarray} y_k &=& frac{e^{a_k}}{sum_{i=1}^n e^{a_i}} end{eqnarray} import numpy as np def softmax1(x): return np.exp(x) / np.sum(np.exp(x)) a = np.array([1010,1000,990]) r1 = softmax1(a) r1 # array([nan, nan, nan]) ただ、この通りにするとオーバーフローを起こすので、 式変形してオーバーフローを回避する。(C)は入力ノードの最大値。 begin{eqnarray} y_k &=& frac{e^{a_k}}{sum_{i=1}^n e^{a_i}} \\ &=& frac{C e^{a_k}}{Csum_{i=1}^n e^{a_i}} \\ &=& frac{e^{a_k+log{C}}}{sum_{i=1}^{n} e^{a_i + log{C}}} \\ &=& frac{e^{a_k+C\'}}{sum_{i=1}^{n} e^{a_i + C\'}} end{eqnarray} import numpy as np def softmax2(x): c = np.max(x) return np.exp(x-c) / np.sum(np.exp(x-c)) a = np.array([1010,1000,990]) r2 = softmax2(a) r2 # array([9.99954600e-01, 4.53978686e-05, 2.06106005e-09]) SoftMax後の和は1 前述の通り、SoftMax後の和は1。 z = np.sum(r2) z # 1.0
行列の積和と順伝播
[mathjax] 順伝搬(forward)はNumpyの積和計算を使って超絶簡単に記述できる。 行列の積和計算 何はともあれNumpyで行列の積を計算する方法について。 この仕組みがあるから、forなどの制御構造を使わずに行列の積を実現できる。 begin{eqnarray} A &=& begin{pmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} end{pmatrix} \\ B &=& begin{pmatrix} b_{11} & b_{21} \\ b_{12} & b_{22} end{pmatrix} \\ A cdot B &=& begin{pmatrix} a_{11} b_{11} + a_{21} b_{12} && a_{11} b_{21} + a_{21} b_{22} \\ a_{12} b_{11} + a_{22} b_{12} && a_{12} b_{21} + a_{22} b_{22} end{pmatrix} end{eqnarray} 積はNumpyのdot関数で一発。 左行列と右行列の\"対応する行\"は一致している必要がある。 つまり(N_1 times M_1 cdot N_2 times M_2)という積計算において(M_1 = N_2)である必要がある。 a11,a21,a12,a22 = (1,2,3,4) b11,b21,b12,b22 = (1,2,3,4) a = np.array([[a11,a21],[a12,a22]]) b = np.array([[b11,b21],[b12,b22]]) np.dot(a,b) # array([[ 7, 10], # [15, 22]]) 和はNumpyの+オペレータで一発。 begin{eqnarray} A &=& begin{pmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} end{pmatrix} C &=& begin{pmatrix} c_{1} \\ c_{2} end{pmatrix} end{eqnarray} begin{eqnarray} A+C &=& begin{pmatrix} a_{11} + c_1 & a_{21} + c_1 \\ a_{12} + c_1 & a_{22} + c_2 end{pmatrix} end{eqnarray} c1, c2 = (5,5) c = np.array([[c1],[c2]]) a + c # array([[6, 7], # [8, 9]]) 各層における信号伝達の実装 Numpyで積和が簡単にかける。 すなわち、各層における伝達が積和で表されるとすると、 入力層から各層における信号伝達の実装を簡単にかける。 左側の層が(X=begin{pmatrix}x_1 x_2end{pmatrix})、重みが(W=begin{pmatrix}w_{11} & w_{21} \\ w_{12} & w_{22} end{pmatrix})、 右側の層が(A=begin{pmatrix} a_{11} & a_{21} end{pmatrix})とする。 すると、右層と左層の関係は以下の通り。 begin{eqnarray} A &=& XW \\ begin{pmatrix} a_{11} & a_{21} end{pmatrix} &=& begin{pmatrix}x_1 x_2end{pmatrix} begin{pmatrix}w_{11} & w_{21} \\ w_{12} & w_{22} end{pmatrix} \\ &=& begin{pmatrix} w_{11}x_1 + w_{21} x_2 \\ w_{12}x_1 + w_{22}x_2end{pmatrix} end{eqnarray} 左層にバイアス(B=begin{pmatrix}b_1 \\ b_1 end{pmatrix})があった場合、線形和で表現できる。 行と列がどちらの方向に対応するかに注意。 1つの層に(N)個のノードが存在することを(Ntimes 1)の行列で表現する。 begin{eqnarray} A &=& XW + B \\ begin{pmatrix} a_{11} & a_{21} end{pmatrix} &=& begin{pmatrix}x_1 x_2end{pmatrix} begin{pmatrix}w_{11} & w_{21} \\ w_{12} & w_{22} end{pmatrix} + begin{pmatrix}b_1 \\b_1 end{pmatrix} \\ &=& begin{pmatrix} w_{11}x_1 + w_{21} x_2 + b_1 \\ w_{12}x_1 + w_{22}x_2 + b_1end{pmatrix} end{eqnarray} x1,x2 = (0.4,0.9) w11,w21,w12,w22 = (0.8,0.5,0.1,0.3) b1 = 0.1 x = np.array([x1,x2]) w = np.array([[w11,w21],[w12,w22]]) b = np.array([b1,b1]) # これで一発 a = np.dot(x,w) + b # array([0.51, 0.57]) # 活性化関数としてsigmoid関数をかます z = sigmoid(a) # array([0.62480647, 0.63876318]) 順伝播 入力層から出力層まで上記のような計算を繰り返していく。これを\"順伝播\"とか言う。 順伝\"搬\"ではなく、順伝\"播\"。ネットワークが多層になったとしても全く同様。 以下みたいにNumpyの積和計算を繰り返すだけで出来る。 def init_network(): network = {} network[\'W1\'] = np.array([[0.1,0.3,0.5],[0.2,0.4,0.6]]) network[\'b1\'] = np.array([0.1,0.2,0.3]) network[\'W2\'] = np.array([[0.1,0.4],[0.2,0.5],[0.3,0.6]]) network[\'b2\'] = np.array([0.1,0.2]) network[\'W3\'] = np.array([[0.1,0.3],[0.2,0.4]]) network[\'b3\'] = np.array([0.1,0.2]) return network def forward(network, x): W1, W2, W3 = network[\'W1\'], network[\'W2\'], network[\'W3\'] b1,b2,b3 = network[\'b1\'], network[\'b2\'], network[\'b3\'] a1 = np.dot(x, W1) + b1 z1 = sigmoid(a1) a2 = np.dot(z1, W2) + b2 z2 = sigmoid(a2) a3 = np.dot(z2, W3) + b3 y = a3 return y network = init_network() x = np.array([1.0,0.5]) y = forward(network,x) # array([0.31682708, 0.69627909])
活性化関数の実装。Step,Sigmoid,ReLU
[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つの線形関数を通しただけ... ということになる。つまり加算と定数倍は何回実行したとしても、 一つの定数倍、加算の計算にまとめることができる。 これだと、層を重ねる意味がなくなってしまう。 活性化関数が非線形だと、重ねた活性化関数をまとめることはできず、 複雑な入出力を表現できるため、活性化関数として非線形関数を使用する。
SNS Count Cache… WP_DEBUG=TRUEでinfoをerrorに吐くのは仕様です
あまりこういうことは書かないのだけれども、あんまりだったので記事にしておく。 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データ構造2
引き続き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: シーケンスの比較 同じシーケンス型の変数を比較できる。辞書順で要素を比較していき全てが同じであれば真を返す。 両者が異なっていれば、異なっていた箇所の要素の辞書順の大小を返す。 いずれかが片方のサブセットである場合、短い方が小、長い方が大となる。
Pythonデータ構造
リスト 覚えるものでもないのだけれども一度は通しておきたいリストのメソッド。 大方の言語と違って、基本的に破壊的メソッドだらけ。 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)
Python制御構造2
任意引数 大方の言語と同じように任意引数を定義できる。任意引数の位置に渡した引数はタプルに変換される。 複数の変数を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を使う。
python制御構造1
個人的には、言語に入門する際にはオンラインより書籍が適切と考えている。 そして最初の一冊はなるべく薄いものを選ぶべきと考えている。 読んだ後、あえて文章として書き出すと良さそう。行間なりツボを自分の言葉でおまとめ。 あそこにこう書いてあった..というように思い出せると良いかなと。 まずは制御構造。 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
やってみた Markov chain Monte Carlo methods, MCMC , gibbs sampling
[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\" ) 実行結果