![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
* Pythonによるブラウン運動データ解析 [#c9a28c9e] このページではPythonを用いてブラウン運動実験で得られたデータを解析する方法について解説をする。データ解析も実験課題の一部なので、わざと不完全な解説になっている。不明な点は、各種参考文献やWeb siteを参照して自習すること。 #contents ** インストール [#p520d870] Python自体は、[[PythonのWebsite>http://www.python.org/]]や[[日本Pythonユーザー会>http://www.python.jp/Zope]]などからダウンロードできる。~ 一般に、Pythonは、[[PythonのWebsite>http://www.python.org/]]や[[日本Pythonユーザー会>http://www.python.jp/Zope]]などからダウンロードできる。 以下に、各種OSにおけるインストール方法を概説する。 Linuxの場合は、Distributionのパッケージシステムからインストールするのが楽で確実である。最初からインストール済みのDistributionも多い。 また、下で説明する[[Enthought Python Distribution (EPD)>http://www.enthought.com/products/epd.php]]を利用する方法もある。 この場合、科学技術計算環境が一気に揃うので便利である。 Windowsの場合、[[Python(x,y)>http://www.pythonxy.com/]]というプロジェクトから、Pythonと科学技術計算に必要な拡張モジュール類が一気にダウンロード、インストールできる。Python(x,y)はLinux版もある。 Windowsの場合、[[Python(x,y)>http://www.pythonxy.com/]]というプロジェクトから、Pythonと科学技術計算に必要な拡張モジュール類が一気にダウンロード、インストールできる。ただし、Python(x,y)は現在のところ、32bit版しか存在しない。下に説明するEPDを利用すると、64bit版も手に入る。 MacOSXには最初からPythonがインストールされているようである。また、最近のバージョンではNumpy(数値計算用拡張モジュール)も最初からインストールされている。 しかし、その他の科学技術計算用拡張モジュールを一気にインストールするならば、 [[Enthought Python Distribution (EPD)>http://www.enthought.com/products/epd.php]]をインストールするのが一番手っ取り早い。EPDは基本的に有料のディストリビューションだが、アカデミック利用は無料である。 しかし、その他の科学技術計算用拡張モジュールを一気にインストールするならば、EPDをインストールするのが一番手っ取り早い。 *** [[Enthought Python Distribution (EPD)>http://www.enthought.com/products/epd.php]] [#jd71399a] Enthoughtという会社が開発している、科学技術計算用のPython環境。 EPDは基本的に有料のディストリビューションだが、アカデミック利用は無料である。 [[ダウンロードページ>http://www.enthought.com/products/getepd.php]]の最下部にアカデミック版へのリンクがあるので、先に進むと、アカデミックメールアドレスの入力を求められる。ここに、ac.jpで終わるアドレスを入れれば、アカデミック版のダウンロードリンクが送られてくる。 *** 追加モジュール [#x29ac246] Pythonは汎用言語であり、本体には基本的なプログラミング言語としての機能しかない。~ 様々な応用的機能は、モジュールとして提供されており、標準ライブラリ(本体に付属)からサードパーティによるモジュールまで多岐に渡っている。~ ここで解説されている、ブラウン運動の解析を行うためには、以下のモジュールが必要である。 - [[Numpy>http://numpy.scipy.org/]]: 数値計算ライブラリ - [[matplotlib>http://matplotlib.sourceforge.net/]]: プロットライブラリ - [[Scipy>http://www.scipy.org/]]: 科学技術計算ライブラリ これらのモジュールは、個別にダウンロードしてインストールすることも可能だが、Linuxの場合Distributionのパッケージ管理システムからインストール可能な場合が多い。~ EPDを利用する場合、これらは自動的にインストールされるので、なにもする必要が無い。 この他、データ解析に必須ではないが、Pythonを対話的に利用する際非常に便利なのが、[[IPython>http://ipython.scipy.org/moin/]]である。EPDでは自動的にインストールされる。Linuxの場合、Distributionのパッケージシステムに含まれているはずである。 ** Pythonの基本 [#s788f0a2] Python自体の入門には、以下のようなサイトが役に立つ。 - [[Python Documentation>http://www.python.org/doc/]]: Pythonの公式ドキュメント。ここのTutorialは入門に適している。基本は英語。[[日本語訳>http://www.python.jp/Zope/links/python_documents]]もあり。 - [[Dive into Python>http://diveintopython.org/]]: 非常に良く書けたpython入門。英語 - [[インスタント Python>http://www.python.jp/Zope/intro/instant_python_jp]]: 手短なPython入門。リンク先は日本語訳。 科学技術関連のモジュールについては、日本語の解説はほとんど無い。~ 英語だが、Scipyの[[チュートリアル集>http://www.scipy.org/Topical_Software#head-bdb651f088978a8cc75acf5fa2e642713e92ceb8]]のページは多くのドキュメントへのリンクを集めている。特に、[[Additional Documentation>http://www.scipy.org/Wiki/Additional_Documentation?action=show]]のページには、numpy, scipyの詳細な説明が載っている。 Matplotlibに関しては、http://matplotlib.sourceforge.net/contents.html に詳細な解説が載っている。 その他、このWiki上にもPythonを用いた科学技術計算の[[解説>https://granite.phys.s.u-tokyo.ac.jp/wiki/Lab/index.php?PythonForScience]]がある。 *** モジュールと名前空間 [#l0ab4430] Pythonを使う上で一つ注意しなければならないのは、名前空間の問題である。 Pythonでは、次のようにモジュールを読み込む。 import numpy この記述によって、numpyモジュールが読み込まれた。numpyには様々な関数が定義されているが、これらの関数は、numpy.関数名 という形式で呼び出さなければならない。例えば、numpyのexp()関数を呼び出すには、 numpy.exp(3) のように呼ばなければならない。これは、様々なモジュールを読み込んだ際に、モジュール内で定義されている関数名がぶつからないように、名前空間を分離しているからである。いちいちnumpyと打つのが面倒な場合、 from numpy import * とすることができる。この場合、numpyの名前空間がトップレベルの名前空間と統合されるため、exp()と記述するだけで、numpy.exp()を呼ぶことができる。しかし、numpyで定義されている関数とトップレベル名前空間にあるオブジェクトの名前が衝突する場合、問題が起こることもあるので、注意が必要である。妥協案として、 import numpy as np のように、numpyに短い別名を付けることもある。この場合、np.exp()とすれば、numpy.exp()を呼び出すことができる。このページの解説では、この方法を採っている。 *** ヘルプ [#t54b0ff8] コンソールにipythonを使っている場合、 np.exp? というように、オブジェクト名(関数名)の後に?を付けてEnterを押すと、このオブジェクトに関するヘルプが表示される。 *** Matlab, Scilabとの相違点 [#qd64e935] PythonはMatlab系言語とはかなり違った文法を持つが、以下の点は特に注意が必要である。 - リストのIndexは0から始まる。 >>>a=[1,2,3,4] >>>a[0] 1 添字でリスト一部を読み出す際の規則(PythonではSlicingと呼ぶ)も若干違う。 >>>a[1:3] [2,3] これについては、Python Tutorialの[[Stringsに関する部分>http://docs.python.org/tutorial/introduction.html#strings]]にあるように、Indexはリスト要素の境界線を指すと考えると分かりやすい。 - べき乗は**で表す。 >>>2**3 8 ** モジュールの読み込み [#j6f3067a] 以下の解説で利用されるモジュールをあらかじめ読み込んでおく。 import numpy as np import matplotlib.pyplot as plt import fileinput import scipy.optimize as opt ** データの読み込み [#bda76670] ASCIIファイルからデータを読み込むには次のようにする。 d=np.loadtxt('brown.dat',delimiter=',') t=d[:,0] v=d[:,1] plt.plot(t,v) *** Matlab風dlmread() [#nccd4a21] 上述の方法とは別に、ここではMatlab風のdlmread()を定義してみる。 def dlmread(filename, delimiter=None,comment="#"): """Load contents of an ascii file""" lineCount=0 for line in fileinput.input(filename): if line.strip()[0] == comment: continue #Skip lines starting from comment if lineCount == 0: data=np.transpose(np.array([map(float,line.split(delimiter))])) lineCount=lineCount+1; else: data=np.hstack((data, np.transpose(np.array([map(float,line.split(delimiter))])))) lineCount=lineCount+1; return data この関数を使って、brown.datからセーブされた時系列データを読み込む。 d=dlmread('brown.dat',',') t=d[0] #時間 v=d[1] #電圧 plt.plot(t,v) ** 平均、分散、自己相関関数 [#p9c2d81f] 平均と分散は以下のコマンドで求まる。 np.mean(v) np.var(v) 自己相関関数は、np.correlate()関数を使ってもよいが、以下のようにプログラム的に求めても良い。 dt=t[1]-t[0] numPoints=100 a=np.zeros(numPoints) tau=np.zeros(numPoints) for ii in range(numPoints): tau[ii]=dt*ii; a[ii]=np.sum(v[ii:]*v[:np.size(v)-ii])/np.size(v[ii:]) plt.plot(tau,a) ** ヒストグラム [#pd6f4f30] vのヒストグラムを求めるには、 (n,bin,pat)=plt.hist(v,20) x=(bin[:-1]+bin[1:])/2 plt.plot(x,n) のようにすれば良い。~ binには各binの境界値が入っているので、x=...の行でbinの中央値を求めている。 ** Fitting [#d62fd105] 得られたヒストグラムをGaussianでFitする。~ まず、Fitする関数型を定義する。(lambda形式については[[ここ>http://www.diveintopython.org/power_of_introspection/lambda_functions.html]]を参照) fitFunc = lambda p,x: p[0]*np.exp(-(x-p[1])**2/p[2]**2) Fitting関数とデータの差を求める関数。 errFunc = lambda p,x,y: fitFunc(p, x) - y 初期パラメータ p0=[2500,0,2.5] 最適パラメータの探索 (p,success) = opt.leastsq(errFunc, p0[:], args=(x,n)) 得られた結果をチェックする。 nfit = fitFunc(p,x) plt.plot(x,n,x,nfit)
