あるとき、統計分析のため、
重みに応じたリサンプリングをする必要ができました。
floatのベクトルvがあったとすると
p=v/sum(v)
を確率ベクトルとして、復元抽出するコードを考えてみました。
探しては見たのですが、すぐには見つかりませんでした。
そこで、したのような処理で実現してみました。
import numpy as np
n = 100
def resamp( v, n ):
vsum = np.sum(v)
vcum = np.cumsum(v)
indices = np.array([ np.argmax(1.0/(vcum - r)) for r in np.random.uniform(0,vsum,n) ])
return indices
確率ベクトルvと出力サンプル数nを入力し、インデックスのnumpy.arrayを返す関数です。
備忘録のため、ブログに載せてみましたが、
関数一つでできる方法がないのだろうか。
2015年11月25日水曜日
regression OLS
pythonで回帰分析(regression)を行うことができます
もし、上のコードでstatsmodelsが無いと言われたら、 easy_installで入れてください。 結果の表示は下のようになりました。
あとはAIC・BICや基本統計量が表示されています。 AIC・BICは値が低いほうがよりデータに近いことを示す情報量規準というものです。 それぞれの統計量についてはまた回を改めて説明します。 回帰分析関数クラスのolsへ入力する'y~x'というのはGNU Rと同様の書式です。 つまり、この場合、ｙが目的変数で右側のxが説明変数です。
回帰分析や検定については他をググってみてください。 もっと詳しく丁寧に書いてあるページがたくさんあります。
回帰分析は古くから知られており、汎用的で世間で一番用いられている統計分析法かもしれません。
回帰分析というのは、簡単に言うならば、説明変数から目的変数をダイレクトに計算する線形式を決めることと言えます。
もちろん、その間に因果関係がなさそうなデータでも、入力（説明）変数と出力（目的）変数のデータさえそろえば、計算できます。
なので、風の風速と桶（おけ）店の振り上げの関係式も分析できます。
推定性能についてはR二乗値などの値で評価します。 ここでは「通常の方法」についてpythonで行う方法についてご紹介します。 OLS法というのはOrdinary Least Square（普通最小二乗法）です。 回帰分析は最小二乗法で「行うべきもの」と考えている人もいるかもしれませんが、 他にもいくつか選択肢はあります。 例えば、重み付き最小二乗とか、ロバスト推定とか。 なので、一番使われている「通常の方法」がOrdinary Least Squareです。
import numpy
import statsmodels.formula.api as sfa
d = {'x':numpy.random.randn(10,3), 'y':numpy.random.randn(10) }
o = sfa.ols( 'y~x', data=d )
r = o.fit()
print r.summary()
もし、上のコードでstatsmodelsが無いと言われたら、 easy_installで入れてください。 結果の表示は下のようになりました。
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.403
Model: OLS Adj. R-squared: 0.104
Method: Least Squares F-statistic: 1.348
Date: Sat, 14 Nov 2015 Prob (F-statistic): 0.345
Time: 23:01:08 Log-Likelihood: -13.518
No. Observations: 10 AIC: 35.04
Df Residuals: 6 BIC: 36.25
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept 1.2127 0.745 1.628 0.155 -0.610 3.035
x[0] 0.6436 0.379 1.697 0.141 -0.284 1.572
x[1] -1.1225 0.837 -1.341 0.229 -3.171 0.926
x[2] -1.2621 1.742 -0.724 0.496 -5.526 3.001
==============================================================================
Omnibus: 0.620 Durbin-Watson: 1.691
Prob(Omnibus): 0.733 Jarque-Bera (JB): 0.417
Skew: -0.425 Prob(JB): 0.812
Kurtosis: 2.472 Cond. No. 8.32
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
出力がたくさんありすぎて、わからなくなってしまうかもしれません。
真ん中の少し下辺りに係数(x[*] coef)や切片(Intercept coef)の推定値が表示されています。
あとはAIC・BICや基本統計量が表示されています。 AIC・BICは値が低いほうがよりデータに近いことを示す情報量規準というものです。 それぞれの統計量についてはまた回を改めて説明します。 回帰分析関数クラスのolsへ入力する'y~x'というのはGNU Rと同様の書式です。 つまり、この場合、ｙが目的変数で右側のxが説明変数です。
回帰分析や検定については他をググってみてください。 もっと詳しく丁寧に書いてあるページがたくさんあります。
2015年11月14日土曜日
python 科学計算関係のライブラリWindows版
pythonは多くのライブラリをimportすることが多いです。
一つのことをするにしてもnumpyやscipy、またグラフを描くならmatplotlibとか。
PIPがうまく動作しない場合も多く、Windowsの場合インストーラを使わなければならないばあいが多いような気がします。 Linuxなら該当するパッケージ名さえわかれば、aptとpipでたいてい何とかなります。 Linxuの場合の話はまた今度します。
それぞれ、pipやeasy_installでインストールできるものもあれば、できないものもありややこしいので、現状でのインストール方法をまとめてみました(2015年10月現在）。
numpyやscipyはeasy_installやpipを使うとlapackなどのライブラリと接続できないなどというエラーがでて面倒なので、おとなしくサイトからsuperpackというものをダウンロードしてインストールしました。
ウィンドウズでバイナリをダウンロードしてくる場合はpythonのバージョン、32/64ビットの違いに注意してください。
pandasはeasy_installでうまくいかず、pipでインストールできました。 なぜかeasy_installではvcvarsall.batが無いと怒られます。そのときのエラーをコピペしておきます。
UPDATING build\lib.win32-2.7\pandas/_version.py set build\lib.win32-2.7\pandas/_version.py to '0.17.0' error: Setup script exited with error: Unable to find vcvarsall.bat
統計分析のstatsmodelsはeasy_installでインストールしました。
ただし、scipyやpandas、patsyを使っているようなので、それぞれをインストールしてからにしましょう。
patsyはpipでできました。
そのほかにもいくつか科学技術計算ライブラリがありますが、 長くなりそうなので、残りはまた今度ご紹介します。
2014年11月25日火曜日
DICOM
今回はかなり専門的な話です。
DICOMというものをご存知でしょうか。
たぶん、ほとんどの人は見る機会さえないので知らないと思いますが、医療用に用いられる画像形式（のようなもの）です。
例えば、レントゲンとかCT、MRIとかです。
昔はフィルムで撮影したものも技術の発展とともにデジタル化され、
今では撮影・撮像機械と病院内のコンピュータシステムが連携しています。
最近は都市部の病院は小さな開業医でも電子カルテを使っていたりして、急速に電子化が進んでいます。
DICOM形式はアメリカで開発されたもので、今では世界的な標準仕様になっています。
なぜ、普通のJPEGやPNGではだめなのかというと、医療用の画像には画像ピクセルそのものの情報だけではなく、撮影した状況、部位、患者情報、日付などが記されています。
仕事の関係でDICOMのヘッダ部分を分析する必要があったため
備忘録として簡単な導入部分をまとめてみました。
pythonはライブラリが豊富なので、DICOMライブラリも読み込むだけで使えます。
その前に、まずPIPでPydicomを インストールします。
$ pip install pydicom
例： 引数で与えるファイル名のDICOMファ イルを読み込み、患者ＩＤなどの情報を表示するスクリプト
ファイル名：showDicomInfo.py
#!/usr/bin/env python
import sys,dicom
d = dicom.read_file( sys.argv[1] )
print d
実行例：
$ ./showDicomInfo.py example.dcm
患者ＩＤのみを表示する場合、
#!/usr/bin/env python
import sys,dicom
d = dicom.read_file( sys.argv[1] )
print d.PatientID
ほかの値についてはdirメ ソッドで見てください。
つまり、上の例でいうと
print dir(d)
とすると
メソッド・メンバ変数の一覧が表示されます。
2014年11月10日月曜日
JSON文法チェック方法
データの送受信または記述のためよく使われるJSON(ジェイソン)という形式があります。
テキストなので、windowsなら標準で入っているメモ帳のようなテキストエディタで編集できます。
例えば、
{
"name" : ["shimotomai", "takayuki"],
"degree":"PhD"
}
こんな感じです。
入れ子構造にもでき、XML同様柔軟性が高い記述方式として最近よく使われています。
ただ、コンマ忘れやすいので、よく文法エラーになります。
もしくは最後の要素の後にコンマを入れたり、
コロンがなかったり、ダブルクォーテーションがシングルになっていたり。
ロボットの音声認識システムで、採用していたのでJSON形式よく利用していましたが、
記述が長いと、目で見てエラーをチェックするとか拷問に近いです。
あり得ません。
ということで、
自動的にチェックしてくれる方法をいくつか利用していました。
そのうちPythonを使って行う方法をご紹介します。
簡単です。
$ python -mjson.tool dummy.json
これだけです。
今回はプログラムさえ書きません。
エラーの場合
Expecting property name: line 12 column 2...
という風に行番号が表示されます。つまりこの場合は、12行目にエラーがありますよということです。
demjsonがインストールされていれば
jsonlint -v hoge.json
でもチェック可能です。
demjsonのインストールはubuntuの場合、
sudo aptitude install python-demjson
でできます。
pythonを使ったjson形式の取り扱いについてはそのうちまた詳しく書きます。
テキストなので、windowsなら標準で入っているメモ帳のようなテキストエディタで編集できます。
例えば、
{
"name" : ["shimotomai", "takayuki"],
"degree":"PhD"
}
こんな感じです。
入れ子構造にもでき、XML同様柔軟性が高い記述方式として最近よく使われています。
ただ、コンマ忘れやすいので、よく文法エラーになります。
もしくは最後の要素の後にコンマを入れたり、
コロンがなかったり、ダブルクォーテーションがシングルになっていたり。
ロボットの音声認識システムで、採用していたのでJSON形式よく利用していましたが、
記述が長いと、目で見てエラーをチェックするとか拷問に近いです。
あり得ません。
ということで、
自動的にチェックしてくれる方法をいくつか利用していました。
そのうちPythonを使って行う方法をご紹介します。
簡単です。
$ python -mjson.tool dummy.json
これだけです。
今回はプログラムさえ書きません。
エラーの場合
Expecting property name: line 12 column 2...
という風に行番号が表示されます。つまりこの場合は、12行目にエラーがありますよということです。
demjsonがインストールされていれば
jsonlint -v hoge.json
でもチェック可能です。
demjsonのインストールはubuntuの場合、
sudo aptitude install python-demjson
でできます。
pythonを使ったjson形式の取り扱いについてはそのうちまた詳しく書きます。
2014年11月8日土曜日
ラン ダム文字列生成
ご無沙汰しております。
最近仕事で、任意の長さのランダム文字列を生成するスクリプトを作りました。
CGIなどでセッション情報などをpythonで管理する場合、重複する可能性の低いランダムな文字列が必要になることがあります。
それで、英数字大小文字混合の２６＋２６＋１０＝６２種類をまぜて任意の長さの文字列を作ります。
下のソースが作成したスクリプトです。
import string
import numpy
m=100
s = ''.join(numpy.random.choice( list( string.digits + string.letters ), m ))
print s
pythonはたくさんのライブラリがあって便利な反面、簡単なことをするにしても、ライブラリを複数読み込む必要があります。
今回もnumpyとstringを読み込んでいます。
numpyは以前も説明した通り、python標準ではありませんので、別にインストールする必要があります。
ポイントとしては
listで文字列をリストに変換しているところと、
numpy.random.choiceで配列要素から100個サンプリングしているところです。
string.lettersは'abcd.....'というようなすべての英数字を含んだ文字列が欲しかったので使いました。
もちろん、他の方法でもできますが、これが一番理解しやすくて簡潔だと思います。
是非使ってみてください。
最近仕事で、任意の長さのランダム文字列を生成するスクリプトを作りました。
CGIなどでセッション情報などをpythonで管理する場合、重複する可能性の低いランダムな文字列が必要になることがあります。
それで、英数字大小文字混合の２６＋２６＋１０＝６２種類をまぜて任意の長さの文字列を作ります。
下のソースが作成したスクリプトです。
import string
import numpy
m=100
s = ''.join(numpy.random.choice( list( string.digits + string.letters ), m ))
print s
pythonはたくさんのライブラリがあって便利な反面、簡単なことをするにしても、ライブラリを複数読み込む必要があります。
今回もnumpyとstringを読み込んでいます。
numpyは以前も説明した通り、python標準ではありませんので、別にインストールする必要があります。
ポイントとしては
listで文字列をリストに変換しているところと、
numpy.random.choiceで配列要素から100個サンプリングしているところです。
string.lettersは'abcd.....'というようなすべての英数字を含んだ文字列が欲しかったので使いました。
もちろん、他の方法でもできますが、これが一番理解しやすくて簡潔だと思います。
是非使ってみてください。
2014年7月22日火曜日
mysql
今回はMySQLとの連携です。
PythonからMySQLに接続するためのライブラリはいくつかあります。
ここではMySQL標準の
mysql-connector-python
を使います。
MySQLのConnectorのページに行けばインストーラがあるはずです。
http://dev.mysql.com/downloads/connector/
ここからインストーラをダウンロードしてインストールします。
呼び出すコードは
import mysql.connector
補足ですが、
MySQLに保存したデータの実体は
C:\ProgramData\MySQL\MySQL Server 5.6\data
というところにあります。
この下にSchema名のフォルダがあります。そのフォルダ下にtable名のファイルがあります。
commitすると更新されるようです。
逆に言うと、commitを実行しないと、更新されません。
PythonからMySQLに接続するためのライブラリはいくつかあります。
ここではMySQL標準の
mysql-connector-python
を使います。
MySQLのConnectorのページに行けばインストーラがあるはずです。
http://dev.mysql.com/downloads/connector/
ここからインストーラをダウンロードしてインストールします。
呼び出すコードは
import mysql.connector
con = mysql.connector.connect( user='x', password='xyz', host='127.0.0.1', database='test' )
cur = con.cursor()
sql = 'CREATE TABLE hoge( id int, value double, comment varchar(32) )'
cur.execute( sql )
con.commit()
cur.close()
con.close()
という感じです。
connectの時の引数は適宜入れ替えて使って下さい。
特にuserとpasswordです。ここでは適当な値を入れていますが、それぞれ、設定した通り入力してください。
特にuserとpasswordです。ここでは適当な値を入れていますが、それぞれ、設定した通り入力してください。
sqlの中身に関しても同様です。今回はテーブルを作成していますが、SQLクエリを作成してexecuteに入れます。
ここはsqliteとだいたい一緒です。
データを引き出してくるときは
cur.execute( 'SELECT * from hoge;' )
for r in cur:
print r
という感じで取ってこれます。
これもsqliteの時とだいたい一緒です。
sqliteとは違うのは、DBのサーバ設定が必要なところです。
そっちはMySQLのサイトなどを参照してください。
補足ですが、
MySQLに保存したデータの実体は
C:\ProgramData\MySQL\MySQL Server 5.6\data
というところにあります。
この下にSchema名のフォルダがあります。そのフォルダ下にtable名のファイルがあります。
commitすると更新されるようです。
逆に言うと、commitを実行しないと、更新されません。
2014年6月23日月曜日
OpenCV 4
Contourとは等高線を指しますが、輪郭という意味もあります。
OpenCVではcontourを輪郭情報という意味で使用しています。
cv2の中にはfindContoursというcontourを計算する関数があります。
つまり、画像から輪郭を抽出することができます。
たとえば、
ret, thresh = cv2.threshold( imgray, 127, 255, 0 )
c = cv2.findContours( thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE )
とすることによりimgrayの画像を閾値処理して輪郭抽出します。
ただし、imgrayというのは２値化された画像で、通常のカラー画像を入れるとエラーになります。
カラー画像imを２値画像gimに変換する場合は、
gim = cv2.cvtColor( im, cv2.cv.CV_BGR2GRAY )
とします。
contourを描画する場合、
cv2.drawContours( im, c, -1, (0,255,0), 3 )
とすると緑色で輪郭情報cを画像imに描画します。
輪郭情報を最小の長方形領域に納める場合
x,y,w,h = cv2.boundingRect(c)
とします。
この長方形領域は
cv2.rectangle(im,(x,y),(x+w,y+h),(0,255,0),2)
とすることによって描画できます。
以上を踏まえて、
WEBカメラから取得した画像から輪郭を抽出して描画するプログラムを作ってみました。
ENTERキーを押すと終了します。
import cv2
import numpy as nmp
cv2.namedWindow('CAMERA')
v = cv2.VideoCapture(0)
v.grab()
im = v.read()
### to binary
gim = cv2.cvtColor( im[1], cv2.cv.CV_BGR2GRAY )
### thresholded
ret, thresh = cv2.threshold( gim, 127, 255, 0 )
c = cv2.findContours( thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE )
for i, cc in enumerate(c[0]):
cv2.drawContours( im[1], cc, -1, (50*i/255,255,i%255), 3 )
print '%d,255,%d' % ( 20*i/255, i%255)
cv2.imshow( 'CAMERA', im[1] )
cv2.waitKey(0)
ラボの画像から輪郭抽出してみました。
たくさんの小さい輪郭ができてしまいました。
実際に使う場合は、いかにして目的とするものの輪郭を計算するかという工夫が重要だと思います。
OpenCVではcontourを輪郭情報という意味で使用しています。
cv2の中にはfindContoursというcontourを計算する関数があります。
つまり、画像から輪郭を抽出することができます。
たとえば、
ret, thresh = cv2.threshold( imgray, 127, 255, 0 )
c = cv2.findContours( thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE )
とすることによりimgrayの画像を閾値処理して輪郭抽出します。
ただし、imgrayというのは２値化された画像で、通常のカラー画像を入れるとエラーになります。
カラー画像imを２値画像gimに変換する場合は、
gim = cv2.cvtColor( im, cv2.cv.CV_BGR2GRAY )
とします。
contourを描画する場合、
cv2.drawContours( im, c, -1, (0,255,0), 3 )
とすると緑色で輪郭情報cを画像imに描画します。
輪郭情報を最小の長方形領域に納める場合
x,y,w,h = cv2.boundingRect(c)
とします。
この長方形領域は
cv2.rectangle(im,(x,y),(x+w,y+h),(0,255,0),2)
とすることによって描画できます。
以上を踏まえて、
WEBカメラから取得した画像から輪郭を抽出して描画するプログラムを作ってみました。
ENTERキーを押すと終了します。
import cv2
import numpy as nmp
cv2.namedWindow('CAMERA')
v = cv2.VideoCapture(0)
v.grab()
im = v.read()
### to binary
gim = cv2.cvtColor( im[1], cv2.cv.CV_BGR2GRAY )
### thresholded
ret, thresh = cv2.threshold( gim, 127, 255, 0 )
c = cv2.findContours( thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE )
for i, cc in enumerate(c[0]):
cv2.drawContours( im[1], cc, -1, (50*i/255,255,i%255), 3 )
print '%d,255,%d' % ( 20*i/255, i%255)
cv2.imshow( 'CAMERA', im[1] )
cv2.waitKey(0)
ラボの画像から輪郭抽出してみました。
たくさんの小さい輪郭ができてしまいました。
実際に使う場合は、いかにして目的とするものの輪郭を計算するかという工夫が重要だと思います。
2014年6月9日月曜日
フラクタル 第4回
フラクタルシリーズの4回めです。
今回も、マンデルブロのアレンジを作ってみました。
ただ、毎回関数を別に作るのは無駄なので、P乗で一般化してみました。
つまり、
z(n+1) = z(n)**p + c
z(0)=0
とした時のcと最終のｚの値を出力します。
import numpy as np
def makeMandelbrotP( itrn, p=2, n=100, xrange=None, yrange=None ):
if xrange == None:
xrange = (-2.5, 2.5)
if yrange == None:
yrange = (-2.5, 2.5)
r = np.random.uniform( xrange[0], xrange[1], n) + 1j*np.random.uniform( yrange[0], yrange[1],n)
c = r.copy()
z = 0
for i in range(itrn):
z = z**p + c
return c, z
今回はp=4として4乗にした時、このｚの値が発散しているかどうかを複素平面上に表示しました。
予想通り、120度の回転対称です。
前回の3乗で180度対称だったので、p乗で306/(p-1) 度の回転対称になるものと想像できます。
この図は16回の繰り返しです。
繰り返し回数については、8回ぐらいでだいたいフラクタルっぽくなります。
一枚あたりのデータ量が多く、リスト内表記を使うとメモリアロケーションのエラーになるので、
仕方なくforループを使いました。
通常書類作成用に使っているPC(Core2Quad 2.4GHz)で、上の画像作成時間は一枚あたり6分位かかりました。
今回も、マンデルブロのアレンジを作ってみました。
ただ、毎回関数を別に作るのは無駄なので、P乗で一般化してみました。
つまり、
z(n+1) = z(n)**p + c
z(0)=0
とした時のcと最終のｚの値を出力します。
import numpy as np
def makeMandelbrotP( itrn, p=2, n=100, xrange=None, yrange=None ):
if xrange == None:
xrange = (-2.5, 2.5)
if yrange == None:
yrange = (-2.5, 2.5)
r = np.random.uniform( xrange[0], xrange[1], n) + 1j*np.random.uniform( yrange[0], yrange[1],n)
c = r.copy()
z = 0
for i in range(itrn):
z = z**p + c
return c, z
今回はp=4として4乗にした時、このｚの値が発散しているかどうかを複素平面上に表示しました。
前回の3乗で180度対称だったので、p乗で306/(p-1) 度の回転対称になるものと想像できます。
この図は16回の繰り返しです。
繰り返し回数については、8回ぐらいでだいたいフラクタルっぽくなります。
繰り返し2回
繰り返し3回
繰り返し4回
繰り返し5回
繰り返し6回
繰り返し7回
繰り返し8回
という感じです。一枚あたりのデータ量が多く、リスト内表記を使うとメモリアロケーションのエラーになるので、
仕方なくforループを使いました。
通常書類作成用に使っているPC(Core2Quad 2.4GHz)で、上の画像作成時間は一枚あたり6分位かかりました。
2014年6月6日金曜日
pythonからsqlite(DB)の使用
今日はsqliteの話です。
簡単なことしか書いていませんが、しばらく使わないと手順を忘れてしまうこともありますので自分自身の備忘録です。
シミュレーションや分析で、パラメータが多いとデータの管理に手間がかかります。
一度のシミュレーションデータが膨大な場合もどのように保存するかで後々の分析の効率が悪くなったり問題が起こったりします。
一度のシミュレーションにかかる時間が長い（1日とか）と次にどのパラメータでシミュレーション（または分析）をするべきか悩むことがあります。
また、観点を変えて分析をしてみたいと思ったときに、データがそろっているとは限りません。判断をするために足りないデータが何なのかを知る必要があります。
そこで、データの管理のため、データベースを使うと便利なことがあります。
python2.6では標準モジュールとしてsqliteが入っています。
sqliteはフリーのデータベースです。
使い方は
データベースに接続
c = sqlite.connect( 'hogehoge.db' )
切断
c.close()
SQLクエリ送信
sql = 'create table ccc ( iteration integer, realpart float, imagpart float );'
c.execute( sql )
結果受信
sql = 'select * from ccc'
cc = c.cursor()
for r in cc.execute( sql ):
print r
という感じで使います。
大事なのが、コミット
c.commit()
です。
とりあえずは「これで保存します」の意味だと思って下さって構わないと思います。
これがないと保存されません。
SQLに関してはネット上を検索すればわかりやすいサイトがたくさんあると思いますので、そちらに任せます。
簡単なことしか書いていませんが、しばらく使わないと手順を忘れてしまうこともありますので自分自身の備忘録です。
シミュレーションや分析で、パラメータが多いとデータの管理に手間がかかります。
一度のシミュレーションデータが膨大な場合もどのように保存するかで後々の分析の効率が悪くなったり問題が起こったりします。
一度のシミュレーションにかかる時間が長い（1日とか）と次にどのパラメータでシミュレーション（または分析）をするべきか悩むことがあります。
また、観点を変えて分析をしてみたいと思ったときに、データがそろっているとは限りません。判断をするために足りないデータが何なのかを知る必要があります。
そこで、データの管理のため、データベースを使うと便利なことがあります。
python2.6では標準モジュールとしてsqliteが入っています。
sqliteはフリーのデータベースです。
使い方は
データベースに接続
c = sqlite.connect( 'hogehoge.db' )
切断
c.close()
SQLクエリ送信
sql = 'create table ccc ( iteration integer, realpart float, imagpart float );'
c.execute( sql )
結果受信
sql = 'select * from ccc'
cc = c.cursor()
for r in cc.execute( sql ):
print r
という感じで使います。
大事なのが、コミット
c.commit()
です。
とりあえずは「これで保存します」の意味だと思って下さって構わないと思います。
これがないと保存されません。
SQLに関してはネット上を検索すればわかりやすいサイトがたくさんあると思いますので、そちらに任せます。
2014年5月29日木曜日
フラクタル 第3回
今度は前回のマンデルブロ集合の変則的な形として、試してみたものです。
今回は
z(n+1) = z(n)^3 + c
という漸化式による収束性について複素平面の図を作ってみました。
できたのは下のような左右対称な図です。
各ｃ（複素数）の点で、ｚの発散しなかった領域が黒い部分です。緑部分は発散領域です。
この場合もフラクタルのような図形になりました。
図の中央が原点で、左右対称な図形となりました。
範囲は実数軸、虚数軸ともに-1.5から1.5までです。
データ点は約14M 個生成して作成しました。
漸化式の繰り返し回数は10回です。
今回はpythonからデータを生成して、MySQLを使ってデータを管理しましたが、30分程度かかりました。
今回の計算結果はWEBを検索しても簡単に見つからなかったので、ちょっとした課題としていいかもしれません。
作った関数はこんなかんじです。
def makeMandelbrot3( itrn, n=100, xrange=None, yrange=None ):
if xrange == None:
xrange = (-2.5, 2.5)
if yrange == None:
yrange = (-2.5, 2.5)
c = np.random.uniform( xrange[0], xrange[1], n) + 1j*np.random.uniform( yrange[0], yrange[1],n)
z = 0
for i in range(itrn):
z = z**3 + c
return c, z
itrnは繰り返し回数で、データ数はnです。
戻り値はcとzなので、
この関数を使用するときは
c, z = makeMandelbrot3( 10, 1000 )
のようにします。
これで繰り返し10回のデータcとzが1000点分生成されます。
メモリ確保の問題さえなければ、nの値を大きくすればその分データが生成されます。
2014年5月21日水曜日
bool
ブール型の演算について
どの言語もそうですが、整数型や浮動小数点型の演算はだいたい同じですが、ブール型の演算は少しずつやり方が異なっているような気がします。
pythonでの扱い方について記録しておきます。
a = True
b = False
a and b
とすると
False
となります。
a+b
とすると
1
となり、
a+a
とすると
2
になります。
つまり、ブール型はTrueが1、Falseが0の値になっていると考えて構わないと思います。
それで、今度はnumpy.arrayでの扱いです。
aa = numpy.array( [ a, a, b] )
とすると
array([ True, True, False], dtype=bool)
となり、bool型として扱われていることがわかります。
で、今度は
aa and aa
とするとエラーになります。
実行した結果のエラーメッセージは下のようになっていました。
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: The truth value of an array with more than one element is ambiguous.
Use a.any() or a.all()
numpy.arrayは他の型ではそれぞれの要素に関して個別に演算するので、ブール型もそうかと思うとそうでもないようです。
で、今度は
aa = numpy.array( [True, True, False], dtype=bool)
bb = numpy.array( [False, True, False], dtype=bool)
aa + bb
aa * bb
とすると、それぞれ下のようになりました。
>>> aa+bb
array([ True, True, False], dtype=bool)
>>> aa*bb
array([False, True, False], dtype=bool)
つまり、要素個別に演算する場合は+と*でandとorの代用する必要があるということです。
andとor,notは単独のTrue, Falseといったブール値にしか使えません。
notで実現していた、反転をさせる場合はどうするかというと、
~aa
とします。
結果は
array([False, False, True], dtype=bool)
となります。
ブールのarrayは結構有用です。
たとえば、
x = numpy.random.randint( 0, 10, 12 )
c = x>5
print x[c]
というように、条件にあったものだけ抽出する際に使います。
cは
array([ True, False, True, False, True, False, True, True, True,
False, False, False], dtype=bool)
といった風にbool型のarrayになっています。
どの言語もそうですが、整数型や浮動小数点型の演算はだいたい同じですが、ブール型の演算は少しずつやり方が異なっているような気がします。
pythonでの扱い方について記録しておきます。
a = True
b = False
a and b
とすると
False
となります。
a+b
とすると
1
となり、
a+a
とすると
2
になります。
つまり、ブール型はTrueが1、Falseが0の値になっていると考えて構わないと思います。
それで、今度はnumpy.arrayでの扱いです。
aa = numpy.array( [ a, a, b] )
とすると
array([ True, True, False], dtype=bool)
となり、bool型として扱われていることがわかります。
で、今度は
aa and aa
とするとエラーになります。
実行した結果のエラーメッセージは下のようになっていました。
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: The truth value of an array with more than one element is ambiguous.
Use a.any() or a.all()
numpy.arrayは他の型ではそれぞれの要素に関して個別に演算するので、ブール型もそうかと思うとそうでもないようです。
で、今度は
aa = numpy.array( [True, True, False], dtype=bool)
bb = numpy.array( [False, True, False], dtype=bool)
aa + bb
aa * bb
とすると、それぞれ下のようになりました。
>>> aa+bb
array([ True, True, False], dtype=bool)
>>> aa*bb
array([False, True, False], dtype=bool)
つまり、要素個別に演算する場合は+と*でandとorの代用する必要があるということです。
andとor,notは単独のTrue, Falseといったブール値にしか使えません。
notで実現していた、反転をさせる場合はどうするかというと、
~aa
とします。
結果は
array([False, False, True], dtype=bool)
となります。
ブールのarrayは結構有用です。
たとえば、
x = numpy.random.randint( 0, 10, 12 )
c = x>5
print x[c]
というように、条件にあったものだけ抽出する際に使います。
cは
array([ True, False, True, False, True, False, True, True, True,
False, False, False], dtype=bool)
といった風にbool型のarrayになっています。
2014年5月14日水曜日
2014年5月12日月曜日
フラクタル
pythonで複素数の練習を兼ねて
フラクタル図形を作ってみました。
マンデルブロ集合というものです。
フラクタルが何かについては専門書または他のサイトなどで調べてみてください。
マンデルブロ集合の定義は簡単な複素数漸化式で記述されます。
つまりマンデルブロ図形というのはある複素数cについて、
z(n+1) = z(n)*z(n) + c
z(0) = 0
とした場合、n→∞で発散するか否かを色分けして複素平面に表示したものです。
複素数平面の点をランダムに選び、n回漸化式を解いてz(n)を作ります。
その時の、
log |z(n)|
の値を画素の輝度値としました。
forループを2重に使うと時間がかかるので、numpy.arrayで一気に計算します。
m = 10000
c = numpy.random.uniform( -2.5,2.5,m ) + 1j * numpy.random.uniform( -2.5,2.5,m )
z = numpy.zeros( m, dtype=numpy.complex )
for i in range(16):
z = z**2 + c
と言う感じです。
opencvを使って画像化します。
しかし、案外時間が掛かるので、zの漸化式計算と、画像化を分けて処理しました。
30Mポイントのサンプリングをして、そこから画像化しましたが、埋め切れていません。画像を拡大してみるとポツポツと黒い点が見えます。
実は細くつながっています。が、データ点数（サンプリング）で埋め切れていないのと画像化の解像度が荒いため、つながっているようには見えません。
フラクタル図形を作ってみました。
マンデルブロ集合というものです。
フラクタルが何かについては専門書または他のサイトなどで調べてみてください。
マンデルブロ集合の定義は簡単な複素数漸化式で記述されます。
つまりマンデルブロ図形というのはある複素数cについて、
z(n+1) = z(n)*z(n) + c
z(0) = 0
とした場合、n→∞で発散するか否かを色分けして複素平面に表示したものです。
複素数平面の点をランダムに選び、n回漸化式を解いてz(n)を作ります。
その時の、
log |z(n)|
の値を画素の輝度値としました。
forループを2重に使うと時間がかかるので、numpy.arrayで一気に計算します。
m = 10000
c = numpy.random.uniform( -2.5,2.5,m ) + 1j * numpy.random.uniform( -2.5,2.5,m )
z = numpy.zeros( m, dtype=numpy.complex )
for i in range(16):
z = z**2 + c
と言う感じです。
opencvを使って画像化します。
しかし、案外時間が掛かるので、zの漸化式計算と、画像化を分けて処理しました。
Mandelbrot set
30Mポイントのサンプリングをして、そこから画像化しましたが、埋め切れていません。画像を拡大してみるとポツポツと黒い点が見えます。
実は細くつながっています。が、データ点数（サンプリング）で埋め切れていないのと画像化の解像度が荒いため、つながっているようには見えません。
2014年5月9日金曜日
複素数
python numpyでの複素数の扱い方についてです。
複素数は型としてはnumpy.complexと定義されています。
これを利用してみます。
つまり、
numpy.array( [3, 2, 1], dtype=numpy.complex )
とすると
array([ 3.+0.j, 2.+0.j, 1.+0.j])
というふうに複素数配列となります。
複素数の演算もできます
a = numpy.array( [3, 2, 1], dtype=numpy.complex )
複素数は型としてはnumpy.complexと定義されています。
これを利用してみます。
つまり、
numpy.array( [3, 2, 1], dtype=numpy.complex )
とすると
array([ 3.+0.j, 2.+0.j, 1.+0.j])
というふうに複素数配列となります。
複素数の演算もできます
a = numpy.array( [3, 2, 1], dtype=numpy.complex )
a ** 2
とすると2乗します。
実行結果は
array([ 9.+0.j, 4.+0.j, 1.+0.j])
となります。
指数計算はそのまま掛けて計算されます。
つまり、python風な数式で表すと
z=x+jy
とすると
z**2 = (x+jy) * (x+jy) = x**2 - y**2 + j 2*x*y
となります。
複素数には色々と実数にはなかった計算方式が出てくるので注意が必要です。
実数部は
numpy.real(z)
虚数部は
numpy.imag(z)
とします。
これらの関数はzがarrayでもOKです。その場合、出力もarrayになります。
複素共役は
numpy.conj(z)
とします。
複素共役が何かわからない場合は検索してみて下さい。
誰かがわかりやすく説明してくれていると思います。
絶対値はnumpy.abs(z)
で計算できます。
複素数の乱数発生は簡単にはできませんでした。
一様分布の場合
z = numpy.random.uniform( -5, 5, 10 ) + 1j * numpy.random.uniform( -5, 5, 10 )
としました。
結果は
array([-4.48658518-4.19169279j, 2.39542899+4.87030019j,
4.03412469-1.00727465j, 1.72313031+0.99345025j,
0.68053779-0.385337j , 2.08333535-0.32240792j,
-0.65638447-0.27364433j, 1.36415791-1.89435328j,
0.19786534-4.51657661j, 2.44853396+3.04850488j])
となりました。
もっとスマートな方法があれば、是非ご教示下さい。
2014年5月7日水曜日
svgwriteのエラー
pythonでSVGファイルを生成する必要があったので、
svgwriteというモジュールをインストールしてみました。
インストールは簡単にできます。
pip install svgwrite
です。
が
サンプルを実行してみようとすると
エラー…
import下だけなのにエラーが出ました。
エラーメッセージを貼り付けると下のようにxml関連らしいです。
>>> import svgwrite
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "c:\Python26\Lib\site-packages\svgwrite\__init__.py", line 51, in <module
>
from svgwrite.drawing import Drawing
File "c:\Python26\Lib\site-packages\svgwrite\drawing.py", line 29, in <module>
from svgwrite.container import SVG, Defs
File "c:\Python26\Lib\site-packages\svgwrite\container.py", line 27, in <modul
e>
from svgwrite.base import BaseElement
File "c:\Python26\Lib\site-packages\svgwrite\base.py", line 12, in <module>
from svgwrite.etree import etree
File "c:\Python26\Lib\site-packages\svgwrite\etree.py", line 33, in <module>
original_serialize_xml = etree._serialize_xml
AttributeError: 'module' object has no attribute '_serialize_xml'
下のサイトを参考に少し修正しました。
https://bitbucket.org/mozman/svgwrite/pull-request/3/fix-etreepy/diff
pythonインストールフォルダに
Lib/site-packages
というフォルダがあります。
この下にいわゆるライブラリがインストールされます。
インストールされていれば、
この下にsvgwriteというフォルダが有ると思います。
その中にetree.py
というファイルがあるので、
33行目
original_serialize_xml = etree._serialize_xml
を削除して、最終行に移動します。
現状のエラーでは循環呼び出しをしてしまったらしいです。
これで、ひとまずimportは出来ました。
使い方などはまた次の機会にレポートします。
svgwriteというモジュールをインストールしてみました。
インストールは簡単にできます。
pip install svgwrite
です。
が
サンプルを実行してみようとすると
エラー…
import下だけなのにエラーが出ました。
エラーメッセージを貼り付けると下のようにxml関連らしいです。
>>> import svgwrite
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "c:\Python26\Lib\site-packages\svgwrite\__init__.py", line 51, in <module
>
from svgwrite.drawing import Drawing
File "c:\Python26\Lib\site-packages\svgwrite\drawing.py", line 29, in <module>
from svgwrite.container import SVG, Defs
File "c:\Python26\Lib\site-packages\svgwrite\container.py", line 27, in <modul
e>
from svgwrite.base import BaseElement
File "c:\Python26\Lib\site-packages\svgwrite\base.py", line 12, in <module>
from svgwrite.etree import etree
File "c:\Python26\Lib\site-packages\svgwrite\etree.py", line 33, in <module>
original_serialize_xml = etree._serialize_xml
AttributeError: 'module' object has no attribute '_serialize_xml'
下のサイトを参考に少し修正しました。
https://bitbucket.org/mozman/svgwrite/pull-request/3/fix-etreepy/diff
pythonインストールフォルダに
Lib/site-packages
というフォルダがあります。
この下にいわゆるライブラリがインストールされます。
インストールされていれば、
この下にsvgwriteというフォルダが有ると思います。
その中にetree.py
というファイルがあるので、
33行目
original_serialize_xml = etree._serialize_xml
を削除して、最終行に移動します。
現状のエラーでは循環呼び出しをしてしまったらしいです。
これで、ひとまずimportは出来ました。
使い方などはまた次の機会にレポートします。
2014年4月25日金曜日
arrayデータの16進パターン
numpyのarrayにはいくつかのメソッドがあります。
覚えておくとたまに便利になる機能がいくつかあります。
まず、tostringです。
動作確認しつつまとめてみます。
tostringメソッドは文字通り、arrayのデータをstringデータとして出力します。
数値を文字列に変換するわけではありません。
どういうときに使用するかというと、場合それぞれですが、
たとえば、シリアル通信やTCP/IP通信などで文字列として受け取ったデータを、数値として扱いたい場合などに使えます。
最近では、画像データの変換として使いました。
例を挙げます。
そこで、x
つまり、arangeで作成した段階で、4[byte](32bit)のint型として扱われ、リトルエンディアン方式で保管されているようです。
これを文字列として保存して、arrayとして戻すためには、
浮動小数点の場合はdtypeをfloatに指定します。
その他、32bitや64bit、または整数型か浮動小数点型か、符号有無などの細かい設定が必要な場合は、numpyの下に定義されている型を利用します。
たとえば
符号無整数
numpy.uint8
numpy.uint16
numpy.uint32
numpy.uint64
符号付整数
numpy.int8
numpy.int16
numpy.int32
numpy.int64
浮動小数点
numpy.float16
numpy.float32
...
複素数
numpy.complex
numpy.complex64
...
などです。
ちなみに、数値データをそれぞれ文字列に変換する場合は、いくつか方法がありますが、
listとして出力するなら
今度、それぞれの方法について処理時間などの比較をしてみます。
例で使用したnumpyは1.8.1 (python 2.6)を使用しています。
覚えておくとたまに便利になる機能がいくつかあります。
まず、tostringです。
動作確認しつつまとめてみます。
数値を文字列に変換するわけではありません。
どういうときに使用するかというと、場合それぞれですが、
たとえば、シリアル通信やTCP/IP通信などで文字列として受け取ったデータを、数値として扱いたい場合などに使えます。
最近では、画像データの変換として使いました。
例を挙げます。
import numpyとするとxの中身は
x = numpy.arange(-3,3)
array([-3, -2, -1, 0, 1, 2])となります。
そこで、x
s = x.tostring()とすると、sの中身は
'\xfd\xff\xff\xff\xfe\xff\xff\xff\xff\xff\xff\xff\x00\x00\x00\x00\x01\x00\x00\x00\x02\x00\x00\x00'となり、16進のバイナリパターンが文字列として出力されています。
つまり、arangeで作成した段階で、4[byte](32bit)のint型として扱われ、リトルエンディアン方式で保管されているようです。
これを文字列として保存して、arrayとして戻すためには、
numpy.frombuffer( s, dtype=int )とします。結果は、
array([-3, -2, -1, 0, 1, 2])となります。
浮動小数点の場合はdtypeをfloatに指定します。
その他、32bitや64bit、または整数型か浮動小数点型か、符号有無などの細かい設定が必要な場合は、numpyの下に定義されている型を利用します。
たとえば
符号無整数
numpy.uint8
numpy.uint16
numpy.uint32
numpy.uint64
符号付整数
numpy.int8
numpy.int16
numpy.int32
numpy.int64
浮動小数点
numpy.float16
numpy.float32
...
複素数
numpy.complex
numpy.complex64
...
などです。
ちなみに、数値データをそれぞれ文字列に変換する場合は、いくつか方法がありますが、
a.astype( str )とすればできます。
listとして出力するなら
map( str, a )としてもよいかもしれません。
今度、それぞれの方法について処理時間などの比較をしてみます。
例で使用したnumpyは1.8.1 (python 2.6)を使用しています。
2014年4月24日木曜日
数理統計学・データサイエンスとpython
pythonを使って統計アルゴリズムを記述することができます。
現在、統計アルゴリズムをコーディングする上でいくつかの比較的簡便な方法が開発されてきました。統計言語Sを基にしたR言語もその一つです。
ここでは最近急速に発展してきたpythonを使った方法について考察してみます。
最近の流れとしては大規模データ、いわゆるビッグデータの分析に注目が集まっているようです。
pythonはこれまでにも書いてきたようにスクリプト言語で基本インタプリタですので、書き方によっては計算時間がネックになってしまうことがあります。
さらには、もともと統計分析を目的として開発されたSまたはR言語のような機能の豊富さには負けます。
しかし、pythonはそれを上回るほどの開発効率やソースコードの簡潔さが特徴で、ライブラリも多岐にわたるため応用範囲が広いという利点があります。
科学技術計算で使うためのライブラリとしては、有名なものだけでも
- Numpy
- Scipy
- Sklearn
- Pandas
- Jubatus
などがあります。
Numpyについてはこれまでも何回か説明してきましたが、統計に限らず数値演算の基本的な機能を提供します。
Numpyのもっとも重要な機能はArrayだと思います。
Arrayというのはつまり配列ですが、python標準のlistより、ベクトル的な計算を意識したもので、これを効果的に使うことにより、効率的な計算が可能です。
NumpyにはMatrixという型もあります。
Matrixを使うと行列演算が簡単にかけるので便利ではありますが、
Arrayの方が使い勝手が良いと思います。
ScipyはNumpyではカバーしきれなかった科学技術寄りの機能を実現しています。
たとえば、統計アルゴリズムの例でいうと、乱数発生はNumpyの機能が使えます。つまり、
numpy.random.randn(100)
などとすると100個の正規乱数がnumpy.array形式で生成されます。
しかし、正規分布の確率密度関数のグラフが書きたいとなれば、
scipy.statsの機能が使えます。
import numpy
import scipy.stats
p = scipy.stats.norm.pdf( numpy.arange(-3,3,0.1) )
とすれば、pに-3から3までの正規密度関数の値が入ります。
正規分布のほかにもcauchy分布やgammaやvon misesなどの分布もあります。
統計の論文に出てきそうな標準的確率分布はあらかたあります。
残りのライブラリの説明はまた今度します。
興味があれば、それぞれググってみてください。
2014年3月17日月曜日
C言語との連携 ctypes
pythonは高機能で大変便利ですが、処理が遅いので困る場合があります。
そういう場合は、いろいろと工夫して高速化する方法があります。
ctypesによるDLLの利用はその高速化手法の一つです。
つまり、C言語（など）でDLLを作成して、pythonから読み込みます。
そのDLLを呼びだすためのライブラリがctypesです。
import ctypes
d = ctypes.cdll.LoadLibrary("libx.dll")
C言語と言いましたが、実際C++で作成しても構いません。
但し、呼びだす関数部分はC言語の方式で書かなければなりません。
つまり裏側ではC++のclass定義やインスタンスを生成しても構いませんが、
pythonから呼び出す部分は
extern "C"
などとしてＣ言語方式で、関数にしなければならないということです。
2014年3月13日木曜日
グラフ
グラフはプレゼンするうえでは必須で、見やすいグラフをどうやって作るかは重要な課題です。
Excelは最近のバージョンで少し見やすくなりましたが、昔はかなりひどいものでした。
これまでMATLABを使ったりRを使ったりしていましたが、
現在は
pythonのグラフ描画ライブラリであるmatplotlibを使っています。
結構グラフがきれいに描画できるのでお勧めです。
import matplotlib.pyplot as ppl
import numpy
t = numpy.arange(100)
ppl.plot( t, numpy.sin( t ) )
ppl.show()
こんな感じで簡単にグラフが描画できます。
行列の値に応じて色分けして表示する場合は、
a = numpy.random.randn(10,15)
ppl.pcolor( a )
とすると表示できます。
Excelは最近のバージョンで少し見やすくなりましたが、昔はかなりひどいものでした。
これまでMATLABを使ったりRを使ったりしていましたが、
現在は
pythonのグラフ描画ライブラリであるmatplotlibを使っています。
結構グラフがきれいに描画できるのでお勧めです。
import matplotlib.pyplot as ppl
import numpy
t = numpy.arange(100)
ppl.plot( t, numpy.sin( t ) )
ppl.show()
こんな感じで簡単にグラフが描画できます。
行列の値に応じて色分けして表示する場合は、
a = numpy.random.randn(10,15)
ppl.pcolor( a )
とすると表示できます。
登録:
投稿 (Atom)