sampling

あるとき、統計分析のため、
重みに応じたリサンプリングをする必要ができました。


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を返す関数です。

備忘録のため、ブログに載せてみましたが、
関数一つでできる方法がないのだろうか。

regression OLS

pythonで回帰分析(regression)を行うことができます

回帰分析は古くから知られており、汎用的で世間で一番用いられている統計分析法かもしれません。
回帰分析というのは、簡単に言うならば、説明変数から目的変数をダイレクトに計算する線形式を決めることと言えます。
もちろん、その間に因果関係がなさそうなデータでも、入力（説明）変数と出力（目的）変数のデータさえそろえば、計算できます。 なので、風の風速と桶（おけ）店の振り上げの関係式も分析できます。

推定性能については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が説明変数です。

回帰分析や検定については他をググってみてください。 もっと詳しく丁寧に書いてあるページがたくさんあります。

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でできました。

そのほかにもいくつか科学技術計算ライブラリがありますが、 長くなりそうなので、残りはまた今度ご紹介します。

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)

とすると

メソッド・メンバ変数の一覧が表示されます。

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形式の取り扱いについてはそのうちまた詳しく書きます。

ラン ダム文字列生成

ご無沙汰しております。
最近仕事で、任意の長さのランダム文字列を生成するスクリプトを作りました。

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.....'というようなすべての英数字を含んだ文字列が欲しかったので使いました。

もちろん、他の方法でもできますが、これが一番理解しやすくて簡潔だと思います。
是非使ってみてください。

mysql

今回はMySQLとの連携です。
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です。ここでは適当な値を入れていますが、それぞれ、設定した通り入力してください。

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を実行しないと、更新されません。