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

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

↑

インストール †

Python自体は、PythonのWebsiteや日本Pythonユーザー会などからダウンロードできる。
Linuxの場合は、Distributionのパッケージシステムからインストールするのが楽で確実である。最初からインストール済みのDistributionも多い。 WindowsやMacの場合、Enthought Python Distribution (EPD)をインストールするのが一番手っ取り早い。EPDには、科学技術計算に必要なPythonモジュールが一通りインストールされている。

↑

追加モジュール †

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

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

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

↑

Pythonの基本 †

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

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

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

↑

モジュールと名前空間 †

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
↑

データの読み込み †

pythonにはMatlabのdlmread()のような便利な関数がないので、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