Pythonによるブラウン運動データ解析

このページではPythonを用いてブラウン運動実験で得られたデータを解析する方法について解説をする。データ解析も実験課題の一部なので、わざと不完全な解説になっている。不明な点は、各種参考文献やWeb siteを参照して自習すること。

インストール

一般に、Pythonは、PythonのWebsite日本Pythonユーザー会などからダウンロードできる。 以下に、各種OSにおけるインストール方法を概説する。

Linuxの場合は、Distributionのパッケージシステムからインストールするのが楽で確実である。最初からインストール済みのDistributionも多い。 また、下で説明するEnthought Python Distribution (EPD)を利用する方法もある。 この場合、科学技術計算環境が一気に揃うので便利である。

Windowsの場合、Python(x,y)というプロジェクトから、Pythonと科学技術計算に必要な拡張モジュール類が一気にダウンロード、インストールできる。ただし、Python(x,y)は現在のところ、32bit版しか存在しない。下に説明するEPDを利用すると、64bit版も手に入る。

MacOSXには最初からPythonがインストールされているようである。また、最近のバージョンではNumpy(数値計算用拡張モジュール)も最初からインストールされている。 しかし、その他の科学技術計算用拡張モジュールを一気にインストールするならば、EPDをインストールするのが一番手っ取り早い。

Enthought Python Distribution (EPD)

Enthoughtという会社が開発している、科学技術計算用のPython環境。 EPDは基本的に有料のディストリビューションだが、アカデミック利用は無料である。 ダウンロードページの最下部にアカデミック版へのリンクがあるので、先に進むと、アカデミックメールアドレスの入力を求められる。ここに、ac.jpで終わるアドレスを入れれば、アカデミック版のダウンロードリンクが送られてくる。

追加モジュール

Pythonは汎用言語であり、本体には基本的なプログラミング言語としての機能しかない。
様々な応用的機能は、モジュールとして提供されており、標準ライブラリ(本体に付属)からサードパーティによるモジュールまで多岐に渡っている。
ここで解説されている、ブラウン運動の解析を行うためには、以下のモジュールが必要である。

これらのモジュールは、個別にダウンロードしてインストールすることも可能だが、Linuxの場合Distributionのパッケージ管理システムからインストール可能な場合が多い。
EPDを利用する場合、これらは自動的にインストールされるので、なにもする必要が無い。

この他、データ解析に必須ではないが、Pythonを対話的に利用する際非常に便利なのが、IPythonである。EPDでは自動的にインストールされる。Linuxの場合、Distributionのパッケージシステムに含まれているはずである。

Pythonの基本

Python自体の入門には、以下のようなサイトが役に立つ。

科学技術関連のモジュールについては、日本語の解説はほとんど無い。
英語だが、Scipyのチュートリアル集のページは多くのドキュメントへのリンクを集めている。特に、Additional Documentationのページには、numpy, scipyの詳細な説明が載っている。

Matplotlibに関しては、http://matplotlib.sourceforge.net/contents.html に詳細な解説が載っている。

その他、このWiki上にもPythonを用いた科学技術計算の解説がある。

モジュールと名前空間

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()を呼び出すことができる。このページの解説では、この方法を採っている。

ヘルプ

コンソールにipythonを使っている場合、

np.exp?

というように、オブジェクト名(関数名)の後に?を付けてEnterを押すと、このオブジェクトに関するヘルプが表示される。

Matlab, Scilabとの相違点

PythonはMatlab系言語とはかなり違った文法を持つが、以下の点は特に注意が必要である。

添字でリスト一部を読み出す際の規則(PythonではSlicingと呼ぶ)も若干違う。

>>>a[1:3]
[2,3]

これについては、Python TutorialのStringsに関する部分にあるように、Indexはリスト要素の境界線を指すと考えると分かりやすい。

モジュールの読み込み

以下の解説で利用されるモジュールをあらかじめ読み込んでおく。

import numpy as np
import matplotlib.pyplot as plt 
import fileinput
import scipy.optimize as opt

データの読み込み

ASCIIファイルからデータを読み込むには次のようにする。

d=np.loadtxt('brown.dat',delimiter=',')
t=d[:,0]
v=d[:,1]
plt.plot(t,v)

Matlab風dlmread()

上述の方法とは別に、ここでは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)

平均、分散、自己相関関数

平均と分散は以下のコマンドで求まる。

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)

ヒストグラム

vのヒストグラムを求めるには、

(n,bin,pat)=plt.hist(v,20)
x=(bin[:-1]+bin[1:])/2
plt.plot(x,n)

のようにすれば良い。
binには各binの境界値が入っているので、x=...の行でbinの中央値を求めている。

Fitting

得られたヒストグラムをGaussianでFitする。
まず、Fitする関数型を定義する。(lambda形式についてはここを参照)

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)

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2013-07-09 (火) 14:50:24