* Matlabを用いたブラウン運動データ解析 [#gfdde69e]
このページでは、Matlabを用いてブラウン運動実験で取得したデータを解析、プロットする方法を紹介する。ただし、データ解析自体がブラウン運動実験の課題であるため、わざと不完全な解説に留めてある。各自、Matlabのマニュアルなどを参照して、不明部分を補うこと。
#contents
** Matlabの基本 [#he7fa615]
Matlabではベクトルや行列が基本オブジェクトとなる。例えば、
>>a=[1,2,3,4]
a=
1 2 3 4
と入力すれば、aという変数に(1,2,3,4)という数列が代入される。これは、行ベクトルである。なお、>>はMatlabのプロンプトを表す。
列ベクトルを作るには以下のようにする。
>>b=[1;2;3;4]
b=
1
2
3
4
行列は
>> c=[1,2,3,4; 5,6,7,8]
c =
1 2 3 4
5 6 7 8
のように作る。
aの要素を参照するには
>>a(1)
ans = 1
>>a(2:3)
ans = 2 3
のようにする。
*** ヘルプ [#f28ec259]
Matlabでは、
>>doc command_name
と打ち込むと、command_nameに関するヘルプが表示される。以下の例において不明なコマンドについては、docを使ってヘルプを参照すれば良い。
*** スクリプトファイルとセル [#p59095c4]
Matlabでは、上記のようにコマンドプロンプトからコマンドを一行一行打ち込むことも可能であるが、複雑な処理を何度も繰り返す場合、これでは効率が悪い。コマンドはファイルに一括して記述し、それをまとめて実行することが可能である。このようなファイルをMatlab scriptと呼ぶ。Matlabスクリプトはさらに、セルという実行単位に分割することが可能である。例えば、
%% Cell 1
a=1;
b=2*a;
%% Cell 2
c=sin(3);
というようなファイルがある場合、行頭の%%から、次のまでが一つのセルとなる。ちなみに、%はコメントを表すので、"Cell 1"等の文字列は無視される。
セルを用いると、各セルだけを実行させることが可能になる。Matlabのエディタ上で、セルの中にカーソルを置いて、Ctrl+Enterを押すと、そのセルだけがMatlabインタプリタに送られ、実行される。これは、スクリプトの一部だけを実行したい場合に便利である。
** データのインポート [#n9df54bc]
まずは、解析すべきデータをMatlabに読み込まなければならない。ブラウン運動でセーブしたデータは、カンマ区切りのテキストファイルになっている。データファイルのファイル名がdatafile.datであった場合、
>>d=dlmread('datafile.dat');
で、データファイルの中身が変数dに読み込まれる。datafile.datが抵抗の熱雑音を測った結果だとすると、ファイルの中身は2列であり、1列目が時間、2列目が電圧である。変数dもN行x2列の行列で、Nは測定されたデータ点数である。
時間と電圧のデータを取り出す。
>>t=d(:,1); %dの一列目を取り出す(時間データ)
>>v=d(:,2); %dの２列目を取り出す(電圧データ)
データのプロット
>>plot(t,v);
** 平均、分散、自己相関関数 [#w4fd2fc2]
平均と分散はMatlabのmean()とvar()コマンドを使えば良い。
>>mean(v) %平均
>>var(v) %分散
自己相関関数は例えば次のようなスクリプトで計算できる
%% Auto correlation
c=[]; tau=[];
dt=t(2)-t(1); %サンプリング間隔
for n=0:100,
tau(n+1)=dt*n;
c(n+1)=sum(v(1:end-n).*v(n+1:end))/(length(v)-n);
end
plot(tau,c);
** ヒストグラムとフィット [#qa4a4ddf]
vのヒストグラムは、
>>hist(v,30)
でプロットできる。ここで、30はBin数である。vのサンプル数に応じて変えてみると良い。
ヒストグラムのデータは、
>>[n,x]=hist(v,30);
とすることで、nとxに代入される。xは各Binの中央値、nは度数である。
>>plot(x,n);
を実行すればnとxの意味がわかるはずである。
*** データのフィッティング [#hf57445f]
熱雑音の電圧は、一般にガウス分布にしたがう。よって、得られたヒストグラム[n,x]をガウス分布でフィッテイングしたくなる。~
Matlabにおいて、フィッティングを行うにはいくつかの方法がある。線形フィッティングの場合、polyfit()やpolyval()を使えば良い。しかし、ガウス分布は非線形関数である。この場合、Statistics Toolboxのnlinfit()関数や、Curve Fitting Toolboxを使うのが一般的であるが、ここではこうしたToolboxを必要としない方法を紹介する。
Matlabにおいて、フィッティングを行うにはいくつかの方法がある。線形フィッティングの場合、polyfit()やpolyval()を使えば良い。しかし、ガウス分布は非線形関数である。この場合、Statistics Toolboxのnlinfit()関数や、Curve Fitting Toolboxを使うのが一般的であるが、ここではこうしたToolboxを必要としない方法を紹介する。以下がサンプルプログラムである。
%% Histogram
[n,x]=hist(v,30);
p0=[2500,0,2.5]; %Initial parameter
n1=p0(1)*exp(-(x-p0(2)).^2/p0(3)^2); %Function values with the initial parameters
%% Test plot
plot(x,n,x,n1);
%% Fitting
testfun=@(p)sum((n-p(1)*exp(-(x-p(2)).^2/p(3)^2)).^2);
popt=fminsearch(testfun,p0);
%% Confirm results
n2=popt(1)*exp(-(x-popt(2)).^2/popt(3)^2);
%% Confirm results
plot(x,n,x,n2)
ガウス分布は、全体のスケール、平均値、標準偏差、という三つのパラメータで特徴付けられる。p0という変数にはこの三つのパラメータの初期値がベクトルの形で代入されている。n1はp0をパラメータとしたときのガウス分布である。
ガウス分布は、全体のスケール、平均値、標準偏差、という三つのパラメータで特徴付けられる。p0という変数にはこの三つのパラメータの初期値がベクトルの形で代入されている。n1はp0をパラメータとしたときのガウス分布を計算したものである。
testfun=@(p)sum((n-p(1)*exp(-(x-p(2)).^2/p(3)^2)).^2);
という部分は、最小化したい関数を定義している。@(p)...という形式は、anonymous functionと呼ばれて、簡単に関数を定義するための書式である。
という部分は、最小化したい関数を定義している。@(p)...という形式は、anonymous functionと呼ばれて、簡単に関数を定義するための書式である。~
この関数では、xの各点で実測データnと、パラメータpを用いて計算されたガウス分布の差を取って、二乗した後に全点の和を取っている。この関数を最小化するようなpを見つけることが我々の目的である。
popt=fminsearch(testfun,p0);
この部分では、Matlabのfiminsearch関数を用いてtestfun(p)を最小化するようなpを、p0を初期値として求めている。これがFittingの本体である。~
あとは、求められた最適パラメータpoptの検証を行っているだけである。
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
