default eye-catch image.

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)

default eye-catch image.

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を使う。

default eye-catch image.

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

default eye-catch image.

PHP7.4の新仕様 … PHP RFC: FFI – Foreign Function Interface

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);

default eye-catch image.

PHPで統計アプリを作れるか否か

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のクラスを使うとか...

default eye-catch image.

多次元正規分布に従うデータを生成する

[mathjax] そろそろ適当なデータを見つけてきて手法を試すのとは別に、 自力でデータを作って試してみたいと思い、NumPyを使った生成法を調べてみた。 一口に乱数といっても、正規分布に従う標本の生成のこと。 多次元正規分布に従う標本をmultivariate_normalで生成して表示してみる。 1次元正規分布に従う標本 その前に普通の乱数。 平均(mu)、標準偏差(sigma)の正規分布(N(mu,sigma))に従う標本の生成。 numpy.random.normal((mu,sigma,n))。 以下、(mu=50,sigma=10)として1000個の標本を作って、 階級数100のヒストグラムとして表示する。 import numpy as np import matplotlib.pyplot as plt values = np.random.normal(50, 10,1000) plt .hist(values,bins=100) 多次元正規分布に従う標本 多次元の乱数を作るのに必要なのは、 各次元における平均(mu=begin{pmatrix}mu_1 \\ mu_2 end{pmatrix})と、分散共分散行列(sum=begin{pmatrix}sigma_{11} & sigma_{21} \\ sigma_{12} & sigma_{22}end{pmatrix})。 ちなみに、(sum=begin{pmatrix}sigma_{11} & sigma_{21} \\ sigma_{12} & sigma_{22}end{pmatrix} = begin{pmatrix} sigma_1^2 & sigma_{21} \\ sigma_{12} & sigma_2^2 end{pmatrix} = begin{pmatrix} sigma_1^2 & Cov(X_1,X_2) \\ Cov(X_2,X_1) & sigma_2^2 end{pmatrix} )。 対角に分散、それ以外にクロスする共分散が入る。 共分散は昔書いてた。 begin{eqnarray} r_{xy} &=& frac{C_{xy}}{S_x S_y} \\ &=& frac{sum_{i=1}^n(x_i-bar{x})(y_i-bar{y})/n}{sqrt{sum_{i=1}^n{(x_i-bar{x})^2}/n} sqrt{sum_{i=1}^n{(y_i-bar{y})^2}/n}} \\ &=& frac{sum_{i=1}^n(x_i-bar{x})(y_i-bar{y})}{sqrt{sum_{i=1}^n{(x_i-bar{x})^2}} sqrt{sum_{i=1}^n{(y_i-bar{y})^2}}} \\ end{eqnarray} [clink url=\"https://ikuty.com/2018/10/17/covariance_zero/\"] 平均(mu=begin{pmatrix} 0 \\ 0end{pmatrix})、分散共分散行列(sum=begin{pmatrix}30 & 20 \\ 20 & 50end{pmatrix})の2次元正規分布に従う標本を1000個作る。 import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [0,0] sigma = [[30,20],[20,50]] values = np.random.multivariate_normal(mu, sigma, 1000) sns.jointplot(x=values[:,0],y=values[:,1]) 散布図とヒストグラムを良い感じに表示してくれるseaborn.jointplotを使って表示してみる。 第1軸、第2軸とも平均は0っぽい。 第1軸は-15から15、第2軸は-25から25くらいが入っていて、対角とあってそう。

default eye-catch image.

sklearnに頼らずRidge回帰を自力で書いてみて正則化項の影響を考えてみるテスト

[mathjax] タイトルの通り。Losso回帰と違って損失関数を偏微分するだけで出来そうなのでやってみる。 Ridge回帰は線形回帰の1種だけれども、損失関数として最小二乗法をそのまま使わず、 (L_2)ノルムの制約を付けたものを使う((L_2)正則化)。 データとモデル 教師データ(boldsymbol{y})、訓練データ(boldsymbol{x})があるとする。 (または目的変数(boldsymbol{y})、説明変数(boldsymbol{x})があるとする。) 例えば(p)次の属性データが(n)個あり、それらと結果の対応が分かっている状況。 begin{eqnarray} boldsymbol{y} &=& begin{pmatrix} y_1 \\ y_2 \\ vdots \\ y_p end{pmatrix} , boldsymbol{x} &=& begin{pmatrix} x_{11} & x_{21} & cdots & x_{n1} \\ x_{12} & x_{22} & cdots & x_{n2} \\ vdots & vdots & ddots & vdots \\ x_{1p} & x_{2p} & cdots & x_{np} end{pmatrix} end{eqnarray} モデルは以下。特徴ベクトル(boldsymbol{w})は訓練データの重み。 特徴空間において損失を最小化する特徴ベクトルを求める問題。 begin{eqnarray} boldsymbol{y} &=& boldsymbol{w} boldsymbol{x} + k \\ boldsymbol{w} &=& begin{pmatrix} w_1 & w_2& cdots &w_p end{pmatrix} end{eqnarray} 損失関数 普通の2乗損失に正則化項((L_2)ノルムを定数倍した値)を付けたものを損失関数として利用する。 正則化項の係数はハイパーパラメータとして調整する値。逆数なのはsklearnに従う。 begin{eqnarray} L(boldsymbol{w}) = |boldsymbol{y} - boldsymbol{w} boldsymbol{x}|^2 +C |boldsymbol{w}|^2 end{eqnarray} 特徴ベクトルは以下。(mathjaxでargminが出せない...) begin{eqnarray} newcommand{argmin}[1]{underset{#1}{operatorname{arg},operatorname{min}};} boldsymbol{w} = argmin w L(boldsymbol{w}) = argmin w |boldsymbol{y} - boldsymbol{w} boldsymbol{x}|^2 + C |boldsymbol{w}|^2 end{eqnarray} 特徴ベクトルを求める 勾配=0と置けば上の式の解を得られる。 損失関数が微分可能だからできる技。 begin{eqnarray} frac{partial L(boldsymbol{w})}{partial boldsymbol{w}} &=& 2 boldsymbol{w}^T (boldsymbol{y} - boldsymbol{w} boldsymbol{x}) + C boldsymbol{w} \\ &=& 0 end{eqnarray} 変形する。 begin{eqnarray} 2 boldsymbol{x}^T (boldsymbol{x}boldsymbol{w}-boldsymbol{y}) + C boldsymbol{w} &=& 0 \\ boldsymbol{x}^T (boldsymbol{x}boldsymbol{w}-boldsymbol{y}) + C boldsymbol{w} &=& 0 \\ boldsymbol{x}^T boldsymbol{x} boldsymbol{w} -boldsymbol{x}^T boldsymbol{y} + Cboldsymbol{w} &=& 0 \\ (boldsymbol{x}^T boldsymbol{x} +C E) boldsymbol{w} &=& boldsymbol{x}^T boldsymbol{y} \\ boldsymbol{w} &=& (boldsymbol{x}^T boldsymbol{x} + C E)^{-1} boldsymbol{x}^T boldsymbol{y} end{eqnarray} テストデータを作る 練習用にsklearnのbostonデータを使ってみる。 ボストンの住宅価格が目的変数、属性データが説明変数として入ってる。 import pandas as pd import numpy as np from pandas import Series,DataFrame import matplotlib.pyplot as plt from sklearn.datasets import load_boston boston = load_boston() boston_df = DataFrame(boston.data) boston_df.columns = boston.feature_names print(boston_df.head()) boston_df[\"PRICE\"] = DataFrame(boston.target) # CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT PRICE # 0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0 15.3 396.90 4.98 24.0 # 1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0 17.8 396.90 9.14 21.6 # 2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0 17.8 392.83 4.03 34.7 # 3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0 18.7 394.63 2.94 33.4 # 4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 18.7 396.90 5.33 36.2 散布図行列を表示してみる。 PRICEと関係がありそうなZN,RM,AGE,DIS,LSTATの5個を使ってみる。 pg = sns.pairplot(boston_df) plt.show() pg.savefig(\'boston_fig.png\') 特徴ベクトルを自力で計算する これを自力で計算してみる。(C=0.01)、(C=0)、(C=100)としてみた。 begin{eqnarray} boldsymbol{w} &=& (boldsymbol{x}^T boldsymbol{x} + C E)^{-1} boldsymbol{x}^T boldsymbol{y} end{eqnarray} X_df = boston_df.drop(columns=[\'CRIM\',\'INDUS\',\'CHAS\',\'NOX\',\'RAD\',\'TAX\',\'PTRATIO\',\'B\',\'PRICE\']) X = X_df.values y = boston.target.T C1 = 0.01 C2 = 0 C3 = 100 e = np.identity(5) w1 = np.dot( np.linalg.inv(np.dot(X.T , X) + C1 * e), np.dot(X.T,y)) w2 = np.dot( np.linalg.inv(np.dot(X.T , X) + C2 * e), np.dot(X.T,y)) w3 = np.dot( np.linalg.inv(np.dot(X.T , X) + C3 * e), np.dot(X.T,y)) print(w1) # [ 0.05338557 5.40396159 -0.01209002 -0.83723303 -0.63725397] print(w2) # [ 0.05338539 5.40403743 -0.01209427 -0.83728837 -0.63725093] print(w3) # [ 0.05612977 4.76664789 0.02374402 -0.38576708 -0.66137596] (C=0)のとき、つまり最小二乗法のとき。 sklearnを使う sklearnのridge回帰モデルを使うと以下みたいになる。 from sklearn.linear_model import Ridge from sklearn.model_selection import train_test_split Xf_train,Xf_test,yf_train,yf_test = train_test_split(X,y,random_state=0) ridge = Ridge().fit(Xf_train,yf_train) print(f\"accuracy for training data:{ridge.score(Xf_train,yf_train):.2}\") print(f\"accuracy for test data:{ridge.score(Xf_test,yf_test):.2f}\") # accuracy for training data:0.68 # accuracy for test data:0.58 print(ridge.coef_) # [ 0.06350701 4.3073956 -0.02283312 -1.06820241 -0.73188192] 出てきた特徴ベクトルを並べてみる 自力で計算したものとsklearnに計算してもらったものを並べてみる。 似てるのか似ていないのかよくわからない .. けど、RMの寄与度が高いというのは似ている。 # 自力で計算 (C=100) # [ 0.05612977 4.76664789 0.02374402 -0.38576708 -0.66137596] # sklearnで計算 # [ 0.06350701 4.3073956 -0.02283312 -1.06820241 -0.73188192] 自力で計算したモデルの正答率を求めてみないとなんとも... そして、正規化項の係数の大小がどう影響するのか、あまり良くわからなかった..。 (L_2)ノルムの制約を付けると、パラメタの大小が滑らかになると言いたかったのだけども。 あと、訓練データに対して68%、テストデータに対して58%という感じで、 大して成績が良くない...。 

default eye-catch image.

NumPy uniqe, File I/O

集合関数 集合関数。ndarrayから重複を取り除きsortした結果を返す。 2dであってもその中から要素を抜き出して1dにする。 hoges = np.array([\"hoge\",\"fuga\",\"hoge\",\"fuga\"]) print(np.unique(hoges)) # [\'fuga\' \'hoge\'] fugas = np.array([[\"hoge\",\"fuga\",\"hoge\",\"fuga\"],[\"hoge2\",\"fuga2\",\"hoge2\",\"fuga2\"]]) print(np.unique(fugas)) # [\'fuga\' \'fuga2\' \'hoge\' \'hoge2\'] ファイルI/O pandasを使わずとも、NumPyだけでファイルI/Oができる。 以下でhoges.npyという無圧縮バイナリファイルが作られる。 それを読み込んで出力する。 hoges = np.array([\"hoge\",\"fuga\",\"hoge\",\"fuga\"]) np.save(\'hoges.npy\',hoges) fugas = np.load(\'hoges.npy\') print(fugas) # [\'hoge\' \'fuga\' \'hoge\' \'fuga\'] 複数の配列を同時に書き込むこともできる。 キーワードを指定して書き込む。キーワードを指定して1つずつ読み込む。 読み込む時はキーワードを指定して参照したときに遅延ロードされる。 hoges = np.array([\"hoge1\",\"fuga1\",\"hoge1\",\"fuga1\"]) fugas = np.array([\"hoge2\",\"fuga2\",\"hoge2\",\"fuga2\"]) np.savez(\'hogefuga.npz\', hoges=hoges, fugas=fugas) hogefugas = np.load(\'hogefuga.npz\') hoges_l = hogefugas[\'hoges\'] fugas_l = hogefugas[\'fugas\'] print(hgoes_l) # [\'hoge1\' \'fuga1\' \'hoge1\' \'fuga1\'] print(fugas_l) # [\'hoge2\' \'fuga2\' \'hoge2\' \'fuga2\']

default eye-catch image.

線形サポートベクトル分類器で画像認識するテスト

線形サポートベクトル分類器で画像認識する流れを理解したので、 定着させるために記事にしてみます。 当然、モデルの数学的な理解がないとモデルを解釈することは不可能だし、 正しいハイパーパラメータを設定することも不可能なので、数学的な理解は不可欠。 NumPy、pandas、matplotlibに慣れないと、そこまで行くのに時間がかかります。 こちらはPythonプログラミングの領域なので、数こなして慣れる他ないです。 機械学習用のサンプル画像で有名なMNISTを使ってNumPy、pandasの練習。 手書き文字認識用の画像データを読み込んでみる。サイズは28x28。各々1byte。 MNISTの手書き文字認識画像の読み込み まず読み込んでみて、データの形を出力してみる。 X_trainは、要素が3個のTupleが返る。3次。 1番外が60000。28x28の2次のndarrayが60000個入っていると読む。 1枚目の画像データはX_train[0]によりアクセスできる。 import tensorflow as tf minst = tf.keras.datasets.mnist (X_train,y_train),(X_test,y_test) = mnist.load_data() print(X_train.shape) # (60000, 28, 28) y_trainは要素が1個のTupleが返る。1次。 1枚目から60000枚目までの画像が0から9のいずれに分類されたかが入っている。 y_train[0]が4なら、1枚目の画像が4に分類された、という意味。 print(y_train.shape) # (60000,) データセットの選択 X_train,y_train、X_test,y_testから、値が5または8のものだけのViewを取得する。 そのために、まず値が5または8のものだけのインデックスを取得する。 NumPyのwhereはndarrayのうち条件を満たす要素のインデックスを返す。 X_trainに入っている60000件の2d arrayのうち、 値が5または8のインデックス(0-59999)を取得するのは以下。 index_train = np.where((X_train==5)|(X_train==8)) print(index_train) # (array([ 0, 11, 17, ..., 59995, 59997, 59999]),) index_test = np.where((X_test==5)|(X_test==8)) print(index_test) # (array([ 8, 15, 23, ..., 9988, 9991, 9998]),) インデックスを使って絞り込む。 X_train,y_train = X_train[index_train],y_train[index_train] X_test,y_test = X_test[index_test],y_test[index_test] print(X_train.shape) # (11272, 28, 28) print(X_test.shape) # (1866, 28, 28) 前処理 0-255の間の値を0-1の間の値に変換する(正規化)。 28x28の画像(2darray)を1x784(1darray)に整形する(平坦化)。 X_train,X_test = X_train / 255.0, X_test / 255.0 X_train = X_train.reshape(X_train.shape[0], X_train.shape[1] * X_train.shape[2]) X_test = X_test.reshape(X_test.shape[0], X_test.shape[1] * X_test.shape[2]) ベストなハイパーパラメータの選択 線形サポートベクトル分類器を作成する。 from sklearn.svm import LinearSVC linsvc = LinearSVC(loss=\"squared_hinge\",penalty=\"l1\",dual=False) 線形サポートベクトル分類器のハイパーパラメータCの選択 逆正則化パラメータCをGridSearchCVで探す。MBP2013Laterで学習(fit)に5分くらいかかった。 GridSearchCVからはC=0.2がbestと返ってくる。 from sklearn.model_selection import GridSearchCV param_grid = {\"C\":[0.025,0.05,0.1,0.2,0.4]} model = GridSearchCV(estimator=linsvc, param_grid=param_grid,cv=5,scoring=\"accuracy\",return_train_score=True) model.fit(X_train,y_train) print(model.cv_results_[\"mean_train_score\"]) # array([0.96291693, 0.96775192, 0.97059085, 0.97340754, 0.97626859]) print(model.cv.results_[\"mean_test_score\"]) # array([0.95626331, 0.95990064, 0.96158623, 0.9625621 , 0.96105394]) print(model.best_params_) # {\'C\': 0.2} 学習、精度評価 C=0.2を使って新しく学習させる。 linsvc = LinearSVC(loss=\"squared_hinge\",penalty=\"l1\",dual=False,C=0.2) linsvc.fit(X_train,y_train) 訓練データ、テストデータに対して正答率を求める。 訓練データについて97.2%、テストデータについて96.2%。 過学習すると訓練データが高くテストデータが低くなる。 from sklearn.metrics import accuracy_score pred_train = linsvc_best.predict(X_train) acc = accuracy_score(y_true = y_train,y_pred = pred_train) print(acc) # 0.9723207948899929 pred_test = linsvc_best.predict(X_test) acc = accuracy_score(y_true = y_test,y_pred = pred_test) print(acc) # 0.9619506966773848 モデルの解釈可能性 [mathjax] 線形SVMの決定境界(f(x))の係数をヒートマップっぽく表示して、どの係数を重要視しているかを確認する。 基本的に真ん中に画像が集まっているので、28x28の隅は使わないのが正しそう。 正則化パラメータによって係数の大きさを制御しているため、正則化パラメータを変えると係数が変わる。 今回のは(L_1)正則化なので、係数が0のものが増える..らしい(..別途調べる..)。 (f(x) = w_0 + w_1 x_1 + w_2 x_2 + cdots w_{784} x_{784}) import matplotlib.pyplot as plt weights = linsvc_best.coef_ plt.imshow(weights.reshape(28,28)) plt.colorbar() plt.show()

default eye-catch image.

NumPy vector operations, universal functions, matplotlib, 3項演算, 次元削減

universal functions ndarrayの全ての要素に対して基本的な計算を実行する。 以下オペランドが1つの単項universal functions。 abs,sqrt,square,exp,log,sign,ceil,floor,rint,modf,isnan,sin,cos,arcsin,arccosなどがある。 array = np.arange(10) print(array) # [0 1 2 3 4 5 6 7 8 9] sqrt = np.sqrt(array) print(sqrt) # [0. 1. 1.41421356 1.73205081 2. 2.23606798 # 2.44948974 2.64575131 2.82842712 3. ] exp = np.exp(array) print(exp) # [1.00000000e+00 2.71828183e+00 7.38905610e+00 2.00855369e+01 # 5.45981500e+01 1.48413159e+02 4.03428793e+02 1.09663316e+03 # 2.98095799e+03 8.10308393e+03] 以下、オペランドが2つの2項universal functions。 いずれかのうち最大の値を残すmaximum()。 add,subtract,divide,power,maximum,minimum,copysign,greater,lessなどがある。 x = np.random.randn(10) y = np.random.randn(10) print(x) # [ 1.3213258 0.12423666 -1.45665939 -1.49766467 -0.6129116 2.00056744 # -0.00816571 0.63247747 0.29497652 0.80000291] print(y) # [-0.76739214 0.95151629 0.03208859 0.40641677 0.82635027 1.01773826 # 0.75601178 0.25200147 1.59929321 0.6251983 ] z = np.maximum(x,y) print(z) # [1.3213258 0.95151629 0.03208859 0.40641677 0.82635027 2.00056744 # 0.75601178 0.63247747 1.59929321 0.80000291] [mathjax] matplotlibにndarrayを引数として渡せば簡単にプロットできる。 (z=sqrt{x^2+y^2})をプロットしてみる。 import numpy as np import matplotlib.pyplot as plt points = np.arange(-5,5,0.01) xs,ys = np.meshgrid(points,points) z = np.sqrt(xs**2 +ys**2) plt.imshow(z, cmap=plt.cm.gray) plt.colorbar() plt.title(\"Image plot\") plt.show() 3項演算子 where マスクの論理値に従って2つのndarrayのうちいずれかの値を選択してリストに書く。 3項演算子を使ってPythonのlistに入れる方法は以下。 xa,xbはndarrayだが最終的なr1はPythonオブジェクト。 import numpy as np xa = np.array([1,2,3,4,5]) xb = np.array([6,7,8,9,10]) cnd = np.array([True,True,False,False,False]) r1 = [(x if c else y) for x,y,c in zip(xa,xb,cnd)] print(r1) 対して、ndarrayに対して直に3項演算子を実行するwhereがある。 import numpy as np xa = np.array([1,2,3,4,5]) xb = np.array([6,7,8,9,10]) cnd = np.array([True,True,False,False,False]) r2 = np.where(cnd,xa,xb) print(r2) 数学関数,統計関数,次元削減 (n)次のndarrayをある軸について集計して(n-1)次のndarrayにする。 集計方法としていくつかの数学関数、統計関数が用意されている。 以下5x4(2次)のndarrayについて、それぞれの列について平均を取り4列(1次)のndarrayにしている。 さらに列の平均を取りスカラーにしている。 import numpy as np ary = np.random.randn(5,4) print(ary) # [[-1.84573174 1.84169514 1.43012623 -0.5416877 ] # [-1.03660701 0.63504086 -0.12239017 -0.77822113] # [ 0.1711323 -0.16660851 -0.7928288 1.17582814] # [-0.29302267 -0.23316282 1.70611457 0.53870384] # [-0.46513289 -1.12207588 0.01930695 0.49635739]] print(ary.mean(axis=0)) # [-0.6938724 0.19097776 0.44806576 0.17819611] print(ary.mean(axis=1)) # [ 0.22110048 -0.32554436 0.09688078 0.42965823 -0.26788611] print(ary.mean()) # 0.030841804893752683